Python中的可逆STFT和ISTFT答案

【问题标题】：Invertible STFT and ISTFT in PythonPython中的可逆STFT和ISTFT
【发布时间】：2011-01-28 09:26:52
【问题描述】：

是否有任何通用形式的short-time Fourier transform 以及相应的逆变换内置到 SciPy 或 NumPy 或其他任何东西中？

matplotlib中有pyplotspecgram函数，调用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。还有什么，或者我应该翻译these MATLAB functions之类的东西吗？

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

【问题讨论】：

我正在寻找完全相同的东西，类似于Matlab的“频谱图”功能。
@khpeek 见matplotlib.org/api/mlab_api.html#matplotlib.mlab.specgram
SciPy 现在有这个：scipy.github.io/devdocs/generated/scipy.signal.stft.html

标签： python scipy fft signal-processing

【解决方案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

注意事项：

列表理解 是我喜欢用来模拟 numpy/scipy 中的信号块处理的一个小技巧。这就像 Matlab 中的blkproc。我没有使用for 循环，而是将命令（例如fft）应用于列表推导式中信号的每一帧，然后scipy.array 将其转换为二维数组。我用它来制作频谱图、色谱图、MFCC-gram 等等。
对于这个例子，我在istft 中使用了一种简单的重叠和相加方法。为了重建原始信号，顺序窗口函数的总和必须是常数，最好等于单位 (1.0)。在这种情况下，我选择了 Hann（或 hanning）窗口和 50% 的重叠，效果很好。请参阅this discussion 了解更多信息。
可能有更多原则性的方法来计算 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)')

【讨论】：

网上有没有简化版可以链接的？
不是我的头顶。但是上面的代码有什么问题吗？如有必要，您可以对其进行修改。
不，但你说“简化了这个答案”，所以我认为这是你写的其他东西的删节版
很抱歉给您带来了困惑。是的，从我的应用程序特定需求简化。示例功能：如果输入是立体声信号，则先使其成为单声道；在给定的频率和时间范围内绘制频谱图；绘制对数谱图；将framesamp 舍入到最接近的二的幂；将 stft 嵌入到 Spectrogram 类中；等等。您的需求可能会有所不同。但是上面的核心代码仍然可以完成工作。
感谢此代码。只是一句话：如果 x 不是 hop 长度的倍数，stft 会发生什么？最后一帧不应该补零吗？

【解决方案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

【讨论】：

你能解释一下结果是什么吗？简而言之。我使用了你的代码，它可以工作，但不知道如何解释它......

【解决方案3】：

librosa.core.stft 和 istft 看起来与我正在寻找的非常相似，尽管它们当时并不存在：

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

不过，它们并没有完全反转；末端是锥形的。

【讨论】：

【解决方案4】：

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

【讨论】：

是的，它是最近添加的。 github.com/scipy/scipy/pull/6058 不过我想这应该是公认的答案。

【解决方案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 是重叠，以样本为单位

【讨论】：

我已经测试过这段代码。它冻结了我的计算机以获取大型数据集。 Steve Tjoa 的实现效果要好得多。

【解决方案6】：

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

【讨论】：

【解决方案7】：

我也在 GitHub 上找到了这个，但它似乎是在管道而不是普通数组上运行的：

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

【讨论】：

【解决方案8】：

我认为 scipy.signal 有你正在寻找的东西。它具有合理的默认值，支持多种窗口类型等...

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)

【讨论】：

虽然github.com/scipy/scipy/issues/5757#issuecomment-191516652没有反函数

【解决方案9】：

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

【讨论】：

只是为了确定，它究竟修复了什么？（有什么错误吗？）

【解决方案10】：

如果您有权访问执行所需操作的 C 二进制库，则使用 http://code.google.com/p/ctypesgen/ 生成该库的 Python 接口。

【讨论】：