FFTのおさらい

プチコンでFFTを使うときの注意事項と結果についての覚書。

• プチコンのFFTは入出力ともに複素数
• 入出力の要素数は2のn乗である必要がある
• 窓関数が使える。矩形窓、ハミング窓、ハニング窓、ブラックマン窓が使える。指定しない場合は矩形窓ということになる

'FFT TEST ' OPTION STRICT ACLS VAR SAMPLES=512 DIM Ire[SAMPLES],Iim[SAMPLES] ' 入力:実部・虚部 DIM Ore[SAMPLES],Oim[SAMPLES] ' 出力:実部・虚部 DIM Win[SAMPLES] ' 窓関数用 VAR T=64 ' 周期 VAR A=100 ' 振幅 VAR V ' 波形の値 VAR W ' 角周波数 VAR I ' カウンタ ' ' 入力波形を作る（矩形波） FOR I=0 TO SAMPLES-1 W = 2*PI()*(I MOD T)/T V = A* ( SIN(W) + SIN(3*W)/3 + SIN(5*W)/5 + SIN(7*W)/7 + SIN(9*W)/9 + SIN(11*W)/11 ) Ire[I] = V Iim[I] = 0 ' 波形描画（画面サイズに正規化） IF I>0 THEN GLINE (I-1)*400/SAMPLES, 128-Ire[I-1], I*400/SAMPLES, 128-Ire[I] ENDIF NEXT WAIT 60 ' ' 窓関数を設定 ' FFTWFN Win, #WFRECT ' フーリエ変換を実行 FFT Ore,Oim,Ire,Iim ' ' 実部のスペクトル表示(赤) GCLS ? "re" FOR I=0 TO SAMPLES-1 IF ABS(Ore[I])>3 THEN ?I,Ore[I] ' そこそこの振幅がある場合は数値を表示 GLINE I*400/SAMPLES, 128, I*400/SAMPLES, 128-Ore[I], RGB(255,0,0) NEXT WAIT 60 ' ' 虚部のスペクトル表示(緑) GCLS ? "im" FOR I=0 TO SAMPLES-1 IF ABS(Oim[I])>3 THEN ?I,Oim[I] ' そこそこの振幅がある場合は数値を表示 GLINE I*400/SAMPLES, 128, I*400/SAMPLES, 128-Oim[I], RGB(0,255,0) NEXT WAIT 60 ' GCLS FOR I=0 TO SAMPLES-1 Ire[I]=0 Iim[I]=0 NEXT ' ' フーリエ逆変換を実行、Ire,Iimが出力 IFFT Ire,Iim,Ore,Oim ' FOR I=0 TO SAMPLES-1 ' 波形描画（画面サイズに正規化） IF I>0 THEN GLINE (I-1)*400/SAMPLES, 128-Ire[I-1], I*400/SAMPLES, 128-Ire[I] ENDIF NEXT
入力波形は Ire に数値を設定する。虚部 Iim は 0で埋める。
ここでは周期 T=64 の矩形波を作成している。サンプル数が512なので、丁度8周期分が含まれることになる。
これをFFTすると、スペクトルが Ore と Oim に返る。簡単に言うと実部は直流成分とcos成分。虚部はsin成分。
この場合の出力は以下のようになる。
re im 8 -50 24 -16.66666603 40 -10 56 -7.14285564 72 -5.55555534 88 -4.54545307 424 -4.54545307 440 -5.55555534 456 -7.14285564 472 -10 488 -16.66666603 504 -50
この例ではsin成分しかないので、実部の出力は0。
SAMPLES のデータで検出可能な周波数成分は標本化定理により SAMPLES/2 周期まで。つまり512サンプルの場合256サンプル周期の波までが検出できる。
Ore と Oim の 0 から SAMPLES/2-1 までが各周期の成分を示している。8=512/64 なので64サンプル周期の波が -50 の強さで含まれているという意味。ここでのマイナスは位相が現れていることになる。実際は複素平面上の原点からの距離を考える必要があるので、パワースペクトル SQRT(Ore[I]*Ore[I]+Oim[I]*Oim[I]) を求めないといけない。SAMPLES/2からSAMPLES-1 までは周期成分軸上の負の値で、SAMPLES-1 が軸上の-1、SAMPLES/2 が -SAMPLES/2。sin 成分は点対称、cos成分は線対称に現れる。そのため、周波数のパワーのみを考える場合は半分だけ見ればよい。パワースペクトルの総和が元の波形の振幅。
逆変換はスペクトルの数値をそのまま入力にいれてやると、出力は振幅も復元される。