STFT et ISTFT inversibles en Python

Question

Existe-t-il une forme polyvalente de transformée de Fourier à court terme avec la transformation inverse correspondante intégrée dans SciPy ou NumPy ou autre ?

Il y a le pyplot specgram dans matplotlib, qui appelle ax.specgram() qui appelle mlab.specgram() qui appelle _spectral_helper() :

#The checks for if y is x are so that we can use the same function to
#implement the core of psd(), csd(), and spectrogram() without doing
#extra calculations.  We return the unaveraged Pxy, freqs, and t.

pero

Il s'agit d'une fonction d'aide qui implémente le point commun entre le 204 #psd, csd, et spectrogramme. Elle est PAS destiné à être utilisé en dehors de mlab

Je ne suis pas sûr que cela puisse être utilisé pour faire une STFT et une ISTFT, cependant. Y a-t-il autre chose, ou dois-je traduire quelque chose comme ces fonctions MATLAB ?

Je sais comment écrire ma propre implémentation ad-hoc ; je recherche simplement quelque chose de complet, qui puisse gérer différentes fonctions de fenêtrage (mais qui ait une valeur par défaut saine), qui soit totalement inversible avec COLA Windows ( istft(stft(x))==x ), testé par plusieurs personnes, pas d'erreurs de type "off-by-one", gère bien les extrémités et le remplissage du zéro, implémentation RFFT rapide pour une entrée réelle, etc.

Answer 1

Voici mon code Python, simplifié pour cette réponse :

import scipy, pylab

def stft(x, fs, framesz, hop):
    framesamp = int(framesz*fs)
    hopsamp = int(hop*fs)
    w = scipy.hanning(framesamp)
    X = scipy.array([scipy.fft(w*x[i:i+framesamp]) 
                     for i in range(0, len(x)-framesamp, hopsamp)])
    return X

def istft(X, fs, T, hop):
    x = scipy.zeros(T*fs)
    framesamp = X.shape[1]
    hopsamp = int(hop*fs)
    for n,i in enumerate(range(0, len(x)-framesamp, hopsamp)):
        x[i:i+framesamp] += scipy.real(scipy.ifft(X[n]))
    return x

Notes :

1. El compréhension de la liste est une petite astuce que j'aime utiliser pour simuler le traitement en bloc de signaux dans numpy/scipy. C'est comme blkproc dans Matlab. Au lieu d'un for boucle, j'applique une commande (par exemple, fft ) à chaque trame du signal dans une liste de compréhension, puis scipy.array le transforme en un tableau 2D. Je l'utilise pour faire des spectrogrammes, des chromagrammes, des MFCC-grammes, et bien plus encore.
2. Pour cet exemple, j'utilise une méthode naïve de chevauchement et d'ajout en istft . Afin de reconstruire le signal original, la somme des fonctions de la fenêtre séquentielle doit être constante, de préférence égale à l'unité (1,0). Dans ce cas, j'ai choisi la fonction de Hann (ou hanning ) et un chevauchement de 50% qui fonctionne parfaitement. Voir cette discussion pour plus d'informations.
3. Il existe probablement des méthodes de calcul de l'ISTFT plus respectueuses des principes. Cet exemple se veut principalement éducatif.

Un test :

if __name__ == '__main__':
    f0 = 440         # Compute the STFT of a 440 Hz sinusoid
    fs = 8000        # sampled at 8 kHz
    T = 5            # lasting 5 seconds
    framesz = 0.050  # with a frame size of 50 milliseconds
    hop = 0.025      # and hop size of 25 milliseconds.

    # Create test signal and STFT.
    t = scipy.linspace(0, T, T*fs, endpoint=False)
    x = scipy.sin(2*scipy.pi*f0*t)
    X = stft(x, fs, framesz, hop)

    # Plot the magnitude spectrogram.
    pylab.figure()
    pylab.imshow(scipy.absolute(X.T), origin='lower', aspect='auto',
                 interpolation='nearest')
    pylab.xlabel('Time')
    pylab.ylabel('Frequency')
    pylab.show()

    # Compute the ISTFT.
    xhat = istft(X, fs, T, hop)

    # Plot the input and output signals over 0.1 seconds.
    T1 = int(0.1*fs)

    pylab.figure()
    pylab.plot(t[:T1], x[:T1], t[:T1], xhat[:T1])
    pylab.xlabel('Time (seconds)')

    pylab.figure()
    pylab.plot(t[-T1:], x[-T1:], t[-T1:], xhat[-T1:])
    pylab.xlabel('Time (seconds)')

