49

SciPyやNumPyなどに組み込まれた対応する逆変換を備えた汎用形式の短時間フーリエ変換はありますか?

specgrammatplotlib には、 を呼び出す pyplot 関数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 実装など。

4

10 に答える 10

65

これが私の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

ノート:

  1. リスト内包表記は、numpy/scipy でシグナルのブロック処理をシミュレートするために使用するのが好きなちょっとしたトリックです。それblkprocはMatlabのようなものです。forループの代わりに、コマンド (例: fft) をリスト内包表記内の信号の各フレームに適用し、scipy.arrayそれを 2D 配列にキャストします。これを使用して、スペクトログラム、クロマグラム、MFCC グラムなどを作成します。
  2. この例では、単純なオーバーラップ アンド アド メソッドを で使用しistftます。元の信号を再構築するために、連続するウィンドウ関数の合計は一定でなければならず、できれば 1.0 に等しい必要があります。この場合、ハン (またはhanning) ウィンドウと 50% のオーバーラップを選択しましたが、これは完全に機能します。詳細については、このディスカッションを参照してください。
  3. 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)')

440 Hz 正弦波の STFT 440 Hz 正弦波の始まりの ISTFT 440 Hz 正弦波の終わりの ISTFT

于 2011-07-31T19:30:16.770 に答える
9

これが私が使用する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
于 2013-12-05T19:38:28.050 に答える
1

別の STFT が見つかりましたが、対応する逆関数はありません:

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

def stft(x, w, L=None):
    ...
    return X_stft
  • wは配列としてのウィンドウ関数です
  • Lはサンプル単位のオーバーラップです
于 2010-03-18T13:54:53.303 に答える
0

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)
于 2016-04-06T18:55:32.120 に答える
0

これも 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))
于 2010-03-17T16:13:47.733 に答える
-3

必要な処理を実行するCバイナリライブラリにアクセスできる場合は、http://code.google.com/p/ctypesgen/を使用してそのライブラリへのPythonインターフェイスを生成します。

于 2011-08-02T01:31:17.250 に答える