【Excel】フーリエ解析(FFT)

ここではExcelの標準の機能であるフーリエ解析について紹介したいと思います。

 

フーリエ変換は、ざっくり言うと信号の中から各周波数ごとの大きさと位相を求める変換です。

この”信号”と言っている部分が音や画像のような測定データの場合、この”信号（データ）”に対してフーリエ変換を行うのが離散フーリエ変換となります。

しかし、離散フーリエ変換は処理速度が遅いので、離散フーリエ変換を行うんだけど、データ個数を２のｎ乗（2, 4, 8, 16, 32・・・）個に制限することで、高速に処理することができる処理アルゴリズムが高速フーリエ変換[Fast Fourier Transform(FFT)]となります。

 

ここで、フーリエ変換の公式ですが、主に以下の３つの公式があります。下図の公式を見て頂くと分かりますが、異なるのはシグマの頭についているデータ個数Ｎに関する係数のみ異なります。
いくつも公式があるのは、おそらくプログラム的な都合で、

１番目の公式が本来の公式？
２番目の公式は、フーリエ変換では各周波数の大きさが相対的にどの周波数が大きいのか？を解析的に用いる場合が多く、離散フーリエ変換しかしない（逆変換をしない）場合も多いので、1/Ｎの計算を省き処理を軽くしたかった？エクセルでは、この２番目の公式を用いています。
３番目の公式は離散フーリエ変換と逆離散フーリエ変換のプログラムをできるだけ共通にしたかったため？異なるのは e の乗数の符号のみとなります。

 

このエクセルで用いられている式はフーリエ変換後の値に意味合いを持たないので、エクセルでフーリエ変換を行い、その結果を1/Nする（データの個数で割る）ところまで紹介します。

 

まずは、データの準備を行います。

エクセルではFFTなので、データの個数は２のｎ乗個で最大４０９６個まで処理が可能です。

今回はデータの個数は３２個用意しました。

 

 

結果が分かりやすくなるように、全く振動していないデータと、２周期分のデータと４周期分のデータの３つの波形を合計したデータに対してフーリエ変換を行います。

 

 

まず、Excelでフーリエ解析を行うには分析ツールというのを表示する必要があります。

この分析ツールの表示方法は下記のページを参照下さい。

 

分析ツールの表示

 

メニューのデータ→分析ツールを選択します。

 

 

表示されたウィンドウから解析を選択しOKを押します。

 

 

表示されたウィンドウの入力範囲の部分にある上矢印の部分をクリックします。

 

 

この状態で、フーリエ解析を行うデータを選択（２のｎ乗個の最大４０９６個まで）し、それでよければ下矢印をクリックします。

 

 

フーリエ解析の結果の出力先はいくつか選べるのですが、今回は新規ワークシートにするとして、OKボタンをクリックします。

 

 

そしてこれ↓がフーリエ変換された結果です。

 

 

３行目のセルとかは　22.6274169979696-22.6274169979695i と表示されていますが、上から順に

周波数０の成分

周波数1の成分

周波数2の成分

となっているので、３行目の22.6274169979696-22.6274169979695i は周波数が２の成分を表してします。

そもそも i って何？

と以前は思っていたのですが、グラフに表示すると少しはわかりやすいかと思います。

22.6274169979696-22.6274169979695i は実数成分が22.627・・・・で、虚数成分が ー22.627・・・・の複素数となります。

上図の赤い矢印の大きさがパワースペクトルで、矢印の傾きが位相となります。

そのため、

$$パワースペクトル\quad =\quad \sqrt { { (実数成分) }^{ 2 }+{ (虚数成分) }^{ 2 } } $$
$$位相\quad ={ tan }^{ -1 }\frac { 虚数成分 }{ 実数成分 } $$

 

であることが分かってしまえば話が早いかと思います。

 

詳細は

 

離散フーリエ変換

 

を参照ください。

 

このフーリエ変換の出力のままでは、見づらいので、セルの左側に０から始まる連番を追加しておきます。

 

 

この連番と言った部分が、周波数（全部のデータに山が何個含まれるか？）になります。

このままの複素数を用いてパワースペクトルなどを計算してもよいのですが、値そのものに意味合いが無いので、この複素数をデータの個数で割ります。

