研究と教育と追憶と展望

ちまたには高速フーリエ変換（FFT）のプログラムが出回っているが、これをExcel VBA（いわゆるマクロ）で使えるようにしたところ、机上作業の効率が飛躍的にアップした。下記にプログラム等をノートしておく。ちまたに出回っているFFTのプログラムは利用者の知らないうちに窓関数がかけられていたり、周波数を自動的に計算したりするので、あまり教育的ではない。一番シンプルなFFTのプログラムで、生の処理を体験すべきである。

ちなみに、フーリエ変換とは、理工系の大学で習う数学的処理の方法の一つで端的に言えば、横軸が時間となっているグラフをフーリエ変換すると、横軸が周波数のグラフとなるもの。

例えば左図は横軸が時間で、周波数の異なるサインカーブを３つ足し合わせたグラフである。0.001秒毎に（サンプリング周波数1 kHz）、1024個のデータを取得したもの。横軸が時間だと、この波の周波数を読み取るのに一苦労である。

これをフーリエ変換すると左のようになる。横軸が周波数となって、何Hzの波が含まれているかがクリアにわかる。時間領域と周波数領域の間を互いに変換するのがフーリエ変換である。

このフーリエ変換のプログラムを、CooleyとTukeyの発見による感動的なアルゴリズムで高速化したのが高速フーリエ変換（Fast Fourier Transform, FFT）である。「高速」のつかないフーリエ変換のプログラムもあるにはあるが普通は用いられない（データ数が2の整数乗個でなくてもよいという点だけがメリット）。

注意事項

• 高速フーリエ変換を行う元のデータも、変換後のデータもすべて複素数。実際には元のデータが複素数であることはまずないので、虚数部分を0とする。変換後のデータは、複素数の絶対値を求めて使用する。
• データの数は2の整数乗個でなければならない。例 1024, 2048, 4096
• 変換後の周波数領域のデータは左右対称となるので左半分を使用する。右半分は削除する。横軸の右端はサンプリング周波数の半分となる。
• サンプリング周波数の半分以上の周波数成分が信号に含まれていたとしても、観察されるはずがないことを直観で理解して欲しい。同じようにサンプリングした時間よりも周期が大きい波、つまり時間領域のグラフに波１個さえ入れない低周波成分も出てくるはずがない。特に後者は要注意で上の周波数領域のグラフで原点に向かって立ち上がりがあるが、別に低周波成分が含まれていることを示すわけではない。サンプリングしている間に、波が複数入ってないと周期的変動として検知できないのである。

使用法

1. プログラム中にデータ数を指定する（5行目のnの部分と7行目の括弧の中）。下のプログラム例ではデータ数は1024個になっている。データの個数は2の整数乗でなければならない。
2. ExcelのA1～A1024に元データの実数部分、B1～B1024に虚数部分を入力する。ただし、虚数部分があるデータはまれだと思われるので、B1～B1024にはゼロを入力しておけばよい。実数部分で、データ数が2の整数乗に足りない場合は余ったところにゼロを入力しておく人もいれば、同じデータを繰り返しておく人もいる（厳密にはよろしくないと思われる）。
3. Excelの［表示］→［ツールバー］→［Visual Basic］で、VBAのメニューバーを表示、Visual Basic Editorのアイコンをクリックする。Visual Basic Editor上で［挿入］→[標準モジュール]、さらに［挿入］→［プロシージャ］でSub、Publicを選択し、名前を適当に入力して、下記プログラムをコピー＆ペースト（最初と最後の１行は重複するので削除）
4. メニューバーからプログラムを実行する（▲が横向きのアイコン）。
5. 4列目（Dの列）にフーリエ変換結果の実数部分、5列目に虚数部分、6列目に（複素数の）絶対値が表示される。
6. 普通は6列目を使用する。グラフを描くと、左右対称になっていることがわかるが、左半分のみを使用する。横軸の左端は0 Hzとなるが、左半分の右端（左右対称のグラフだと中央、512個目のデータ）はサンプリング周波数の半分に対応する。例えば、元のデータが0.001秒毎に取得したデータであれば、サンプリング周波数は1000 Hzなので、グラフの右端は500Hzとなる。この500Hz、つまりサンプリング周波数の半分のことをナイキスト周波数と呼ぶこともある。

