Python 中的可逆 STFT 和 ISTFT
-
20-09-2019 - |
题
有没有通用的形式 短时傅立叶变换 与内置于 SciPy 或 NumPy 或其他什么中的相应逆变换?
这是pyplot specgram
matplotlib 中的函数,它调用 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 窗口完全可逆(istft(stft(x))==x
),经过多人测试,无差一错误,处理结束和零填充良好,针对真实输入的快速 RFFT 实现等。
解决方案
我有点晚了,但意识到 scipy 已经内置了 伊斯特夫 从 0.19.0 开始运行
其他提示
这里是我的代码,简化为这个回答:
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
注:
- 的 列表理解 是一个小小的把戏我喜欢用它来模拟框处理的信号,在顽固/这.这就像
blkproc
在Matlab。而不是一个for
循环,我应用的命令(例如,fft
)对各个框架的信号内的一个列表中的理解,然后scipy.array
转换到一个2D阵。我用这个来做谱图,chromagrams,图克,等等。 - 对这个例子中,我使用一个天真的重叠和添加在方法
istft
.为了重建原始信号的总和顺序窗口的功能,必须不断,最好等于统一(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)')
下面是我使用STFT代码。 STFT + ISTFT这里给出的完美的重构(甚至对于第一帧)。我稍微修改由史蒂夫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
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'>)
不过,它们并没有完全颠倒。末端是锥形的。
找到另一个STFT,但没有相应的反功能:
http://code.google.com/p/pytfd/source/browse/trunk/pytfd/stft.py
def stft(x, w, L=None):
...
return X_stft
- w 是一窗口的功能作为一个阵列
- 我 是重叠的,样本
无论是上述的回答运作良好开箱即用的我。所以我修改史蒂夫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...
我还发现此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))
我觉得scipy.signal有你在找什么。它具有合理的默认值,支持多个窗口类型,等等...
HTTP://文档。 scipy.org/doc/scipy-0.17.0/reference/generated/scipy.signal.spectrogram.html
from scipy.signal import spectrogram
freq, time, Spec = spectrogram(signal)
一个固定的版本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()
如果你有机会到你想要做什么,然后用一个C二进制库的http:// code.google.com/p/ctypesgen/ 以产生一个Python接口到该库。