sumisumi

必要があってFFTを使ったので、
メモとして記録しておく。

FFTをネットで検索すると、ファイナルファンタジーなんとか、が
沢山ひっかかり、本当に調べたいものへの障害になる。


■20120726追記　参考ファイルを修正した．

「FFT-memo.pdf」をダウンロード に、メモをまとめたpdfを置く。
「FFT-test.xls」をダウンロードに、サンプルのエクセルファイルを置く。
サンプルシートを2枚にしました．

まだ誤りがあるかもしれない。
■


＝＝＝　フーリエ解析

ひょんなことからフーリエ解析なるものを使うことになった。
理系（工学系？）の大学生ならば一度は通る道である。

一番身近なのは、音の分析だろう。
・高周波成分（高音域）
・低周波成分（低音域）
がそれぞれどれくらい含まれているのか、
を見るときに使うのである。

人間の耳は、
20ヘルツ（1秒間に20回振動）の低音から、
2万ヘルツまで認識できるそうだが、
周波数成分に分解して、それぞれの成分がどれくらい強いのか、という分布をみて音の特徴を見出すことができる。
その周波数ごとの成分の強さをスペクトルと呼び、
横軸に周波数（あるいは周期）をとってスペクトルの強さを縦軸にとったグラフを、スペクトル分布と呼ぶ。

余談だが、
２つの音がオクターブの関係にあるときは、周波数や周期で互いに倍半分の関係になるから、スペクトルをとると、倍半分の周波数でピークを見ることになる。
また、和音も、互いに整数の比（1～5が多い）をとる。
（平均律をとるピアノなどは正確にはそうはならないが、近似的にそうなる。）

フーリエ解析は、時間的に変化する音のような信号を元データから、
各周波数の波の成分ごとの強さを計算する方法である。

ーーー

フーリエ解析の理論と、プログラムの知識があれば、計算するのは難しくはない。また、ケース数が少なく、計算時間さえ気にしなければ、周波数帯が限られていれば、ひたすら積分すれば求められる。

プログラムさえいらない場合もある。

例えば、データは6万くらいまで、周波数が120くらいまでであればエクセル(2003)でも１シートで計算できる。データ列を一番左に、そしてその右に、周波数ごとのSIN関数、COS関数を計算しておき、SUMPRODUCT関数を使って、積の積分
　∫x(t)sin(2πkt/T)dt
をすれば必要な量が求められるのである。

この方法でやるならば何の問題もないのであるが、
長いデータでかつ広い周波数成分にわたる計算が必要ならば、専用のツールを使った方がよい。

ーーー　FFTを使ってみる。

よく使われるのがFFT、高速フーリエ変換、と呼ばれる計算方法とのこと。
私はFFTのアルゴリズムを理解していないし、そのつもりもないのであるが、
出力結果が何を表すのかくらいは知っていないとお話にならない。

4069個（=2の12乗）までの短い信号、短いデータであれば、エクセルで解析することができる。専用ツールでやればよいのだが、「本当になるのか？」と思って使ってみる。

インストールしたままでは使えないので、
アドインに「分析ツール」を追加して使う。
2007の場合は、画面左上のマークから「EXCELのオプション」、2003の場合は、「ツール」「オプション」から、
（場合によってはインストールディスクが必要かもしれない。）

試しに使ってみる。2007の場合は、「データ」から「データ分析」をみてフーリエ解析を選ぶ。

入力(信号・データ）の範囲と、出力の範囲を指定する。
この図では、同じシートの、データの４列右に出力している。

この結果をみてみる、
入力データは、平均値2cm、波の高さ2cm、周期16秒の波で、データ間隔は2秒である。

FFTで出力されたデータは、スペクトル（振幅ではなく、積分値）の値である。
（クリック拡大）　■　表中の説明を再修正しました2012/07/26　■



1行目（k=0）は、周波数成分を持たない平均定常値（正確には平均値xデータ長16秒）である。（周期無限大）
だから、周期に対するスペクトルとして扱ってはいけないはずである。
周波数も0である。これは結構忘れがちだ。
Web上の解説でも、k=1としているものがかなりある。

２行目以降は、その16秒の間の波数kに対してのスペクトル（積分）値であるが、sin,cosではなく、複素数exp(-i2πk/N)で計算している。その関係で、
計算結果③はkのマイナス領域まである。
kの範囲は、－N/2～0～N/2である。

実数のデータを入力した場合、複素数で現れるスペクトルX(k)とX(-k)は共役の関係になる。どちらにもsineとcosineの両方の情報が含まれている。
どちらも複素数の絶対値④は同じになり、その値は、合成した波の振幅x(N/2）として現れている。

FFTの計算では時間刻みは１として計算されている。
本当はΔT=2秒づつなので、
その周波数成分のスペクトル|X(k)|は⑤のように、
FFTの結果にΔTをかけたものになる。

パワースペクトルP(k)を計算する時には、これ⑤を2乗してT0で割ればよい。

また、工業数学の初期に習う離散フーリエ級数の係数⑥を求めるには、
・定常項　a0　は、k=0のスペクトル値X(k=0)をT0で割る。
・そのほかの項は、
　コサインの振幅akに対しては③の実数成分（IMREAL関数）にΔTをかけて（T0/2)で割ればよく、
　サイン　の振幅bkに対しては③の虚数成分（IMAGINARY関数)に-ΔTをかけて（T0/2)で割ればよい。（マイナスを忘れていましたので修正しました．2012/7/26）

ーーー
エクセルのフーリエ解析の実際への応用については、ウェブ上で詳細まで掲載されているものはあまりないようである。（検索方法が悪いだけかもしれないが）

特にデータ刻みや単位の問題、パワースペクトルへの変換、定義との関係など、よくわからなかった。　また、なぜN/2の項は一つしかないのか、なども考えてみればわかるものの、そのあたりの「かゆいところに手の届く」情報は見つからなかった。

ーーー
もちろんツールは万能ではない。
限界もあるが、それより何より、それをわかっているか、という問題である。
私は、わからない場合は、
ツールを使わずに自分で計算し、既知のデータで検算してしまうので、そう言ったことはおこらない代わりに作業時間がかかる。

道具を使う方が早いか、自分で組んだ方が早いか。
いつも迷うところである。