python - Python中的可逆STFT和ISTFT

Question

是否有short-time Fourier transform 的通用形式？将相应的逆变换内置到 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.

但是

This is a helper function that implements the commonality between the 204 #psd, csd, and spectrogram. It is NOT meant to be used outside of mlab

不过，我不确定这是否可用于进行 STFT 和 ISTFT。还有什么，或者我应该翻译一下these MATLAB functions之类的吗？ ?

我知道如何编写自己的临时实现；我只是在寻找功能齐全的东西，它可以处理不同的窗口功能(但具有合理的默认值)，与 COLA 窗口完全可逆 (istft(stft(x))==x) ，经过多人测试，没有一个错误，很好地处理结束和零填充，快速实现真实输入的 RFFT 等。

Answer 1

这是我的 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 中信号 block 处理的一个小技巧。这就像 Matlab 中的 blkproc。代替 for 循环，我将命令(例如，fft)应用于列表理解中信号的每一帧，然后是 scipy.array 将其转换为二维数组。我用它来制作频谱图、色谱图、MFCC-gram 等等。
2. 对于这个例子，我在 istft 中使用了一种简单的重叠和相加方法。为了重构原始信号，顺序窗口函数的总和必须是常数，最好等于单位 (1.0)。在这种情况下，我选择了 Hann(或 hanning)窗口和 50% 的重叠，效果很好。见 this discussion了解更多信息。
3. 计算 ISTFT 的方法可能更有原则。这个例子主要是为了教育。

测试:

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