STFT e ISTFT invertibili in Python
-
20-09-2019 - |
Domanda
Esiste una forma generica di Trasformata di Fourier di breve periodo con la corrispondente trasformazione inversa incorporata in SciPy o NumPy o altro?
C'è il pyplot specgram
funzione in matplotlib, che chiama ax.specgram()
, che chiama mlab.specgram()
, che chiama _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.
Ma
Questa è una funzione di supporto che implementa la comunanza tra 204 #PSD, CSD e spettrogramma.È NON pensato per essere utilizzato al di fuori di mlab
Tuttavia, non sono sicuro che questo possa essere utilizzato per eseguire STFT e ISTFT.C'è qualcos'altro o dovrei tradurre qualcosa del genere queste funzioni MATLAB?
So come scrivere la mia implementazione ad hoc;Sto solo cercando qualcosa di completo, che possa gestire diverse funzioni di finestra (ma abbia un'impostazione predefinita sana), sia completamente invertibile con finestre COLA (istft(stft(x))==x
), testato da più persone, nessun errore singolo, gestisce bene le estremità e il riempimento zero, implementazione RFFT rapida per input reale, ecc.
Soluzione
Sono un po 'tardi per questo, ma SciPy realizzato ha integrato istft funzione di 0.19.0
Altri suggerimenti
Ecco il mio codice Python, semplificato per questa risposta:
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
Appunti:
- IL comprensione delle liste è un piccolo trucco che mi piace usare per simulare l'elaborazione a blocchi dei segnali in numpy/scipy.È come
blkproc
in Matlab.Invece di afor
loop, applico un comando (ad esempio,fft
) a ciascun frame del segnale all'interno di una lista di comprensione, e poiscipy.array
lo lancia in un array 2D.Lo uso per creare spettrogrammi, cromagrammi, grammi MFCC e molto altro. - Per questo esempio, utilizzo un ingenuo metodo di sovrapposizione e aggiunta
istft
.Per poter ricostruire il segnale originale la somma delle funzioni finestra sequenziali deve essere costante, preferibilmente pari all'unità (1,0).In questo caso, ho scelto l'Hann (orhanning
) finestra e una sovrapposizione del 50% che funziona perfettamente.Vedere questa discussione per maggiori informazioni. - Probabilmente esistono metodi più basati su principi per calcolare l’ISTFT.Questo esempio è principalmente pensato per essere educativo.
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)')
Ecco il codice STFT che uso. STFT + ISTFT qui dà perfetta ricostruzione (anche per i primi fotogrammi). Ho leggermente modificato il codice qui data da Steve Tjoa:. Qui l'ampiezza del segnale ricostruito è la stessa di quella del segnale di ingresso
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
librosa.core.stft
e istft
guardare piuttosto simile a quello che stavo cercando, anche se non esistessero al tempo:
librosa.core.stft(y, n_fft=2048, hop_length=None, win_length=None, window=None, center=True, dtype=<type 'numpy.complex64'>)
Non si invertono esattamente, però; le estremità sono rastremate.
Trovato un'altra STFT, ma nessuna corrispondente funzione inversa:
http://code.google.com /p/pytfd/source/browse/trunk/pytfd/stft.py
def stft(x, w, L=None):
...
return X_stft
- w è una funzione finestra come una matrice
- L è la sovrapposizione, in campioni
Nessuna delle risposte di cui sopra ha funzionato bene OOTB per me. Così ho modificato Steve Tjoa di.
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...
Ho trovato anche questo su GitHub, ma sembra operare su tubazioni, invece di array normali:
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))
Credo scipy.signal ha quello che stai cercando. Ha valori di default ragionevoli, supporta più tipi di finestre, ecc ...
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)
Una versione fissa di risposta di 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()
Se si ha accesso a una libreria C binario che fa quello che si vuole, quindi utilizzare http: // code.google.com/p/ctypesgen/ per generare un'interfaccia Python a quella libreria.