fourier - Invertible STFT e ISTFT en Python

numpy fft python example

Question

¿Hay alguna forma de propósito general de transformada de Fourier a corto plazo con la correspondiente transformación inversa incorporada en SciPy o NumPy o lo que sea?

Existe la función specgram de specgram en matplotlib, que llama ax.specgram() , que llama a mlab.specgram() , que llama a _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 los puntos en común entre 204 #psd, csd y el espectrograma. NO está destinado a ser utilizado fuera de mlab

Sin embargo, no estoy seguro de si esto se puede usar para hacer un 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 tiene un buen valor predeterminado), es totalmente invertible con ventanas COLA ( istft(stft(x))==x ), probado por varias personas, no errores off-by-one, maneja bien los extremos y el relleno cero, rápida implementación de RFFT para entradas reales, etc.

Answer 1

Aquí está el código STFT que utilizo. STFT + ISTFT aquí ofrece una reconstrucción perfecta (incluso para los primeros cuadros). Modifiqué ligeramente el código dado 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

Answer 2

Aquí está mi código de 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:

1. La lista de comprensión es un pequeño truco que me gusta usar para simular el procesamiento de bloque de señales en números y caracteres. Es como blkproc en Matlab. En lugar de un bucle for , aplico un comando (por ejemplo, fft ) a cada cuadro de la señal dentro de una lista de comprensión, y luego scipy.array lo scipy.array en un arreglo 2D. Uso esto para hacer espectrogramas, cromagramas, MFCC-gramos y mucho más.
2. Para este ejemplo, uso un método ingenuo de superposición y adición en istft . Para reconstruir la señal original, la suma de las funciones de la ventana secuencial debe ser constante, preferiblemente igual a la unidad (1.0). En este caso, he elegido la ventana Hann (o hanning ) y una superposición del 50% que funciona perfectamente. Vea esta discusión para más información.
3. Probablemente hay más formas de computar el 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)'')

Answer 3

Creo que scipy.signal tiene lo que estás buscando. Tiene valores predeterminados razonables, admite múltiples 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)

Answer 4

Encontré otro STFT, pero ninguna función inversa correspondiente:

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.

Answer 5

Llego un poco tarde a esto, pero me di cuenta de que scipy tiene istft función de istft incorporada istft 0.19.0

Answer 6

Ninguna de las respuestas anteriores funcionó bien OOTB para mí. Así que modifiqué la 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 7

Si tiene acceso a una biblioteca binaria de C que hace lo que quiere, entonces use http://code.google.com/p/ctypesgen/ para generar una interfaz Python para esa biblioteca.

Answer 8

También encontré esto en GitHub, pero parece funcionar 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))

Answer 9

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()

Answer 10

librosa.core.stft y istft parecen bastante a lo que estaba buscando, aunque no existían en ese momento:

librosa.core.stft(y, n_fft=2048, hop_length=None, win_length=None, window=None, center=True, dtype=<type ''numpy.complex64''>)

Sin embargo, no se invierten exactamente; Los extremos son cónicos.