短時間フーリエ変換
Pythonで音声信号処理(2011/05/14)
今回は、短時間フーリエ変換(Short-Time Fourier Transform: STFT)を実装してみます。音声信号スペクトルの時間変化を解析する手法です。ある一定の長さの信号サンプルを切り出し、それに窓関数をかけてからフーリエ変換という手順を切り出す範囲を少しずつずらしながら行います。音声を再生しながらリアルタイムにフーリエ変換する必要があるので高速フーリエ変換(2011/6/18)を使ってみます。最終的には、Windows Media Playerなどの音楽プレイヤーでよく見るスペクトルアナライザ(っぽいもの)を作ります。
窓関数
今まで離散フーリエ変換(2011/6/11)や高速フーリエ変換(2011/6/18)を試したときには、切り出した波形サンプルをそのままフーリエ変換していました。しかし、一般的に、切り出した波形に窓関数をかけてからフーリエ変換する方がよいとのこと。フーリエ変換は切り出した範囲が周期関数であることを仮定していますが、適当に切り出すと端っこがきれいにつながらず周期関数にならないという問題があるためです。窓関数をかけることで端っこがなめらかにつながり周期関数にできます。窓関数を使う理由は、窓関数を用いる理由のページがわかりやすかったです。まず、Numpyの窓関数を試してみます。
#coding:utf-8 import numpy as np from pylab import * N = 512 hammingWindow = np.hamming(N) hanningWindow = np.hanning(N) bartlettWindow = np.bartlett(N) blackmanWindow = np.blackman(N) kaiserWindow = np.kaiser(N, 5) subplot(231) plot(hammingWindow) title("Hamming Window") axis((0, N, 0, 1)) subplot(232) plot(hanningWindow) title("Hanning Window") axis((0, N, 0, 1)) subplot(233) plot(bartlettWindow) title("Bartlett Window") axis((0, N, 0, 1)) subplot(234) plot(blackmanWindow) title("Blackman Window") axis((0, N, 0, 1)) subplot(235) plot(kaiserWindow) title("Kaiser Window") axis((0, N, 0, 1)) show()
実行すると、
のようになります。Window functionsを見ると、Numpyには5種類の窓関数があります。ちなみにScipyにはもっと多くの窓関数があります。どういう使い分けをするのか、いまいちよくわかりませんが、実際に使われるのは、ハミング窓とハニング窓(ハン窓)だそうです。
窓関数の効果
次に、窓関数をかけないときとかけたときでスペクトルがどのように変化するかを試してみます。窓関数をかけるには、元の信号に窓関数の信号を単に掛け算するだけです。
#coding:utf-8 import wave import numpy as np import scipy.fftpack from pylab import * if __name__ == "__main__": wf = wave.open("data/combined.wav", "r") fs = wf.getframerate() # サンプリング周波数 x = wf.readframes(wf.getnframes()) x = frombuffer(x, dtype="int16") / 32768.0 # 0-1に正規化 start = 0 # サンプリングする開始位置 N = 512 # FFTのサンプル数 hammingWindow = np.hamming(N) # ハミング窓 hanningWindow = np.hanning(N) # ハニング窓 blackmanWindow = np.blackman(N) # ブラックマン窓 bartlettWindow = np.bartlett(N) # バートレット窓 originalData = x[start:start+N] # 切り出した波形データ(窓関数なし) windowedData = hammingWindow * x[start:start+N] # 切り出した波形データ(窓関数あり) originalDFT = np.fft.fft(originalData) windowedDFT = np.fft.fft(windowedData) freqList = np.fft.fftfreq(N, d=1.0/fs) originalAmp = [np.sqrt(c.real ** 2 + c.imag ** 2) for c in originalDFT] windowedAmp = [np.sqrt(c.real ** 2 + c.imag ** 2) for c in windowedDFT] # 波形を描画(窓関数なし) subplot(221) # 2行2列のグラフの1番目の位置にプロット plot(range(start, start+N), originalData) axis([start, start+N, -1.0, 1.0]) xlabel("time [sample]") ylabel("amplitude") # 波形を描画(窓関数あり) subplot(222) # 2行2列のグラフの2番目の位置にプロット plot(range(start, start+N), windowedData) axis([start, start+N, -1.0, 1.0]) xlabel("time [sample]") ylabel("amplitude") # 振幅スペクトルを描画(窓関数なし) subplot(223) plot(freqList, originalAmp, marker='o', linestyle='-') axis([0, fs/2, 0, 100]) xlabel("frequency [Hz]") ylabel("amplitude spectrum") # 振幅スペクトルを描画(窓関数あり) subplot(224) plot(freqList, windowedAmp, marker='o', linestyle='-') axis([0, fs/2, 0, 100]) xlabel("frequency [Hz]") ylabel("amplitude spectrum") show()
実行結果は、
のようになります。combined.wavは440Hzと880Hzの2つの正弦波を混ぜた信号です。作り方は、正弦波の合成(2011/6/7)を参照。左側はオリジナルの信号とそのスペクトル、右側がハミング窓をかけた信号とそのスペクトルです。今回は、512サンプルを切り出してFFTしてますが、オリジナルの信号は切り出した区間の左端と右端の位置がずれていて周期関数になってません(ちょっとわかりにくいけど)。このように両端が不連続になるとそのスペクトルは440Hzと880Hz以外の成分も含まれてしまい、スペクトルのピークのすそが広がっていることがわかります。一方、窓関数をかけると左端と右端が0に近くなって滑らかにつながり周期関数になります。窓関数をかけるとスペクトルのピークが鋭くなる様子がわかります。
スペクトルアナライザ
最後に、短時間フーリエ変換を用いてスペクトルアナライザを作ります。matplotlibは、単にグラフを描画するだけではなく、アニメーションさせることもできます。アニメーションさせるにはwxPythonというライブラリが必要なので実行する場合は別途インストールしてください。
#coding:utf-8 import sys import wave import numpy as np import scipy.fftpack import matplotlib matplotlib.use('WXAgg') import matplotlib.pyplot as plt wf = wave.open("data/tamsu02.wav" , "r" ) fs = wf.getframerate() # サンプリング周波数 x = wf.readframes(wf.getnframes()) x = np.frombuffer(x, dtype= "int16") / 32768.0 # -1 - +1に正規化 wf.close() fig = plt.figure() sp1 = fig.add_subplot(211) sp2 = fig.add_subplot(212) print len(x) start = 0 # サンプリングする開始位置 N = 512 # FFTのサンプル数 SHIFT = 128 # 窓関数をずらすサンプル数 hammingWindow = np.hamming(N) freqList = np.fft.fftfreq(N, d=1.0/fs) # 周波数軸の値を計算 def update(idleevent): global start windowedData = hammingWindow * x[start:start+N] # 切り出した波形データ(窓関数あり) X = np.fft.fft(windowedData) # FFT amplitudeSpectrum = [np.sqrt(c.real ** 2 + c.imag ** 2) for c in X] # 振幅スペクトル # 波形を更新 sp1.cla() # クリア sp1.plot(range(start, start+N), x[start:start+N]) sp1.axis([start, start+N, -0.3, 0.3]) sp1.set_xlabel("time [sample]") sp1.set_ylabel("amplitude") # 振幅スペクトルを描画 sp2.cla() sp2.plot(freqList, amplitudeSpectrum, marker= 'o', linestyle='-') sp2.axis([0, fs/2, 0, 20]) sp2.set_xlabel("frequency [Hz]") sp2.set_ylabel("amplitude spectrum") fig.canvas.draw_idle() start += SHIFT # 窓関数をかける範囲をずらす if start + N > len(x): sys.exit() import wx wx.EVT_IDLE(wx.GetApp(), update) plt.show()
tamsu02.wavは、波形を見る(2011/5/19)で紹介した音楽ファイルです。元は、MP3なのでffmpegで8000Hz、16bitのWAVEに変換しています。ここでは、512サンプルを切り出してハミング窓をかけてからFFTしてスペクトルを求めています。その後、128サンプルだけ窓をずらしてハミング窓をかけてからFFTというのを何度も繰り返しています。上のプログラムを実行すると下のようにスペクトルの変化が動画で表示されます。