Python中的可逆STFT和ISTFT

Question

是否有任何通用形式的短时傅里叶变换,其中相应的逆变换内置于SciPy或NumPy或其他任何东西？

specgram在matplotlib中有pyplot 函数,它调用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函数吗？

我知道如何编写自己的临时实现; 我只是在寻找功能齐全的东西,它可以处理不同的窗口函数(但是有一个合理的默认值),完全可以与COLA windows(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的信号块处理.就像blkproc在Matlab中一样.for我将命令(例如fft)应用于列表scipy.array推导内的每个信号帧,而不是循环,然后将其转换为2D数组.我用它来制作光谱图,色谱图,MFCC-gram等等.
2. 对于这个例子,我使用了一个天真的重叠和添加方法istft.为了重建原始信号,顺序窗口函数的总和必须是常数,优选地等于1(1.0).在这种情况下,我选择了Hann(或hanning)窗口和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 2

这是我使用的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 3

我有点晚了，但意识到 scipy从 0.19.0 开始已经内置了istft函数