scipy-Python中的可逆STFT和ISTFT

SciPy或NumPy或其他内置的逆傅立叶变换是否有通用形式的短时傅立叶变换？

matplotlib中有pyplot specgram函数，该函数调用2734972705589589167104，该函数调用mlab.specgram()，该函数调用istft(stft(x))==x：

istft(stft(x))==x

但

这是一个辅助功能，可实现   204＃psd，csd和频谱图。 它是   不打算在mlab之外使用

我不确定这是否可以用于STFT和ISTFT。 还有什么，还是我应该翻译像这些MATLAB函数一样的东西？

我知道如何编写自己的临时实现； 我只是在寻找功能全面的东西，它可以处理不同的窗口功能(但默认设置是健全的)，并且可以通过COLA窗口(istft(stft(x))==x)完全反转，并由多人测试，没有错位的错误，可以处理 末端和零填充很好，用于实际输入的快速RFFT实现等。

10个解决方案

62 votes

这是我的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中的istft一样。 我将命令(例如fft)应用于列表推导中的信号的每一帧，而不是hanning循环，然后scipy.array将其转换为2D数组。 我用它来制作频谱图，色谱图，MFCC图等。

对于此示例，我在istft中使用了幼稚的“重叠与加法”。为了重建原始信号，顺序窗口函数之和必须恒定，最好等于1(1.0)。 在这种情况下，我选择了Hann(或hanning)窗口和50％的重叠，效果很好。 有关更多信息，请参见此讨论。

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

Steve Tjoa answered 2020-02-06T02:07:40Z

9 votes

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

Basj answered 2020-02-06T02:08:00Z

3 votes

2734974703168390390和2734974703168390390144看上去与我正在寻找的外观非常相似，尽管当时它们并不存在：

librosa.core.stft(y, n_fft=2048, hop_length=None, win_length=None, window=None, center=True, dtype=)

但是，它们并不能完全反转。 两端是锥形的。

endolith answered 2020-02-06T02:08:29Z

2 votes

我有点迟了，但意识到scipy自0.19.0起具有内置的istft功能

Miss Palmer answered 2020-02-06T02:06:57Z

1 votes

找到另一个STFT，但没有相应的逆函数：

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

def stft(x, w, L=None):

...

return X_stft

w是作为数组的窗口函数

L是重叠，在样本中

endolith answered 2020-02-06T02:09:02Z

1 votes

以上两种答案都不能很好地满足我的需求。 因此，我修改了史蒂夫·乔亚的。

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

capybaralet answered 2020-02-06T02:09:23Z

0 votes

我也在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))

endolith answered 2020-02-06T02:09:47Z

0 votes

我认为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)

meatcomputer answered 2020-02-06T02:10:12Z

0 votes

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

Денис Колесников answered 2020-02-06T02:10:32Z

-3 votes

如果您可以访问所需的C二进制库，请使用[http://code.google.com/p/ctypesgen/]生成该库的Python接口。

Michael Dillon answered 2020-02-06T02:10:52Z