sin関数のFFT

scilabではfft(高速フーリエ変換)という関数を用いて、離散フーリエ変換を計算することができる。

x=440Hzのsin波, を、fftで変換してプロットしてみよう。

xx=fft(x,-1); plot2d(n,abs(xx));

Q このプロットでは横軸をサンプル番号とし、縦軸もそのままとした。しかし、本来このプロットの横軸と縦軸はどのようにすべきだろうか?

A 右端がfsになるようにし、振幅は2/l倍すると、sin波の振幅と周波数を読み取ることができる。

実際にやってみよう。

plot2d(fs*n/(2*fs),abs(xx*2/l));

ちょうど目盛が440の部分に高さは0.6のサンプルがあることが読み取れる。これは、440Hzで振幅0.6のsin波が検出されたことを意味する。

それでは以下の3つのsin波を合わせた波形をフーリエ変換してプロットしてみよう。

周波数	振幅
440	0.6
660	0.3
880	0.4

 

x=sin(2*%pi*440*t)*0.6+sin(2*%pi*660*t)*0.3+sin(2*%pi*880*t)*0.4; xx=fft(x,-1); plot2d(fs*n/(2*fs),abs(xx*2/l));

各波形の周波数と振幅が、正確にプロットから読み取れる。

このように、離散フーリエ変換の結果は、各周波数のsin波の振幅を表している。つまり、時間軸のサンプルが、周波数軸のサンプルに変換された、と見ることができる。

そこで、周波数軸に変換されたサンプルを、「周波数サンプル」と呼ぶ。

周波数範囲、周波数間隔

離散フーリエ変換の変換前と変換後の関係を、調べておこう。

tをサンプル間隔、fをサンプリング周波数、Nをサンプル数とする。

変換前と変換後のサンプル数は同じだ。

変換後の周波数範囲は、変換前のサンプリング間隔の逆数になる。

したがって、

周波数サンプルのサンプル間の間隔(周波数間隔) = f/N

となる。

離散フーリエ変換の絶対値の 2/l倍を求めると、元の信号に含まれていたsin波の振幅になる。周波数はサンプル番号の fs/n 倍だ。

対称性

適当な波形の離散フーリエ変換の結果の実部をプロットすると。

x1=abs(x) - (sum(abs(x))/length(x)); xx=fft(x1,-1); plot2d(fs*n/(2*fs),real(xx*2/l),rect=[0,-1,fs,1]);

このように、fft結果の実部は、4000Hzを中心として左右対称になっているように見える。

虚数部をプロットすると、

plot2d(fs*n/(2*fs),imag(xx*2/l),rect=[0,-1,fs,1]);

このように実数部は0で、虚数部は4000Hzを中心に反対称になっている。

離散フーリエ変換の式を見て、注意深く計算を行えば、上記が一般的に言えることも確認することができる。

 

 

独立な数値の数

ところで、2N個の独立な数値の組(独立変数)を誤差を生じずにN個の数値(独立変数)に変換し、また元に戻すことはできないはずである。プロットを見ると左右対称になっている。つまりl個のサンプルから計算したDFTの結果は、l/2個しか独立した値はない。これでは計算上、l個の数値を表すことはできないはずだ。

しかし、DFTの結果は「複素数」なので、1サンプル当たり2個の数値を持っている。だから、互いに独立な数値(独立変数)の数はどちらもl個で同じになる。