13. 高速フーリエ変換 (FFT)

注 – 例題ワークシートをダウンロードすることができます. 1. 例題ワークシートの入手と使用方法 をご覧ください.

離散フーリエ変換 (Discrete Fourier Transform (DFT)) は, 離散データに対するフーリエ変換を数値的に求めるものです. DFT の逆変換を離散フーリエ逆変換 (Inverse Discrete Fourier Transform (IDFT)) とよびます. DFT/IDFT を少ない演算回数で求める数値計算法が高速フーリエ変換 (Fast Fourier Transform (FFT)) です.

高速フーリエ変換の XLPack による計算例を説明します.

例題 1

400Hz と 1400Hz の正弦波を混合した波を作ると, 次のような波形になります.

13-01

1/10000 秒ごとのデータを用意し, DFT を求めます(変換結果の格納内容についてはリファレンスマニュアルを参照). そして, さらにその IDFT を求めて変換が正しく行われていることを確認します.

DFT を求めるためには, WRfft1f ワークシート関数を使います. 下に50点のデータを使った例を示します.

DFT が求められたら, その内容に対して WRfft1b ワークシート関数により IDFT を適用します.

下のように, 元のデータが復元されたのが確認できます.

例題 2 (パワースペクトル)

FFT の応用のひとつにパワースペクトル分析があります.

N 個のデータ fj の DFT 出力から, 周波数 fk におけるパワーを次のようにして求めることができます. ただし, fk = (k/N) fs, fs はサンプリング周波数です.

  
  P0 = P(f0) = N |F0|2 = N ak2
  Pk = P(fk) = N (|Fk|2 + |FN-k|2) = (N/2)(ak2 + bk2)  (k=1～N/2-1)
  PN/2 = P(fN/2) = N |FN/2|2 = N aN/22  (Nが偶数の場合)

 
全エネルギーはこれらの総和より求められます. 

  
   E = Σ Pk  (Σはk=0～N/2)

 
上の例題の DFT 結果よりパワースペクトルを求めてみます.

ワークシート上で計算すると複雑になるので, パワーを求める次の VBA 関数を用意します.

Function PsdR(N As Long, VR As Variant) As Variant
    Dim R() As Double, P() As Double, I As Long
    ReDim R(N - 1), P(N / 2, 0)
    For I = 1 To N
        R(I - 1) = VR(I)
    Next
    P(0, 0) = N * R(0) ^ 2
    For I = 1 To N - 3 Step 2
        P((I + 1) / 2, 0) = (N / 2) * (R(I) ^ 2 + R(I + 1) ^ 2)
    Next
    P(N / 2, 0) = N * R(N - 1) ^ 2
    PsdR = P
End Function

この関数を使ってパワースペクトルを求めます.

求められたパワースペクトルより作成されたグラフはと次のようになります(ピリオドグラムといいます).

13-06

f(x) には 400Hz と 1400Hz の成分が含まれていることがわかります.

公開日: 2015/10/11 (更新日: 2022/11/03)