Python の可逆 STFT と ISTFT
-
20-09-2019 - |
質問
汎用的な形式はありますか 短時間フーリエ変換 対応する逆変換は SciPy や NumPy などに組み込まれていますか?
パイプロットがあります 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
)、複数の人によってテストされ、off-by-one エラーがなく、端とゼロ パディングを適切に処理し、実際の入力に対する高速 RFFT 実装などが挙げられます。
解決
少し遅れましたが、scipy が組み込まれていることに気づきました。 イストフト 0.19.0 以降の機能
他のヒント
この回答のために簡略化された私の 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 で信号のブロック処理をシミュレートするために使用したいちょっとしたトリックです。まるで
blkproc
マットラボで。の代わりにfor
ループにコマンドを適用します (例:fft
) をリスト内包表記内の信号の各フレームに適用し、その後scipy.array
それを 2D 配列にキャストします。私はこれを使ってスペクトログラム、クロマグラム、MFCC グラムなどを作成します。 - この例では、単純なオーバーラップアンド追加メソッドを使用します。
istft
. 。元の信号を再構成するには、逐次ウィンドウ関数の合計が一定である必要があり、できれば 1.0 に等しくなければなりません。この場合、私はハン (または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 配列としてのウィンドウ関数です
- L サンプル内のオーバーラップです
上記の回答のどちらもが私のためによくOOTBを働きました。だから私はスティーブ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インタフェースを生成します。