Excel

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

スポンサーリンク

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

Excelのフーリエ解析では解析用のデータを更新しても自動で結果が更新されないので、マクロを使用して自動で更新されるようにしたExcelにつていは、こちら↓を参照ください。

【Excel】高速フーリエ変換(FFT)のマクロ(VBA)

 

C#で作成したフーリエ変換プログラムはこちら↓

【C#】フーリエ変換(FFT, DFT)プログラム

 

そもそも、フーリエ変換って何??

という場合は、下記ページを参照して頂けると幸いです。

フーリエ変換

 

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

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

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

 

ここで、フーリエ変換には処理速度や逆変換を考慮したときに、何を優先させるか?でいくつかの公式があるのですが、エクセルでは下記の真ん中の式が用いられています。

 

Excel フーリエ変換 FFT

 

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

 

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

エクセルではFFTなので、データの個数は2のn乗個で最大4096個まで処理が可能です。

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

 

Excel フーリエ変換 FFT

 

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

 

Excel フーリエ変換 FFT

 

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

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

 

分析ツールの表示

 

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

 

Excel フーリエ変換 FFT

 

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

 

Excel フーリエ変換 FFT

 

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

 

Excel フーリエ変換 FFT

 

この状態で、フーリエ解析を行うデータを選択(2のn乗個の最大4096個まで)し、それでよければ下矢印をクリックします。

 

Excel フーリエ変換 FFT

 

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

 

Excel フーリエ変換 FFT

 

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

 

Excel フーリエ変換 FFT

 

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

周波数0の成分

周波数1の成分

周波数2の成分

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

そもそも i って何?

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

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

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

そのため、

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

 

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

 

詳細は

 

離散フーリエ変換

 

を参照ください。

 

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

 

Excel フーリエ変換 FFT

 

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

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

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

 

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

 

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

Excel フーリエ変換 FFT

 

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

 

COMPLEX(実数、虚数)

 

となります。

 

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

Excel フーリエ変換 FFT

 

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

 

Excel フーリエ変換 FFT

 

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

 

Excel フーリエ変換 FFT

 

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

 

Excel フーリエ変換 FFT

 

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

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

 

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

 

周波数結果パワースペクトル位相(°)
20.70710678118655-0.70710678118647i-45
40.433012701892219+0.25i0.530
280.433012701892219-0.25i0.5-30
300.70710678118655-0.70710678118647i145

 

 

のようになっています。

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

 

Excel フーリエ変換 FFT

 

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

周波数振幅位相
30
245
1-30

 

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

 

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

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

 

Excel フーリエ変換 FFT

 

つまり、複素共役になっている周波数(この例では2と30、4と28)は2つで1つの周波数を表し、振幅も2つを足し合わせることで、振幅が求まります。

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

 

試しにデータの個数32個、周波数(山の数)を0から16まで変化させ、すべて振幅が1の波形を足し合わせた波形

Excel フーリエ変換 FFT
             ↓足し合わせる
Excel フーリエ変換 FFT

 

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

Excel フーリエ変換 FFT

 

これだけを見ると、

 

周波数0と16だけ周波数成分が強くて、他はその半分

 

という見方をしてしまいそうですが、実際の波形は周波数が1と31、2と30、・・・、14と18が複素共役になっていることから、この2つを足し合わせた波形が実際の波形となり、振幅も倍になります。

 

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

周波数N-1の成分を1へ

周波数N-2の成分を2へ

周波数N-3の成分を3へ

周波数N/2+1の成分をN/2ー1へ持っていきます。

Excel フーリエ変換 FFT

 

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

 

Excel フーリエ変換 FFT

 

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

 

前半の部分で、

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

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

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

 

フーリエ変換へ戻る

コメント

  1. 小久保秀之 より:

    私の説明が不十分でした。
    本文全体で「パワースペクトル」と書かれている部分を全部修正した方がいいです。

      「パワースペクトル」→「振幅」
    とすべきところと、
      「パワースペクトル」→「スペクトル」
    とあります。
    ついでですが、折れ線グラフ、棒グラフの縦軸は「振幅」ですね。

    • akira より:

      小久保様、コメントありがとうございます。

      ただ、私の理解が付いていけない部分がありまして。。
      ここでの記事ではデータ個数が32個で周波数を0~16まで変えた例で示していますが、17~31の成分(ここでパワースペクトルと言っている成分)も出てきます。

      入力データが実数の時は、複素共役となる2つの周波数成分を足した時に振幅となりますが、各周波数の成分0~31全てを「振幅」と表現するのには、少し抵抗があるのですが、合っていますでしょうか?

      また、Excelのフーリエ解析で出てくる値は、本記事で最初の方に示した2番目の式が用いられているのですが、フーリエ変換(順変換)の時にはデータの個数で割られないので、やっぱり「振幅」と書くのには抵抗があるもので。。
      (だからと言って、何と言った方が正解かも分かっていないのですが。)

  2. 小久保秀之 より:

    とても参考になりました。ありがとうございます。
    ただ、気になる表現があったのでコメントさせていただきます。
    最後の示された図ですが、フーリエ変換で直接出てくるのは振幅であり、それをパワースペクトルと呼ぶのは適切でないと思います。
    単にスペクトルと言えばいいのでは?

    なお、波のエネルギー(パワー)は振幅の2乗に比例するので、振幅(imabs)の2乗を算出し、さらに複素共役分のパワーを足し合わせればパワースペクトルになります。
    最後の図の例をパワースペクトルで表すと、周波数0の成分は1、その他の成分は0.5になります。

    • akira より:

      小久保様、ご指摘ありがとうございます。
      とりあえず、最後の図に関する「パワースペクトル」の部分を修正しましたが、これで合っていますでしょうか?
      私の記事では、数学的な理解の部分は、自己流なところもあって、学術的にはかなり怪しいかと思います。。

タイトルとURLをコピーしました