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

13.1 概要

注 – フーリエ変換/逆変換, DFT/IDFT などの定義において, 正規化係数 (1/N など) が文献によって異なることがありますので適宜読み替えてください.

13.1.1 フーリエ変換

\(f(x)\) が複素関数のとき, 次をフーリエ変換といいます.
\[
F(u) = (1/2\pi)\int_{-\infty}^{\infty}f(x)e^{-iux}dx
\] この逆変換をフーリエ逆変換といいます.
\[
f(x) = \int_{-\infty}^{\infty}F(u)e^{iux}du
\]

13.1.2 離散フーリエ変換

有限個の入力関数値(離散値)を使用して数値的にフーリエ変換 (の離散近似) を求めることを考えます.

等間隔にとった \(N\) 個の分点 \(x_j = 2{\pi}j/N (j = 0 \sim N – 1)\) における関数値 \(f(x_j)\) を使用することにし, \(f_j\) と表します.

\(f(x)\) が周期 \(2\pi\) の周期関数であると仮定し, 積分区間を \([0, 2\pi]\) として \(F(u)\) の積分を複合台形公式により求めると次のようになります.
\[
F(u) = (1/N)\sum_{j=0}^{N-1}f_je^{-iux_j}
\] \(u_k = 2{\pi}k, F_k = F(u_k)\) と表すと次のようになります.
\[
F_k = (1/N)\sum_{j=0}^{N-1}f_je^{-i2{\pi}jk/N} \space (k = 0 \sim N – 1)
\] これを離散フーリエ変換 (DFT) といいます. これはフーリエ変換の近似になっています.

同様に離散フーリエ逆変換 (IDFT) は次のように表されます.
\[
f_j = \sum_{k=0}^{N-1}F_ke^{i2{\pi}kj/N} \space (j = 0 \sim N – 1)
\]

13.1.3 実数の離散フーリエ変換

入力関数値が実数の場合には \(F_{N-k} = \overline{F_k}\) となることから \(k = 0 \sim N/2\) についてだけ求めればよいことになります (\(\overline{F_k}\) は \(F_k\) の共役複素数を表します). また, \(F_0\) と \(F_{N/2}\) は実数になります.

通常は, \(a_k = 2Re(F_k), b_k = -2Im(F_k)\) と置き換えて, 次のように実数演算で表します (\(Re\) は実数部, \(Im\) は虚数部を表します).
\[
\begin{align}
& a_0 = (1/N)\sum_{j=0}^{N-1}f_j \\
& a_k = (2/N)\sum_{j=0}^{N-1}f_jcos(2{\pi}jk/N) \space (k = 1 \sim N2) \\
& a_{N/2} = (1/N)\sum_{j=0}^{N-1}(-1)^jf_j \\
& b_k = (2/N)\sum_{j=0}^{N-1}f_jsin(2{\pi}jk/N) \space (k = 1 \sim N2) \\
\end{align}
\] ただし, \(N2 = N/2 – 1\) (\(N\) が偶数の場合), \(N2 = (N – 1)/2\) (\(N\) が奇数の場合).

IDFT は次のようになります.
\[
\begin{align}
& N が偶数の場合:
f_j = a_0 + \sum_{k=1}^{N/2}(a_kcos(2{\pi}kj/N) + b_ksin(2{\pi}kj/N)) + (-1)^ka_N/2 \space (j = 0 \sim N – 1) \\
& N が奇数の場合:
f_j = a_0 + \sum_{k=1}^{(N-1)/2}(a_kcos(2{\pi}kj/N) + b_ksin(2{\pi}kj/N)) \space (j = 0 \sim N – 1) \\
\end{align}
\]

13.2 高速フーリエ変換 (FFT)

DFT および IDFT を少ない演算回数で求める数値計算法を高速フーリエ変換 (FFT) といいます.

DFT を定義どおり計算するならば, \(\sim\) の中で \(N\) 回の計算を, \(N\) 個の \(k\) について行うと全体で \(N^2\) 回のオーダーの計算が必要になります.

しかし, \(N\) が小さな素数のべき乗の積で表される場合には FFT のアルゴリズムを使えばずっと少ない計算量で済みます. 例えば, \(N = 2^n\) のように 2 のべき乗で表されれば \(Nlog_2 N\) 回のオーダーで計算することができます (例: \(N = 1024\) であれば, 通常約100万回必要なのが約1万回でよくなり, 100 倍速くなります).

13.3 XLPack を使った高速フーリエ変換 (FFT) の使い方

VBA サブルーチン Rfft1f, Rfft1b, Rfft1i あるいはワークシート関数 WRfft1f, WRfft1b を使い 1 次元の実フーリエ変換および実フーリエ逆変換を求めることができます.

例題

2.5kHz の搬送波を (100Hz の正弦波 + 200Hz の正弦波) で AM 変調すると次のような波形になります.

このデータを使って DFT および IDFT を求めてみます.

13.3.1 ワークシート関数を使用した求め方

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

DFT を求めるためには, WRfft1f ワークシート関数を使います.

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

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

13.3.2 パワースペクトル

FFT の応用のひとつにパワースペクトル分析があります. これは, スペクトラムアナライザを使って入力信号の周波数分布を観測するのに相当します.

\(N\) 個のデータ \(f_j\) の DFT 出力から, 周波数 \(\omega_k\) におけるパワー \(P_k\) を次のようにして求めることができます. ただし, \(\omega_k = (k/N)\omega_s\) (\(\omega_s\) はサンプリング周波数) です.
\[
\begin{align}
& P_0 = P(f_0) = N |F_0|^2 = N a_k^2 \\
& P_k = P(f_k) = N (|F_k|^2 + |F_{N-k}|^2) = (N/2)(a_k^2 + b_k^2) \space (k = 1 \sim N/2 – 1) \\
& P_N/2 = P(f_N/2) = N |F_N/2|^2 = N a_{N/2}^2 \space (N が偶数の場合) \\
\end{align}
\] 全エネルギー \(E\) はこれらの総和により求めることができます.
\[
E = \sum_{k=0}^{N/2} P_k
\] さて, 上の例題で使用したデータは 20kHz でサンプリングした AM 変調波のデータでした. この 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

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

DFT 結果を入力に指定します.

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

2.5kHz の搬送波の上下 ±100Hz と ±200Hz のところに信号がある AM 変調波になっていることが確かめられました.