Public Sub fft()
Dim g, h, i, j, k, l, m, n, o, p, q As Integer
i = 0: j = 0: k = 0: l = 0: p = 0: h = 0: g = 0: q = 0
'データ数nを指定（2の整数乗である必要あり）
n = 1024: m = Log(n) / Log(2)
'xr,xiはデータ数以上、s,cはデータ数の半分以上を指定しておく
Dim xr(1024), xi(1024), xd, s(512), c(512), a, b As Single 
'データの読み込み　１列目を実数部、２列目を虚数部とする
For i = 1 To n
xr(i - 1) = Cells(i, 1): xi(i - 1) = Cells(i, 2)
Next i
'FFTの計算
a = 0: b = 3.14159265359 * 2 / n
For i = 0 To n / 2
s(i) = Sin(a): c(i) = Cos(a): a = a + b
Next i
l = n: h = 1
For g = 1 To m: l = l / 2: k = 0
For q = 1 To h: p = 0
  For i = k To l + k - 1: j = i + l
   a = xr(i) - xr(j): b = xi(i) - xi(j)
   xr(i) = xr(i) + xr(j): xi(i) = xi(i) + xi(j)
   If p = 0 Then
    xr(j) = a: xi(j) = b
   Else
    xr(j) = a * c(p) + b * s(p): xi(j) = b * c(p) - a * s(p)
   End If
   p = p + h: Next i
k = k + l + l: Next q
h = h + h: Next g
j = n / 2
For i = 1 To n - 1: k = n
If j < i Then
xd = xr(i): xr(i) = xr(j): xr(j) = xd: xd = xi(i): xi(i) = xi(j): xi(j) = xd
End If
k = k / 2
Do While j >= k
j = j - k: k = k / 2
Loop
j = j + k
Next i
'データの出力
For i = 1 To n   '４列目に実数部、５列目に虚数部、６列目に絶対値を表示
Cells(i, 4) = xr(i - 1): Cells(i, 5) = xi(i - 1)
Cells(i, 6) = (xr(i - 1) ^ 2 + xi(i - 1) ^ 2) ^ 0.5
Next i
End Sub

先人の業績をExcel VBAで使えるようにしたものなので、ご自由にご利用いただいて構いません。

参考文献

科学計測のための波形データ処理（南茂夫編著、CQ出版）　amazonへリンク　古書としてのみ入手可能・・・本書は名著中の名著。FFTに関しても面白くてためになる。しかし、PCの方が進歩しすぎたためか、残念ながら絶版。

ニューメリカルレシピ・イン・シー　Ｃ言語による数値計算のレシピ（技術評論社） amazonへリンク・・・これも名著。内容が充実している。FFTでは実際にはデータが実数で、虚数部分のゼロが計算の無駄になるので、計算を工夫して、２つ分のフーリエ変換を同時に行うプログラムが同書に提案されている。

補足　逆フーリエ変換（周波数領域→時間領域）について

逆フーリエ変換は、上と同じプログラムで実行できる。ただし、注意点として、元となる周波数領域のデータ（複素数）は、左右対称のデータを用いる必要がある（フーリエ変換直後のデータと同様、右半分と左半分が対称なものを作成する）。また、縦軸の倍率が変わるので、縦軸の数値を全データ数（1024とか4096とか）で割った上、正負を逆にする必要がある。

時間領域のデータから出発して、フーリエ変換→逆フーリエ変換を行って、元のデータと比較すれば理解できる。フーリエ変換→逆フーリエ変換を行った後、時間領域のデータの虚数成分は完全なゼロにはならず、ゼロに近い数値になるが、これは誤差のためである。

補足　C言語のFFTについて

上述の技術評論社のニューメリカルレシピ・イン・シー（C言語による数値計算のレシピ）のp.390に記載されています。