Python中的可逆STFT和ISTFT
Posted
技术标签:
【中文标题】Python中的可逆STFT和ISTFT【英文标题】:Invertible STFT and ISTFT in Python 【发布时间】: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 【参考方案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 接口。
【讨论】:
以上是关于Python中的可逆STFT和ISTFT的主要内容,如果未能解决你的问题,请参考以下文章
如何从 librosa 中的 mel 频谱图重建 STFT 矩阵,以便重建原始音频?