複素数を割るにはIMDIV関数を用いるのですが、この関数は複素数を複素数で割る関数なので、データの個数（実数）で割るためには虚数成分が０の複素数として計算する必要があり、

 

＝IMDIV(セル, COMPLEX(32,0))

 

と入力すると、複素数を割ることができます。

 

COMPLEX関数は実数成分と虚数成分から複素数を作る関数で

 

COMPLEX(実数、虚数）

 

となります。

 

複素数をデータの個数で割った状態↓

 

この複素数の絶対値（パワースペクトル）を計算するにはIMABS関数を用います。

 

 

絶対値の偏角（位相）を求めるにはIMARGUMENT関数を用いますが計算結果がラジアンなので度に変換するためにDEGREES関数でラジアンから度に変換しています。

 

 

その結果が以下のようになります。

 

 

複素数の絶対値の部分がフーリエ変換のパワースペクトルで複素数の偏角が波形の位相となります。

今回はエクセルのFFTの結果をデータの個数で割りましたが、各周波数がどの程度の割合で含まれているのか？を見るだけであれば、データの個数で割る必要はありません。

 

ここまでで、エクセルでフーリエ変換をする説明は終了なのですが、このフーリエ変換の結果を見ると

 

周波数	結果	パワースペクトル	位相（°）
０	３	３	０
2	0.70710678118655-0.70710678118647i	１	-45
4	0.433012701892219+0.25i	0.5	30
28	0.433012701892219-0.25i	0.5	-30
30	0.70710678118655-0.70710678118647i	1	45

 

 

のようになっています。

このパワースペクトルの部分をグラフにするとよく見るグラフとなります。

 

 

ただ、もともとのフーリエ変換する前のデータは

周波数	振幅	位相
０	3	0
２	2	45
４	1	-30

 

で作っているのに周波数２と４のパワースペクトルは振幅の半分になってるし、余計な２８と３０の周波数の成分が追加されています。

 

フーリエ変換の結果を見ると、周波数の２と３０、４と２８は虚数の成分だけが正負が逆になっていますが、この関係を複素共役と言います。

この２つの成分を足し合わせると虚数成分が打ち消しあって実数成分が２倍になり、元の成分に戻るのですが、イメージで表現するとこんな感じ↓になります。

 

 

つまり、複素共役になっている周波数（この例では２と３０、４と２８）は２つで１つの周波数を表し、振幅も２つを足し合わせることで、振幅が求まります。

そのため、パワースペクトルの値は共役になっているのか？いないのか？で実際の波形の振幅と倍、異なるので、注意が必要です。

 

試しにデータの個数３２個、周波数（山の数）を０から１６まで変化させ、すべて振幅が１の波形を足し合わせた波形

↓足し合わせる

 

この波形にフーリエ変換を行いパワースペクトルを求めると、このようになります↓

 

これだけを見ると、

 

周波数０と１６だけ周波数成分が強くて、他はその半分

 

という見方をしてしまいそうですが、実際の波形は周波数が１と３１、２と３０、・・・、１４と１８が複素共役になっていることから、この２つを足し合わせた波形が実際の波形となり、振幅も倍になります。

 

以上のことから、パワースペクトルのグラフを見るときは、データがN個の場合、N/2周期までを見て、

周波数N－１の成分を１へ

周波数N－２の成分を２へ

周波数N－３の成分を３へ

：

周波数N/２＋１の成分をN/２ー１へ持っていきます。

 

そのようにした、この状態↓が、各周波数の振幅と一致した状態となります。

 

 

ちなみに、周波数０は波形の山の個数が０個、つまりデータ全体の平均値となります。

 

前半の部分で、

フーリエ変換の結果をデータの個数で割らないと値そのものに意味合いが無い

と書きましたが、上記グラフのように、データの個数で割るとスペクトルの値が各周波数の値と一致するので、個人的にはデータ個数で割った方が好みです。

ただ、フーリエ変換を行う場合は、パワースペクトルの値そのものは関係なくて、各周波数成分のどの成分が大きいのか？というような相対的な見方しかしない場合も多いので、わざわざ個数でわる必要が無い場合もあるので、実際に個数で割るか割らないかは、好みというか、使い方次第です。

 

←フーリエ変換へ戻る

関連記事