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

離散フーリエ変換 (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 (更新日: 2019/10/11) | XLPack 入門編
Top