Python中的可逆STFT和ISTFT

Question

有没有什么通用的短时傅里叶变换（short-time Fourier transform）和对应的反变换，已经内置在SciPy、NumPy或者其他库里呢？

在matplotlib里，有一个叫做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之外使用而设计的。

我不太确定这个能不能用来做短时傅里叶变换和反变换。不过还有其他选择吗？还是说我应该把像这些MATLAB函数转成Python呢？

我知道怎么写一个临时的实现；我只是想找一个功能齐全的，能够处理不同的窗口函数（但有一个合理的默认值），完全可逆的COLA窗口（istft(stft(x))==x），经过多个人测试过，没有离散误差，能够很好地处理边界和零填充，快速的RFFT实现用于实数输入等等。

Answer 1

这是我使用的短时傅里叶变换（STFT）代码。STFT 加上反短时傅里叶变换（ISTFT）可以实现完美重建（即使是对于第一帧也一样）。我稍微修改了这里 Steve Tjoa 提供的代码：重建后的信号的幅度和输入信号的幅度是一样的。

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

这是我简化过的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图等等。
2. 在这个例子中，我在istft中使用了一种简单的重叠加法方法。为了重建原始信号，连续窗口函数的总和必须是恒定的，最好是等于1.0。在这里，我选择了汉宁窗（Hann window）和50%的重叠，这样效果很好。想了解更多信息，可以查看这个讨论。
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)')

Answer 3

我来得有点晚，但我发现从0.19.0版本开始，scipy里有一个内置的istft函数。