Answer 2

Voici le code STFT que j'utilise. STFT + ISTFT donne ici reconstruction parfaite (même pour les premières images). J'ai légèrement modifié le code donné ici par Steve Tjoa : ici la magnitude du signal reconstruit est la même que celle du signal d'entrée.

import scipy, numpy as np

def stft(x, fftsize=1024, overlap=4):   
    hop = fftsize / overlap
    w = scipy.hanning(fftsize+1)[:-1]      # better reconstruction with this trick +1)[:-1]  
    return np.array([np.fft.rfft(w*x[i:i+fftsize]) for i in range(0, len(x)-fftsize, hop)])

def istft(X, overlap=4):   
    fftsize=(X.shape[1]-1)*2
    hop = fftsize / overlap
    w = scipy.hanning(fftsize+1)[:-1]
    x = scipy.zeros(X.shape[0]*hop)
    wsum = scipy.zeros(X.shape[0]*hop) 
    for n,i in enumerate(range(0, len(x)-fftsize, hop)): 
        x[i:i+fftsize] += scipy.real(np.fft.irfft(X[n])) * w   # overlap-add
        wsum[i:i+fftsize] += w ** 2.
    pos = wsum != 0
    x[pos] /= wsum[pos]
    return x

Answer 3

J'ai trouvé une autre STFT, mais pas de fonction inverse correspondante :

http://code.google.com/p/pytfd/source/browse/trunk/pytfd/stft.py

def stft(x, w, L=None):
    ...
    return X_stft

• w est une fonction fenêtre sous forme de tableau
• L est le chevauchement, en échantillons

Answer 4

Aucune des réponses ci-dessus n'a bien fonctionné pour moi. J'ai donc modifié celle de Steve Tjoa.

import scipy, pylab
import numpy as np

def stft(x, fs, framesz, hop):
    """
     x - signal
     fs - sample rate
     framesz - frame size
     hop - hop size (frame size = overlap + hop size)
    """
    framesamp = int(framesz*fs)
    hopsamp = int(hop*fs)
    w = scipy.hamming(framesamp)
    X = scipy.array([scipy.fft(w*x[i:i+framesamp]) 
                     for i in range(0, len(x)-framesamp, hopsamp)])
    return X

def istft(X, fs, T, hop):
    """ T - signal length """
    length = T*fs
    x = scipy.zeros(T*fs)
    framesamp = X.shape[1]
    hopsamp = int(hop*fs)
    for n,i in enumerate(range(0, len(x)-framesamp, hopsamp)):
        x[i:i+framesamp] += scipy.real(scipy.ifft(X[n]))
    # calculate the inverse envelope to scale results at the ends.
    env = scipy.zeros(T*fs)
    w = scipy.hamming(framesamp)
    for i in range(0, len(x)-framesamp, hopsamp):
        env[i:i+framesamp] += w
    env[-(length%hopsamp):] += w[-(length%hopsamp):]
    env = np.maximum(env, .01)
    return x/env # right side is still a little messed up...

Answer 5

J'ai également trouvé ceci sur GitHub, mais il semble fonctionner sur des pipelines au lieu de tableaux normaux :

http://github.com/ronw/frontend/blob/master/basic.py#LID281

def STFT(nfft, nwin=None, nhop=None, winfun=np.hanning):
    ...
    return dataprocessor.Pipeline(Framer(nwin, nhop), Window(winfun),
                                  RFFT(nfft))

def ISTFT(nfft, nwin=None, nhop=None, winfun=np.hanning):
    ...
    return dataprocessor.Pipeline(IRFFT(nfft), Window(winfun),
                                  OverlapAdd(nwin, nhop))