1分で試すフーリエ変換（FFT）
投稿者: alvin 日付: 日, 2008-11-16 11:10
はじめに、2つの正弦波(100Hz と 150Hz)を合成します（上グラフ）。
次に、この合成された波形をフーリエ変換して、周波数成分を可視化します（下グラフ）。
このように、100Hz と 150Hz にピークがあることが分かりました。
それでは早速、実習に取り掛かりましょう。
使用する環境は次の3つから選んでください。
R による実習
R による実習では、上に掲載した図のような結果が得られます。
$ R
R version 2.8.0 (2008-10-20)
> n = 0:255 # 256個のデータ
> t = n * 0.001 # 横軸＝t(s,秒)
> f = n * 1 / (256 * 0.001) # 横軸＝f(Hz,周波数)
>
> wave = sin(2*pi*100*t) + sin(2*pi*150*t) # 2つの正弦波を合成
> spec = abs(fft(wave)) # フーリエ変換(FFT)
>
> par(bg="grey95", mar=c(4,4,1,1), mfrow=c(2,1)) # 上下2分割表示、他
> plot(t, wave, type="l", col="darkgreen") # 合成波を描画
> plot(f, spec, type="l", col="navy", xlim=c(0,250)) # スペクトルを描画
>
> q()
$
|
Octave による実習
Octave による実習では、下図のような結果が得られます。
$ octave
GNU Octave, version 3.0.1
octave:1> n = 0:255; # 256個のデータ
octave:2> t = n * 0.001; # 横軸＝t(s,秒)
octave:3> f = n * 1 / (256 * 0.001); # 横軸＝f(Hz,周波数)
octave:4>
octave:4> wave = sin(2*pi*100*t) + sin(2*pi*150*t); # 2つの正弦波を合成
octave:5> spec = abs(fft(wave)); # フーリエ変換(FFT)
octave:6>
octave:6> subplot(2,1,1); # 上側に
octave:7> plot(t, wave, "3"); # 合成波を描画
octave:8> subplot(2,1,2); # 下側に
octave:9> plot(f, spec, "1"); # スペクトルを描画
octave:10> axis([0 250]); # 但し、最大250Hz
octave:11>
octave:11> quit
$
|
Ruby による実習
Ruby による実習では、下図のような結果が得られます。
↓ irb の別名として、irbr を用意しておく
$ alias irbr='irb --simple-prompt -rgsl -rpgplot -rnumru/fftw3'
$ irbr
>> include NMath
>> include Pgplot
>> include NumRu
>>
?> n = NArray[0..255].to_f # 256個のデータ
>> t = n * 0.001 # 横軸＝t(s,秒)
>> f = n * 1 / (256 * 0.001) # 横軸＝f(Hz,周波数)
>>
?> wave = sin(2*PI*100*t) + sin(2*PI*150*t) # 2つの正弦波を合成
>> spec = FFTW3.fft(wave, -1).abs # フーリエ変換(FFT)
>>
?> pgbeg "/xwin", 1, 2 # 上下2分割表示
>> pgpap 4.3, 1.0 # 表示面の大きさを設定
>> pgsci 1; pgenv 0, 0.25, -2, 2 # ↓ X=0〜0.25, Y=-2〜2
>> pgsci 3; pgline t, wave # 合成波を描画
>> pgsci 1; pgenv 0, 250, 0, 100 # ↓ X=0〜250, Y=0〜100
>> pgsci 4; pgline f, spec # スペクトルを描画
>> pgend
>>
?> exit
$
|
覚え書き
-
フーリエ変換(FFT)の結果は複素数になります。
上の実習では、これを abs 関数を使用して絶対値に変換しています。
しかし、
FFT(Octave)
によると、厳密にはパワースペクトラム密度(PSD)を計算する必要があります。
実際のところ、
フーリエ変換(Octave)
にあるように計算してみると、鋭く立ち上がるスペクトルが表示されました。
更新履歴
| 日付 |
内容 |
| 2009-01-10 |
追加 |
Ruby による実習 |
| 2008-11-30 |
追加 |
Octave による実習 |
| 2008-11-16 |
初版 |
R による実習 |
|
