SPTKの使い方 (4) フーリエ変換
SPTKの使い方 (3)(2012/7/7)の続き。
今回は、SPTKでフーリエ変換してみます。SPTKには、fftとfftrという2つのコマンドがありました。fftが複素信号用、fftrが実信号用です。音声信号は実信号なので今回はfftrを使います。
サイン波の作成
まずFFTするサイン波を作成してみます。
sin -l 256 -p 10 > sin.dat
このコマンドで長さ256サンプル、周期10サンプルのサイン波を作成してsin.datへ格納しています。中身は、dmpコマンドで見られます。
dmp +f sin.dat 0 0 1 0.587785 2 0.951057 3 0.951057 4 0.587785 5 1.22465e-16 6 -0.587785 7 -0.951057 8 -0.951057 9 -0.587785 10 -2.44929e-16 11 0.587785
sinコマンドの戻り値は、-1.0から1.0のfloatで返るみたい。ちゃんと10サンプルごとでもとに戻ることも確認できます。
FFT
フーリエ変換は簡単にできます。波形データをfftrコマンドに流し込むだけ。
sin -l 256 -p 10 | fftr -l 256 -A | fdrw | xgr
-lはFFTのサンプル数、-Aは振幅スペクトルを出力することを意味します。他にも-Iで位相スペクトル、-Pでパワースペクトルを求めることもできるようです。上の結果のグラフは、
・・・。横軸もわからないし、左右対称のままだし、見た目もいまいちなのでmatplotlibを使って描画してみよう。
#coding:utf-8 # plot_fft.py # SPTKのFFT結果をプロットする import struct import sys from pylab import * if len(sys.argv) != 4: print "usage: python plot_fft.py [fft_file] [N] [fs]" sys.exit() fftfile = sys.argv[1] N = int(sys.argv[2]) fs = int(sys.argv[3]) # FFTファイルをロード fft = [] f = open(fftfile, "rb") while True: # 4バイト(FLOAT)ずつ読み込む b = f.read(4) if b == "": break; # 読み込んだ4バイトをFLOAT型(f)でアンパック val = struct.unpack("f", b)[0] fft.append(val) f.close() # 横軸(Hz) freqList = np.fft.fftfreq(N, d=1.0/fs) # プロット plot(freqList[:N/2], fft[:N/2]) xlim([0, fs/2]) xlabel("frequency [Hz]") ylabel("amplitude spectrum") show()
表示範囲は、0Hzからナイキスト周波数までなので左右対称ではなくなります。横軸をHzで表示するために仮のサンプリング周波数として16000Hzを指定しました。256はFFTのサンプル数です。
sin -l 256 -p 10 | fftr -l 256 -A > fft.result python plot_fft.py fft.result 256 16000
結果は、
1600Hzあたりにピークがあることがわかります。なぜかというと、作成したサイン波の周期は10サンプルなので、サンプリング周波数を16000Hzとすると10サンプルは10 / 16000秒になります。よって、周波数は周期の逆数なので、16000 / 10 = 1600 Hzとなります。
窓関数
次に窓関数をかけてからFFTをしてみよう。窓関数はwindowコマンドをfftrの前に挟むだけでかけられます。windowやfftrに限りませんが、波形データをコマンドのパイプでつなげて処理していくという設計にはすごい感心しました。windowの引数は-lと-Lがありますが、-lは入力のフレーム長、-Lは出力のフレーム長なのでここでは-Lです。
sin -l 256 -p 10 | window -L 256 -w 1 | fftr -l 256 -A > fft.win.result python plot_fft.py fft.win.result 256 16000
-w 1はハミング窓を意味します。他にもハニング窓、ブラックマン窓、バートレット窓、三角窓などが実装されているようです。上のコマンドの実行結果は、
窓関数をかけたことで裾野が狭くなり、ピークが急峻になることがわかります。
音声データのFFT
最後に、sinコマンドの出力ではなく、実際のwavデータをFFTしてみます。使うのは
です。サンプリング周波数は16000Hz、長さは1秒間、440Hzと880Hzのサイン波の混合波です。混合波の作り方は、正弦波の合成(2011/6/7)を参照してください。
wav2raw combined.wav x2x +sf combined.raw | bcut +f -s 0 -e 255 | window -l 256 -w 1 | fftr -l 256 -A > result.fft python plot_fft.py result.fft 256 16000
SPTKのほとんどのコマンドはwavファイルを直接扱えないので、wav2rawコマンドでrawファイルに変換してます。また、rawファイルは -32768から32767のshort型なのに対し、fftrの入力はfloat型なので、x2xコマンドでshort型からfloat型に変換します。dmpするとわかりますが、単にfloat型にキャストするだけで、-1.0から1.0に正規化されるわけではなさそうです。
その後、bcutでFFTをかける範囲の波形を切り出します。ここでは、0サンプル目から255サンプル目までの256サンプルが対象です。あとは、今までと同様に窓関数をかけてからFFTしてます。結果は、
440Hzと880Hzにピークがあることがわかります。