STFT e ISTFT invertibles en Python
-
20-09-2019 - |
Pregunta
¿Existe alguna forma de propósito general de transformada de Fourier de corto tiempo ¿Con la transformación inversa correspondiente integrada en SciPy o NumPy o lo que sea?
Ahí está el diagrama de bits specgram
función en matplotlib, que llama ax.specgram()
, que llama mlab.specgram()
, que llama _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
Esta es una función auxiliar que implementa la comunidad entre los 204 #PSD, CSD y el espectrograma.Es NO destinado a ser utilizado fuera de mlab
Sin embargo, no estoy seguro de si esto se puede usar para hacer STFT e ISTFT.¿Hay algo más o debería traducir algo como estas funciones de MATLAB?
Sé cómo escribir mi propia implementación ad hoc;Solo estoy buscando algo con todas las funciones, que pueda manejar diferentes funciones de ventanas (pero que tenga un valor predeterminado sensato), que sea completamente invertible con ventanas COLA (istft(stft(x))==x
), probado por varias personas, sin errores uno por uno, maneja bien los extremos y el relleno de ceros, implementación RFFT rápida para entradas reales, etc.
Solución
Estoy un poco tarde para esto, pero se dio cuenta de scipy ha incorporado istft función como de 0.19.0
Otros consejos
Aquí está mi código Python, simplificado para esta respuesta:
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
Notas:
- El lista por comprensión es un pequeño truco que me gusta usar para simular la ejecución del bloque de señales en numpy / scipy. Es como
blkproc
en Matlab. En lugar de un buclefor
, aplico un comando (por ejemplo,fft
) para cada trama de la señal dentro de una lista por comprensión, y luegoscipy.array
arroja a un 2D-array. Lo utilizo para hacer espectrogramas, chromagrams, MFCC-gramos, y mucho más. - En este ejemplo, utilizo un solapamiento y suma ingenua método en el
istft
. Con el fin de reconstruir la señal original la suma de las funciones de la ventana secuenciales debe ser constante, preferiblemente igual a la unidad (1,0). En este caso, he elegido la ventana de Hann (ohanning
) y una superposición del 50%, que funciona a la perfección. Ver esta discusión para más información. - Probablemente hay maneras más de principio del cálculo de la ISTFT. Este ejemplo está destinado principalmente a ser educativo.
Una prueba:
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)')
Aquí está el código STFT que uso.STFT + ISTFT aquí da reconstrucción perfecta (incluso para los primeros fotogramas).Modifiqué ligeramente el código proporcionado aquí por Steve Tjoa:aquí la magnitud de la señal reconstruida es la misma que la de la señal de entrada.
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
y istft
parecen bastante similar a lo que estaba buscando, a pesar de que no existían en el tiempo:
librosa.core.stft(y, n_fft=2048, hop_length=None, win_length=None, window=None, center=True, dtype=<type 'numpy.complex64'>)
No invierten exactamente, sin embargo; los extremos se estrechan.
Encontrados otro STFT, pero no correspondiente función inversa:
http://code.google.com /p/pytfd/source/browse/trunk/pytfd/stft.py
def stft(x, w, L=None):
...
return X_stft
- w es una función de ventana como una matriz
- L es la superposición, en muestras
Ninguna de las respuestas anteriores ha funcionado bien para mí OOTB. Así que he modificado Steve 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...
También encontrado esto en GitHub, pero parece operar en tuberías en lugar de matrices normales:
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))
Creo scipy.signal tiene lo que busca. Tiene configuraciones razonables, es compatible con varios tipos de ventanas, 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)
Una versión fija de la respuesta de 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 usted tiene acceso a una biblioteca binaria C que hace lo que quiere, a continuación, utiliza http: // code.google.com/p/ctypesgen/ para generar una interfaz de Python a esa biblioteca.