人工知能に関する断創録

このブログでは人工知能のさまざまな分野について調査したことをまとめています(更新停止: 2019年12月31日)

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でパワースペクトルを求めることもできるようです。上の結果のグラフは、

f:id:aidiary:20120716131950p:plain

・・・。横軸もわからないし、左右対称のままだし、見た目もいまいちなので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

結果は、

f:id:aidiary:20120716132129p:plain

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はハミング窓を意味します。他にもハニング窓、ブラックマン窓、バートレット窓、三角窓などが実装されているようです。上のコマンドの実行結果は、

f:id:aidiary:20120716132423p:plain

窓関数をかけたことで裾野が狭くなり、ピークが急峻になることがわかります。

音声データの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してます。結果は、

f:id:aidiary:20120716134153p:plain

440Hzと880Hzにピークがあることがわかります。