Inversible STFT et ISTFT en Python

Question

Y avec inverse correspondant toute forme d'usage général de Fourier à court terme de transformation transform construit en SciPy ou NumPy ou quoi?

Il y a la fonction specgram de pyplot 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.

  

Ceci est une fonction d'assistance qui met en œuvre les points communs entre la   204 #psd, cdd et spectrogramme. Il est    pas destiné à être utilisé en dehors de mLab

Je ne sais pas si cela peut être utilisé pour faire une STFT et ISTFT, cependant. Y at-il quelque chose d'autre, ou devrais-je traduire quelque chose comme ces fonctions Matlab?

Je sais comment écrire ma propre mise en œuvre ad hoc; Je suis en train de chercher quelque chose complet, qui peut gérer différentes fonctions de fenêtrage (mais a un défaut sain d'esprit), est entièrement inversible avec des fenêtres de vie chère (de istft(stft(x))==x), testé par plusieurs personnes, pas hors par une erreur, gère les ends et bien zéro padding, rapide mise en œuvre de RFFT pour l'entrée réelle, etc.

Answer 1

Je suis un peu en retard à ce sujet, mais scipy réalisé a intégré function istft à partir de 0.19.0

Answer 2

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. la compréhension de la liste est un petit truc que j'aime utiliser pour simuler le traitement du bloc de signaux dans numpy / scipy. Il est comme blkproc dans Matlab. Au lieu d'une boucle de for, j'applique une commande (par exemple, fft) pour chaque trame du signal à l'intérieur d'une compréhension de la liste, puis scipy.array elle met à un tableau 2D. J'utilise ceci pour faire spectrogrammes, chromagrams, MFCC-grammes, et bien plus encore.
2. Pour cet exemple, j'utilise une méthode de chevauchement et d'addition naïve dans istft. Afin de reconstituer le signal d'origine à la somme des fonctions de fenêtre successives doit être constante, de préférence égale à l'unité (1,0). Dans ce cas, j'ai choisi la fenêtre Hann (ou hanning) et un recouvrement de 50%, ce qui fonctionne parfaitement. Voir cette discussion pour plus d'informations.
3. Il y a probablement des moyens plus fondées sur des principes de calcul du ISTFT. Cet exemple est principalement destiné à être é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 3

Voici le code STFT que j'utilise. STFT + ISTFT donne ici reconstruction parfaite (même pour les premières images). Je légèrement modifié le code donné ici par Steve Tjoa. Ici, l'amplitude 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 4

librosa.core.stft et

trouvé un autre STFT, mais aucune fonction inverse correspondant:

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 de fenêtre en tant que tableau
• L est le chevauchement, dans des échantillons

Aucune des réponses ci-dessus a bien fonctionné OOTB pour moi. Donc, je Steve modifié Tjoa de.

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...

J'ai aussi trouvé cela sur GitHub, mais il semble fonctionner sur les 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))

Je pense que scipy.signal a ce que vous recherchez. Il a par défaut raisonnables, prend en charge plusieurs types de fenêtres, etc ...

http: // docs. scipy.org/doc/scipy-0.17.0/reference/generated/scipy.signal.spectrogram.html

from scipy.signal import spectrogram
freq, time, Spec = spectrogram(signal)

Une version corrigée de la réponse basj.

import scipy, numpy as np
import matplotlib.pyplot as plt

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.vstack([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]
    rcs=int(np.ceil(float(X.shape[0])/float(overlap)))*fftsize
    print(rcs)
    x=np.zeros(rcs)
    wsum=np.zeros(rcs)
    for n,i in zip(X,range(0,len(X)*hop,hop)): 
        l=len(x[i:i+fftsize])
        x[i:i+fftsize] += np.fft.irfft(n).real[:l]   # overlap-add
        wsum[i:i+fftsize] += w[:l]
    pos = wsum != 0
    x[pos] /= wsum[pos]
    return x

a=np.random.random((65536))
b=istft(stft(a))
plt.plot(range(len(a)),a,range(len(b)),b)
plt.show()

Si vous avez accès à une bibliothèque binaire C qui fait ce que vous voulez, puis utilisez http: // code.google.com/p/ctypesgen/ pour générer une interface de python de cette bibliothèque.

Answer 5

Answer 6

Answer 7

Answer 8

trouvé un autre STFT, mais aucune fonction inverse correspondant:

Answer 9

Answer 10

Aucune des réponses ci-dessus a bien fonctionné OOTB pour moi. Donc, je Steve modifié Tjoa de.

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 11

J'ai aussi trouvé cela sur GitHub, mais il semble fonctionner sur les 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))

Answer 12

Je pense que scipy.signal a ce que vous recherchez. Il a par défaut raisonnables, prend en charge plusieurs types de fenêtres, etc ...

http: // docs. scipy.org/doc/scipy-0.17.0/reference/generated/scipy.signal.spectrogram.html

from scipy.signal import spectrogram
freq, time, Spec = spectrogram(signal)

Answer 13

Une version corrigée de la réponse basj.

import scipy, numpy as np
import matplotlib.pyplot as plt

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.vstack([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]
    rcs=int(np.ceil(float(X.shape[0])/float(overlap)))*fftsize
    print(rcs)
    x=np.zeros(rcs)
    wsum=np.zeros(rcs)
    for n,i in zip(X,range(0,len(X)*hop,hop)): 
        l=len(x[i:i+fftsize])
        x[i:i+fftsize] += np.fft.irfft(n).real[:l]   # overlap-add
        wsum[i:i+fftsize] += w[:l]
    pos = wsum != 0
    x[pos] /= wsum[pos]
    return x

a=np.random.random((65536))
b=istft(stft(a))
plt.plot(range(len(a)),a,range(len(b)),b)
plt.show()

Answer 14

Si vous avez accès à une bibliothèque binaire C qui fait ce que vous voulez, puis utilisez http: // code.google.com/p/ctypesgen/ pour générer une interface de python de cette bibliothèque.