scilab講座その2 フーリエ変換を行う
前回の時系列波形を得てからの続きです。時系列波形の周波数スペクトルを得るためにフーリエ変換を行います。
サンプルデータは結果がわかりやすい様にsin波に代えました。
[wav, sf] = wavread("sin1k.wav");
L=wav(1,:);//処理したい列データを抜く
FFT=abs(fft(L(1:4096),-1));
subplot(2,1,1)
plot(L(1:4096));
subplot(2,1,2)
plot(FFT);
結果はこうなります。
ここの本題はFFT=abs(fft(L(1:4096),-1));で、関数としてfft(データ,-1)で結果が複素数で返ってきます。データは通常は2の乗数のデータを入れます。ここでは4096個としています。
このscilabのfft関数は2の乗数のデータでなくても結果を返しますが、2の乗数の方が高速フーリエ変換を行う為、結果を出す早さが違います。-1というのはフーリエ変換を行うと言う意味で、逆フーリエ変換を行う場合はデータに複素数を放り込んで1にします。
放り込んだデータの数だけの複素数が帰ってきますが、グラフでも2本立っていますように、フーリエ変換が返すデータは真ん中を中心に半分が折り返した値で、前半分しか使えません。前半分しか使えないにせよ、返ってくる値は複素数なので、振幅を見たい場合はabs()で絶対値にする必要があります。
さて、フーリエ変換にはいくつか行わなくてはならない事があり、上のfft()に放り込むデータは始まりと終わりが0に収束しているものなら問題ないのですが、そうではないデータに関しては窓関数をかけて無理やり0に収束させる必要があります。
mado=;//窓関数
mado=window('hn',4096);
scilabはfftをfft()だけで行ってくれる様に、窓関数もwindow()だけで習得してくれます。この場合、'hn'でハニング窓を指定し、4096個の数になるようにしています。
この窓関数をfftに放り込む前にデータにかけてやれば、このデータは端部が0に収束し、フーリエ変換を行った際に「漏れ誤差」を低減させます。
また、fft()で帰ってきた結果をabs()で絶対値にした結果は、「2/データ数」でかけて規格化する必要があります。また窓関数を使った場合は振幅が小さくされてしまうので、補正が必要になります。ハニング窓の場合は2をかけます。
[wav, sf] = wavread("sin1k.wav");
L=wav(1,:);//処理したい列データを抜く
mado=;//窓関数
mado=window('hn',4096);
FFT=[];
FFT=abs(fft(L(1:4096).*mado,-1))*(2/4096)*2;//2/nで正規化。*2で窓ぶんを補正(ハニングの場合)
先ほど使用したplotではなく、周波数軸を入れた配列を作成してplot2dで作図し、見やすいように波形プロットとスペクトルプロットを一部の表示だけにします。
周波数1kHzの振幅1の波形がフーリエ変換されて1kHzの周波数でピークが立ち振幅も1に近く、周波数成分に変換出来たことが判ります。
このようにscilabはフーリエ変換も数行で行います。今回はデータ長を4096で行いましたが、実際のデータはもっと長く、通常はデータをどんどんシフトさせながらフーリエ変換を何回も行い、結果を加算平均してイレギュラーな結果を馴らしていきます。シフトの方法は、例えば、最初は1から4096のデータをfftしたら、次は4097から8193と言う場合はオーバーラップ率0と言い、次は2049から6145だと、オーバーラップ率50%になります。75%が良く使われるようです。
