Python中可转换的STFT和ISTFT

是否有任何通用的短时傅立叶变换forms，并在SciPy或NumPy中内置相应的逆变换？

在matplotlib中有pyplot specgram函数，它调用ax.specgram() ，它调用mlab.specgram() ，它调用_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. 

但

这是一个帮助函数，实现204 #psd，csd和频谱图之间的通用性。 这并不意味着在mlab之外使用

不过，我不确定这是否可以用来做STFT和ISTFT。 还有什么，或者我应该翻译像这些MATLAB函数 ？

我知道如何编写我自己的临时实现; 我只是寻找一些function齐全的，可以处理不同的窗口function（但有一个理智的默认），完全可逆COLA窗口（ istft(stft(x))==x ），由多人testing，没有逐个错误，处理结束和零填充，为实际input快速RFFT实现等。

Answer 1

我对这个有点晚了，但是实现了0.19.0的scipy内置的istft函数

Answer 2

这是我的Python代码，简化了这个答案：

 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 

笔记：

1. 列表理解是我喜欢用来模拟numpy / scipy中的信号块处理的小技巧。 这就像在Matlab的blkproc 。 而不是一个for循环，我将一个命令（例如， fft ）应用于列表理解中的每个信号帧，然后scipy.array将其转换为二维数组。 我用它来制作谱图，色谱图，MFCC-grams等等。
2. 对于这个例子，我在istft使用了天真的overlap-and-add方法。 为了重构原始信号，顺序窗口函数的和必须是恒定的，最好等于1（1.0）。 在这种情况下，我select了Hann（或汉hanning ）窗口和50％的重叠，完美的作品。 有关更多信息，请参阅此讨论 。
3. 计算ISTFT可能有更多原则性的方法。 这个例子主要是为了教育。

一个testing：

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

这是我使用的STFT代码。 STFT + ISTFT在这里给予完美的重build （即使是第一帧）。 我稍微修改了Steve Tjoa给出的代码：这里重build信号的幅度与input信号的幅度相同。

 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和istft外观与我所寻找的非常相似，尽pipe它们当时并不存在：

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

尽pipe如此，它们并不完全相反。 两端是锥形的。

Answer 5

发现另一个STFT，但没有相应的反函数：

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

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

• w是一个数组的窗口函数
• L是样本中的重叠

Answer 6

上述两个答案对我来说都不是很好。 所以我修改了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

我也发现这在GitHub上，但它似乎在pipe道而不是正常的数组上运行：

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 8

我认为scipy.signal有你在找什么。 它有合理的默认值，支持多种窗口types等…

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 9

如果您有权访问C二进制库，然后使用http://code.google.com/p/ctypesgen/为该库生成一个Python接口。;