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

以前、Excelのマクロを使って、データ個数に応じて高速フーリエ変換(FFT)と離散フーリエ変換(DFT)の処理を自動で切り替えるマクロを作成したのですが、Excelではデータ数が多い時など、使いにくい場合もあるので、今度は、C#でフーリエ変換部分をライブラリ(*.dll)にし、そのライブラリを使ったプログラムを作成しました。

 

実行ブログラムはこちら FourierTransform.zip

ソースコードは、すべてGitHubで公開しています。

https://github.com/ImagingSolution/FourierCSharp

 

使用方法は以下の通りです。

 

プログラムの実行

zipファイルを解凍すると、

FourierTransform.exe

FourierCSharp.dll

sampledata

となっていますので、FourierTransform.exeをダブルクリックして実行してください。

 

フーリエ変換できるデータフォーマット

データは1行あたりに1つのデータのCSVファイルとなります。

例)

1
1.30939551
1.59096344
1.820538079
1.980767038
2.06331351
2.069803727
2.01139008
1.906994542
1.780482109
1.65716389
1.560123288
1.506883103
1.506883103
1.560123288

入力データに実数部と虚数部が含まれる場合は、1行あたりに

実数部、虚数部

とします。

例)

50, 0
0, -25
0, 0
0, -12.5
0, 0
0, 0

サンプルのデータとして、ZIPファイル内の sampledataフォルダ内にいくつかCSVファイルを入れてありますので、そちらを参考にしてください。

 

データの読込

メニューの File → LoadData をクリックし、CSVファイルを指定することで、データが読み込まれ、フーリエ変換の結果(フーリエ変換後の各周波数の大きさ)が表示されます。

 

フーリエ変換はデータの個数が2のn乗個の場合:FFT、その他の場合:DFTを行います。

 

フーリエ変換、逆フーリエ変換の切替

メニューの Fourier direction → Forward もしくは Backward をクリックすることで、離散フーリエ変換、逆離散フーリエ変換が切り替わります。

 

離散フーリエ変換は

$$F(t)=\sum _{ x=0 }^{ N-1 }{ f(x){ e }^{ -i\frac { 2\pi tx }{ N } } } $$

逆離散フーリエ変換は

$$f(t)=\frac { 1 }{ N } \sum _{ t=0 }^{ N-1 }{ F(t){ e }^{ i\frac { 2\pi xt }{ N } } } $$

で計算しています。

 

窓関数

メニューの Window→ Hamming, Hanning, Blackman のいづれかをクリックすることで、それぞれの窓関数を通します。

↓ 窓関数

 

ただし、窓関数を通すと、元に戻せないため、元に戻す場合は、再度、データ読込を行って下さい。

 

免責事項

ここに公開されているプログラムは自由に使って頂いて構いませんが、バグ等による責任は待てませんので、自己責任において参照下さい。

 

フリーウェアへ戻る

フーリエ変換アプリ

こちらのページでは、フーリエ変換のイメージをビジュアル的に見せるためのプログラムになります。

単に、フーリエ変換をしたい場合は、下記のプログラムをご使用下さい。

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

フーリエ変換を教える時に、自分の中では糸(データ)を巻き取るようなイメージ↓(こんな感じ)

 

 

があるのですが、これを分かるように説明するのがなかなか難しく、その様子が分かるように、C#でプログラムを作成してみました。

 

動いている様子はこちら↓

 

実際のプログラムはこちら↓(.NET Framework4.5.2以上で動作すると思います)

FourierTransformAnimation.zip

 

上記ファイルを解凍後、FourierTransformAnimation.exeファイルをダブルクリックすると、プログラムが実行できますが、Windowsの警告画面が表示されるので、「詳細情報」をクリック後、プログラムを実行してください。

また、いくつかサンプルデータを入れてありますので、sampledataフォルダ内のcsvファイルを開いてお試しください。

フーリエ変換で用いている式は、いくつかあろうかと思いますが、下記の式に基づいて処理を行っています。

基本的な機能としては

●CSVファイルによる任意データ入力

●データの離散フーリエ変換

●離散フーリエ変~逆離散フーリエ変換までを動画表示

●窓関数

●離散フーリエ変換の結果をファイル保存

となります。

 

プログラムの説明

データファイルを開くと、離散フーリエ変換を行い、大きさと位相を表示します。

処理のアニメーションの時に用いる複素平面のグラフは、90°反時計周りに回転した状態が初期状態になっています。

離散フーリエ変換の途中途中に出てくる赤い線は、複素平面において、原点からデータの平均の位置までの線で、この線の大きさをMagnitudeグラフへ記載しています。(描画スケールが異なります。)

 

さらにこの線の傾き(12時方向が0°で反時計周りが正)をPhaseグラフへ記載しています。

 

 

逆離散フーリエ変換の時は、各周波数のデータを積算していきますが、前回までのデータの合計を太い青い線で、現在の周波数のデータを赤い線で表示しています。

 

 

使用方法

ファイルを開く

ファイルメニューより、File → Open(*.csv) をクリックし、CSVファイルを選択します。

CSVファイルのフォーマットは1行に1データを縦に記載します。

【例】

1.2
1.107165359
0.914844142
0.80157706
0.872515202
1.061803399
1.193716632
1.145793725
0.962523737
0.814044703
0.838196601
1.012558104
1.175261336

入力データに虚数成分がある場合は、実数と虚数をカンマ(,)でつなげて以下のようにします。

【例】

1,0
0.999390827,0.034899497
0.99756405,0.069756474
0.994521895,0.104528463
0.990268069,0.139173101
0.984807753,0.173648178
0.978147601,0.207911691
0.970295726,0.241921896
0.961261696,0.275637356
0.951056516,0.309016994
0.939692621,0.342020143
0.927183855,0.374606593
0.913545458,0.406736643
0.898794046,0.438371147
0.882947593,0.469471563

データ数がどこまでいけるか?評価していませんが、アニメーション表示をするなら50~100個ぐらいまでが目安です。

アニメーション表示しなければ、そこそこいける?!

 

アニメーションの開始/停止

Youtubeの動画のように、処理の様子を動かすには、ファイルメニューの Animation → Start/Stop をクリックします。スペースキーを押しても同様の動きをします。

 

アニメーションのステップ実行

処理を1つ1つのデータで行う場合は、アニメーションを停止させ、矢印キーの右()を押します。

 

アニメーション速度の調整

アニメーション速度を速くする場合は、矢印キーの上()、遅くする場合は、矢印キーの下()を押します。

 

窓関数処理

基本的に全データを1周期分を想定していますが、この1周期に窓関数(Hamming、Hanning、Blackman)を通します。

ファイルメニューの Window → Hamming、Hanning、Blackman をクリックします。

 

フーリエ変換のイメージ

フーリエ変換では、データを複素平面へ巻き取るようなイメージでになります。

 試しにCosθのデータ全体を1回転で巻き取るようにすると

 

最終的に複素平面の原点から、データが巻き取られた点(複素平面上の点)の平均の位置までの線の大きさが、その周波数の大きさで、線の傾きが位相となります。

 

このデータを巻き取るときの回転の速度が、データ全体を0回転、1回転、2回転・・・で巻き取るようにすると、それぞれの周波数(0,1,2・・・)の大きさと位相が取得できます。

 

試しにCosθのデータ全体を2回転で巻き取るようにしてみると、複素平面上の点の平均は0となります。

もともとのデータにない周波数成分(Cos2θの成分)は、複素平面上で平均を取ると、データが相殺されて0になるのが面白い!

 

逆変換の場合は、取得した各周波数ごとの大きさと位相からコサイン波形を生成し、足し合わせた結果が変換後の実部となります。

ここで、ちょっと違和感があるであろう コサイン波形 と書きましたが、サイン波形を足し合わせた結果は、虚部となります。

 

おそらくこんな書き方をする人はいないのですが、データの個数をN、周波数tの時の大きさをAt、位相をφtとすると、逆変換後の実部(Re)と虚部(Im)は以下のようになります。

 

$$Re(x)\quad =\quad \sum _{ t=0 }^{ N-1 }{ { A }_{ t }Cos(\frac { 2\pi t }{ N } x+{ \varphi }_{ t }) } \\ Im(x)\quad =\quad \sum _{ t=0 }^{ N-1 }{ { A }_{ t }Sin(\frac { 2\pi t }{ N } x+{ \varphi }_{ t }) } $$

 

以下の画像は周波数1の結果を逆変換したものになります。

これだけを見ると、実部のグラフのcosθの振幅は小さくない?!

と思われたかもしれませんが、フーリエ変換を行う元々のデータに虚数部分が含まれない場合、大きさが同じで位相がマイナスの関係になった周波数が存在します。

こんな感じ↓

 

この2つの波形を足し合わせる事で、逆変換では元の波形に戻ります。

この部分の詳細を知りたい場合は、複素共役で調べてもらうと分かるかも?しれません。

 

このフーリエ変換を発明したジョゼフ・フーリエさんは、自分自身をミイラのように包帯でグルグル巻きになるという、ちょっと変わった趣味の持ち主だったそうですが、このグルグル巻きになった時にフーリエ変換が思いついたのではないのかな?(詳細は不明)

 

フーリエ変換へ戻る

 

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

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

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

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

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

 

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

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

N :   データ個数
t:   全データに含まれる波の個数(何周期分か?)
x :   データのインデックス番号 x = 0, 1, 2,… , N-1

データの準備

エクセルのフーリエ解析はFFTなので、データの個数は2のn乗個(2, 4, 8, 16, 32, 64・・・)に限られ、最大4096個まで処理が可能です。

今回はデータの個数は32個用意しました。
実際にフーリエ解析したいデータは各自、用意してください。
データはどこかに1列に書いておけば大丈夫です。(あとで、データの場所を指定します)

結果が分かりやすくなるように、

全く振動していない(t=0)データ(値=3)と、
2周期分(t=2)のデータ(振幅=2、位相=30°)と、
4周期分(t=4)のデータ(振幅=1、位相=-60°)

の3つの波形を合計したデータに対してフーリエ変換を行います。

Excelでフーリエ解析の手順

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

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

【Excel】分析ツールの表示

 

メニューのデータデータ分析を選択します。
ホーム→データ分析にあるデータ分析とは別なので注意してください。

 

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

 

 

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

 

 

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

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

 

 

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

フーリエ変換の結果の見方

※ここから先は、入力データが実数(虚数成分が無い)の場合を前提として説明します。

フーリエ変換の結果は

96
0
27.7128129211021+16i
0
7.99999999999999-13.856406460551i
0
0

のように複素数(実数と虚数成分が含まれる)が並んでいますが、これは上から順にt=0、1、2、…の時の結果になっています。

tの値は、FFTに用いたデータ全体に何周期分の波が含まれているか?を表しています。

tの値から周波数f(Hz)を求めるには、データ全体のサンプリング時間をΔt(sec)とすると、周波数は

$$f=\frac{t}{\triangle t}  (Hz)$$

となります。

ここで、フーリエ変換の結果に、tの値と周波数fを書き加えます。

この結果を見ると、多少の計算誤差はありますが、 t= 2 と t = 30、 t = 4 と t = 28 の値が、虚数成分(i の係数)の正負が異なるだけで、似ている事がわかります。

この実数は同じで虚数の正負が異なる値を複素共役と言います。

フーリエ変換の入力データが実数の場合、必ずこのようになるのですが、t と N – t が共役となるため、フーリエ変換では、 t = 0~N/2 までの結果のみを使う場合がよくあります。

 

フーリエ変換の結果から各周波数の振幅Aを求めるには、フーリエ変換の結果(複素数[実数と虚数])の大きさ |F(t)| を用いると

t = 0 と t = N/2のとき

$$A=|F(t)| / N$$

t = 1 ~ N/2-1 のとき

$$A=|F(t)| / N \times 2$$

となります。

エクセルで複素数の大きさを求めるには IMABS()関数 を用います。

さらに複素数の偏角(位相)を求めるには IMARGUMENT()関数 を用いますが、偏角の結果はラジアン単位なので、度の単位に変換するには DEGREES()関数 を用います。

 

これらの事をエクセルの結果に複素数の大きさ、各周波数ごとの振幅および位相を追加すると以下のようになります。

これで、データをFFTして、各周波数ごとの振幅(A)と位相(φ)を求める事ができています。

各周波数ごと(各tごと)の波形は

$$A_{t} \cos(2\pi t\frac{x}{N}+\phi_{t})$$

となります。

 

今回のFFT解析に用いたExcelファイル(*.xlsx)は、下記リンクをクリックしてダウンロードできます。

フーリエ解析(FFT)_rev1.zip

 

パワースペクトルは??

FFT解析は、各周波数のパワースペクトルを求めること!

みたいな認識もあるかと思いますが、ここまでの説明にパワースペクトルが出てきていません。

パワースペクトルは、FFTで求めた各周波数の結果(複素数)の大きさを振幅スペクトルと言い、振幅スペクトルを2乗したものをパワースペクトルと言います。

パワーという響きから、どうしても力?みたいなのをイメージしてしまいますが、どちらかというと、C言語でn乗の計算をするときに pow関数 を用いると思いますが、このn乗の方に近いでしょうか?(本当のところはよく知りません。。)

 

各周波数ごとに求められたFFTの結果(複素数)を a + bi とすると

$$振幅スペクトル=\sqrt{a^{2}+b^{2}}$$

$$パワーペクトル=(振幅スペクトル)^{2}=a^{2}+b^{2}$$

となることから、プログラム的にはパワースペクトルの計算の方がルートの計算をしていない分だけ計算コストが低くなります。

エクセルでは、IMABS()関数で振幅スペクトルを求め、振幅スペクトルを2乗(^2)することで、パワースペクトルを求める事になり、計算コストが増えてしまい、ちょっと微妙。。

複素共役の性質

入力データが実数のときのFFTの結果で tN – t が共役になると説明しましたが、

複素数 a + bi と、共役の複素数 a – bi と足し合わせると

$$(a+bi)+(a-bi) = 2a$$

となるので、虚数成分は打ち消しあい、実数成分が2倍となります。

この性質から、振幅を求める際に、t = 1 ~ N/2-1 のとき、最後に x 2 を付けているのは、そのためです。

 

 

フーリエ変換へ戻る

参考

【Excel】フーリエ変換

エクセルでフーリエ変換をするには、Webで検索するとほとんどの場合、分析ツールのフーリエ解析で行う方法が紹介されているかと思います。

この方法は以下のページを参照ください。

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

 

しかしながら、分析ツールで行うフーリエ解析では、解析対象のデータを変更しても、自動で結果を再計算しれくれない上に、Excelのフーリエ解析はFFTであるため、データ個数が2のn乗個(2,4,8,16,32,64,128・・・)という制限があります。

 

そのため、エクセルのセルの関数のみを使い、フーリエ変換(離散フーリエ変換[DFT])を行う方法を紹介します。

ただし、結論からすると処理が重いため、データ個数が多い場合はExcelのフーリエ解析(FFT)を使った方が良さそうです。

ここでは、フーリエ変換のお勉強用に使って頂ければ幸いです。

 

エクセルでフーリエ変換を行ったファイルはこちら↓

discrete_fourier_transform_rev1.zip

 

エクセルのイメージはこちら↓

 

データf(x)と書かれた下の黄色いセルの部分にデータを入力すると、フーリエ変換され、振幅スペクトルと位相が計算されます。

 

エクセルの分析ツールで行うフーリエ解析では以下の式を用いていますが、

 

ここでは、フーリエ変換の結果の値に意味合いを持つ下記の式を用いています。

 

 

離散フーリエ変換についてはこちらのページでも紹介していますが、エクセルでフーリエ変換を行うために必要な情報をまとめておきます。

 

まず、フーリエ変換を行うと、各周波数ごとに複素数の成分として変換されます。

 

複素数zは実数成分をa、虚数成分をb、実数Aとすると

      z = a + bi = Ae

 

のように表示され、図示すると下図のようになります。

 

  a = A x Cos θ

  b = A x Sin θ

 

複素数の絶対値は

 

位相(複素数では偏角といいますが)は

 

となります。

このへんの感覚は、ベクトルの成分が(a, b)、ベクトルの長さがA、傾きがθの場合と同じなので、覚えやすいかと思います。

 

 

 エクセルで用いる複素数の関数

エクセルで複素数の計算に用いる関数は関数名の先頭にIMが付く関数を用いますが、今回用いた関数を紹介します。

 

■複素数 COMPLEX関数

実部をa、虚部をbとする複素数は  COMPLEX(a, b)

 

■複素数の絶対値 IMABS関数

複素数(a + bi)の絶対値 IMABS(“a + bi”)

結果は√(a2 + b2) となります。

 

■複素数の実部の取得 IMREAL関数

複素数(a + bi)の実部は  IMREAL(“a + bi”) = a

 

■複素数の虚部の取得 IMAGINARY関数

複素数(a + bi)の虚部は  IMAGINARY(“a + bi”) = b

 

■複素数の偏角の取得 IMARGUMENT関数

複素数(a + bi)の偏角は  IMARGUMENT(“a + bi”)  で-π~πの範囲でラジアン単位で取得

 

■複素数のべき乗 IMEXP関数

zが複素数のべき乗は IMEXP(z)

(例)ei π/3 を計算するには

IMEXP(COMPLEX(0, PI() / 3))

= 0.866 + 0.5i

となります。

 

■複素数の足し算(合計を含む) IMSUM関数

複素数xと複素数yの足し算(片方が実数の場合も含む)は IMSUM(x, y)

(例)

セルA1の複素数とセルA2の足し算は IMSUM(A1, A2)

実数の合計の関数(SUM関数)と同様にセルA1~A10を合計するには IMSUM(A1:A10)となります。

 

■複素数の掛け算 IMPRODUCT関数

複素数xと複素数yの掛け算(片方が実数の場合も含む)は IMPRODUCT(x, y)

 

■複素数の割り算 IMDIV関数

複素数xと複素数yの掛け算(片方が実数の場合も含む)は IMDIV(x, y)

 

 

フーリエ変換の意味合いを読み解く

今回、フーリエ変換の式に

の式を用いたのには、フーリエ変換後の値に意味合いを持つため、あえてこの式を用いました。

 

今回、用いたデータはこちらのページ(http://www.data.jma.go.jp/gmd/risk/obsdl/)より2015年と2016年の東京の2年分(731日分)のデータを用いました。

 

2年分の気温データのため、731個のデータの周波数は2になることを期待しています。

 

そこで、あらためてエクセルのシートを眺めてみると、振幅ペクトルは周波数tが0の時の次に最大になるのはt=2の時で周波数2の成分が大きい事がわかります。

 

このことは、731 / 2 = 365.5日 が1周期分のデータ個数となります。

 

t=0 の時の振幅ペクトルの値(16.4625171)は全体の値の平均値を示しています。

 

t=2 の時の振幅スペクトルの値(5.12963)の2倍の値(10.25926)が周波数が2の時の振幅を示しています。

※実際にはt=2と複素共役と呼ばれる虚数成分がt=2のときの負になるt=729の時の2つの波形を足し合わせることで、虚数成分が打ち消しあい、実数部分の振幅が2倍となります。

 

つまり、今回の2年間の気温データは 平均16.4625171 ± 10.25926 ℃ の傾向で気温変化した事がわかります。

 

t=2 の時の位相(155.3578°)の負の値(-155.3578°)の位置にピークの位置が来ることを意味しています。

 

ここまでわかると、2年間の気温データの変化の主な成分はCOS波形で表現する事ができ

 

16.4625171 + 10.25926 × Cos(x / 365.5 × 360° +155.3578°)

 

となります。

この式を気温データに重ね合わせると↓

 

 

見事、気温変化の傾向がつかめました!!

 

(追記)

フーリエ変換のマクロ処理版のページも作成しました。

合わせてご参照頂けたら幸いです。

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

 

フーリエ変換へ戻る

離散フーリエ変換(DFT)の直感的理解

フーリエ変換に関しては

と説明してきて、ようやく本題の離散フーリエ変換(DFT:Discrete Fourier Transform)となりました。

いきなり離散フーリエ変換って何?と思われる事もあるかと思いますが、画像や音声、AD信号などのデータに対するフーリエ変換は離散フーリエ変換となります。

フーリエ変換と言えば、FFTの方が有名ですよね。FFTは、Fast Fourier Transformの略で、日本語では高速フーリエ変換となります。このFFTは、データの個数を2のn乗個に制限することで、高速に処理を行う処理アルゴリズムであって、基本的には離散フーリエ変換を行っています。
そのため、いきなりFFTを理解するのではなく、まずは離散フーリエ変換を理解すると、FFTの理解も深まります。

そこで以下については、離散フーリエ変換にいての説明です。

 

フーリエ変換に関しては、数々の説明が本やHPでされていますが、なかなか理解するには至らず、ここでは私なりの理解で説明したいと思います。

説明の途中でベクトルという言葉が出てきますが、厳密なベクトルではなく、複素数のイメージで説明しているように複素数を指している場合があるので、ご了承下さい。

e(ネイピア数)を使った説明

※e(ネイピア数)を使わない説明は後半に記載しました。

 

ここで、いきなり離散フーリエ変換の公式ですが、主に以下の3つの公式を目にします。

 

どの公式もΣ以降は同じで離散フーリエ変換はeの指数部分がマイナス、逆離散フーリエ変換のeの指数部分がプラスとなっています。

異なるのはΣの頭に付く係数のみです。

 

3つも公式を並べると、どれを使えばいいの?みたいに思われるかもしれませんが、どれの公式も正解で使う公式によってパワースペクトルの大きさが異なってきますが、各周波数によって得られるパワースペクトルの相対的に大きさには変わりがありません。

しかし、最初の式だけはパワースペクトル(振幅スペクトル)の値に意味を持つので、最初の式を用いて説明します。
※一般には2番目の式が用いられる場合が多いです。(エクセルも2番目の式となります)

【Memo】
離散フーリエ変換を行った結果[F(t)]は複素数となります。
この複素数の絶対値(大きさ)を振幅スペクトルと言います。
振幅スペクトルを2乗したものがパワースペクトルとなります。

まず、公式のf(x)がフーリエ変換を行うデータの配列で、データ個数がN個となります。

が周波数に近い値で全データ中に含まれる波(1周期分の波)の数を表し、データの個数と同じ0~N-1分だけ波が存在します。

フーリエ変換を行った結果F(t)はある特定の周波数tの時の結果で、やはり0~N-1分だけの結果が存在し、この結果は複素数となります。

 

それでは離散フーリエ変換の詳細を見てみたいと思います。

まず、フーリエ変換の公式の複素数(e)の部分に着目します。

 

 

これを見ると複素数のイメージの部分でも説明しましたが、θの部分が

\(\theta=-\frac{2\pi tx}{N}\)

ですが、x はデータのインデックス番号(0, 1, 2…, N-1)で、Nはデータの個数となっています。

これを少々変形し、

\(\theta=-2\pi t\frac{x}{N}\)

とすると、x/Nの部分は、N個のデータ全体(0, 1, 2,…,N-1)が、0~1になるように調整されているので、tの値により、θの範囲が変わります。

t = 0 のとき、θ = 0
t = 1 のとき、0<= θ < -2\(\pi\)
t = 2 のとき、0<= θ < -4\(\pi\)
t = 3 のとき、0<= θ < -6\(\pi\)
:

となります。

つまり、複素空間におけるベクトルがtの値により時計周り(回転角度が負の方向)にt回転するイメージです。

こんな感じ↓

 

tの値が大きくなると矢印の回転も速くなっていくのですが、途中から回転が逆になり、遅くなるように見えます。

この現象はちょうど、走っている車のタイヤをビデオカメラで撮影すると、タイヤが途中で止まって、さらに車の速度が速くなると、タイヤが逆転するように見えるのと似ています。

 

次にf(x)の部分に着目すると複素数(e)の頭の部分にf(x)が付いているので、各xに対応した複素数の大きさ(ベクトル)がf(x)倍される事になります。

 

 

最後にΣでN個の値を合計して、Nで割っているので、複素数を平均している事が分かります。

 

 

ここで、データの個数が16個だけの簡単な例を示します。

f(0)~f(15)のデータを以下のようにします。

 

このデータを周波数t(0~15)ごとに離散フーリエ変換を行います。

(この結果がF(0)~F(15))

 

 

上図のように、複素平面に変換されたデータ(左側の黒い点)の平均へ向かうベクトル(赤い矢印)こそが、離散フーリエ変換した結果となります。

つまり、離散フーリエ変換の結果は各周波数ごとの複素数となります。

 

このベクトル(赤い矢印)の大きさが各周波数における波形の大きさ(振幅スペクトル)でベクトルの向きが各周波数の位相となります。

ただし、一般に、この振幅スペクトルは各周波数の波形の振幅ではない事にご注意下さい。

 

ここで着目して頂きたいのが、t=7 と t=9、t=6 と t=10、t=5 と t=11・・・がベクトル(赤い矢印)の大きさ同じで位相(ベクトルの向き)が±逆となる組み合わせが存在しています。

 

さらにもう少し詳しく見ていきます。

 

フーリエ変換とは?のページで紹介してますが、任意の波形を以下のように定めます。

 

ただし、

N:データ個数

t = 0、1、2、・・・・N-1

n = 0、1、2、・・・・N-1

At:tの時の振幅

φt:tの時の位相

t:データ全体に含まれる1周期分の波の個数。この値で周波数を調整します。

 

この時、波形の位相φが変化すると、フーリエ変換を行うと、どのように変化するのか?

周波数t=2の時、位相を変えながらフーリエ変換を行う例を以下に示します。

 

上図の左側が複素平面で右側が実数平面です。

 

当然ながら、途中に出てくるベクトル(赤い矢印)の角度(位相)が変化します。

 

次に振幅を変えるとどうなるか?やってみると

 

 

途中に出てくるベクトル(赤い矢印)の大きさ(振幅スペクトル)が変化します。

この振幅スペクトルの大きさが、本記事の上部で紹介している離散フーリエ変換の公式の一番最初に出てくる公式を用いると、振幅スペクトルの大きさは信号の振幅の半分の大きさになるのが、好きなとこ。

 

ここまでは、普通な感じでしょうか。

次からが、フーリエ変換を使う大事なところ!?

 

離散フーリエ変換を行う特定の周波数(ここではt = 2の時)に着目し、それ以外の周波数の信号をフーリエ変換するとどうなるのか?見てみると...

 

 

上図のように、信号の周波数とフーリエ変換を行う周波数とが一致した時のみ振幅スペクトルが出現し、それ以外は振幅スペクトルの大きさが0となります。

上図を見ても分かるように、周波数が一致しないと、複素平面上で中心から放射状にデータを打ち消しあう作用が働き、特定の周波数のみの成分を抽出する事が可能となります。

この辺が、フーリエ変換に複素数を使うのが好きなところ。

 

今度は、各データにバイアス成分というか、直流成分というか、ただのオフセットを加えて変化させてみます。

 

 

途中に出てくるベクトル(赤い矢印)の大きさ(振幅スペクトル)が変化しません!

つまり、信号の周波数が分かっていれば、バイアス成分の影響を受ける事なく、信号の位相を抽出する事ができます。

この性質は、フーリエ変換の位相を使った処理では大事なところです。

 

さらに、周波数t=2の時の例を示してきましたが、この信号に周波数t=2以外の信号を加えたら、周波数t=2の時のフーリエ変換がどうなるか?見ていきます。

 

上図右側の緑色の線が周波数t=2の信号、水色の線が周波数t=2以外の信号、

青色の線が足し合わせた信号です。

 

これを見てもやはり、途中に出てくるベクトル(赤い矢印)の大きさ(振幅スペクトル)が変化しません!

 

もはや複素平面の黒い点はバラバラに見えるのですが、他の周波数の影響を受けずに周波数t=2の信号成分のみを抽出してくれています。

この性質が、よくフーリエ変換を行って、各周波数の振幅スペクトルを解析する事で、周波数の特性を使った処理に使われる要因です。

 

ただし、ここで挙げている例は、少し意図的な部分もあって、合成している波形は全て整数倍で示される周波数の波形を足し合わせた例を示しています。

しかし、実際の信号では、必ずしも整数倍の周波数となる訳もなく、整数倍でない周波数が加わると、複素平面の円の放射方向にデータを打消し合わないなので、振幅スペクトルにも影響が出るのでご注意下さい。

 

と、フーリエ変換については、少し前までは、各周波数成分との正規直交との射影をとる事で、各周波数成分を抽出する!ぐらいの認識でいたのですが、これだといまいちビジュアル的にイメージできません。

 

ここに書いてあるイメージだと、なんだかよく分からなかった複素数もビジュアル的にイメージできて、ようやくフーリエ変換を理解できた気になってます。

 

ちなみに、このフーリエ変換のイメージ、自称 ぐるぐるフーリエ変換

 

e(ネイピア数)を使わない説明

eを使うと、数式的にはシンプルに書けますが、実際のプログラムではeを使えない場合も多いので、sinとcosを使った説明をしたいと思います。

 

離散フーリエ変換は、例えば、周波数が1の場合は、全データをちょうど1回転するようにデータをぐるぐる、複素平面へ変換し、変換したデータの平均(平均すると、実数成分と虚数成分の2つ出てきます)が周波数1の時の離散フーリエ変換の結果となります。

周波数が2の時は、全データをちょうど2回転するようにデータを変換し、平均します。

下図は周波数が2の場合を示していますが、この変換を周波数0~N-1(Nはデータ個数)まで行います。

 

ぐるぐる回す、つまりデータを回転させているのですが、回転といえば、画像処理をしているとアフィン変換に出てくる回転行列がありますよね。

 

$$\left( \begin{matrix} { { x }^{ ‘ } } \\ { y }^{ ‘ } \end{matrix} \right) =\begin{pmatrix} cos\theta & -sin\theta \\ sin\theta & cos\theta \end{pmatrix}\left( \begin{matrix} x \\ y \end{matrix} \right) $$

 

離散フーリエ変換における\(x\)と\(y\)は何に相当するか?というと\(x\)が実部、\(y\)が虚部となります。

そのため、\(n\)番目の実部データを\({ Re }_{ n }\)、\(n\)番目の虚部データを\({ Im }_{ n }\)、周波数\(t\)の時の\(n\)番目の角度を\({ \theta t }_{ n }\)とすると、実部と虚部のデータの回転を行列で表すと

 

$$\left( \begin{matrix} { { Re }_{ t,n }^{ ‘ } } \\ { Im }_{ t,n }^{ ‘ } \end{matrix} \right) =\begin{pmatrix} cos{ \theta t }_{ n } & -sin{ \theta t }_{ n } \\ sin{ \theta t }_{ n } & cos{ \theta t }_{ n } \end{pmatrix}\left( \begin{matrix} { Re }_{ n } \\ { Im }_{ n } \end{matrix} \right) $$

 

となります。

周波数\(t\)の時の\(n\)番目の角度\({ \theta t }_{ n }\)は、データの個数Nの全体で1周期分、2周期分…となるように調整し、

 

$${ \theta t }_{ n }=-2\pi t\frac { n }{ N } $$

 

となります。離散フーリエ変換の場合(逆離散フーリエ変換ではなく)、回転方向がマイナス方向(時計回り)になるため、角度の頭にマイナスが付きます。

 

この \({ { Re }_{ t,n }^{ ‘ } }\) と \({ { Im }_{ t,n }^{ ‘ } }\) の\(n\)について、0~N-1の値を平均すると、周波数\(t\)の時の離散フーリエ変換後の値が求まります。

 

周波数\(t\)の時の離散フーリエ変換後の

実部は

$$\frac { 1 }{ N } \sum _{ n=0 }^{ N-1 }{ { Re }_{ t,n }^{ ‘ } } =\frac { 1 }{ N } \sum _{ n=0 }^{ N-1 }{ \left\{ cos\left( -2\pi t\frac { n }{ N } \right) \times { Re }_{ n }-sin\left( -2\pi t\frac { n }{ N } \right) \times { Im }_{ n } \right\} } $$

虚部は

$$\frac { 1 }{ N } \sum _{ n=0 }^{ N-1 }{ { Im }_{ t,n }^{ ‘ } } =\frac { 1 }{ N } \sum _{ n=0 }^{ N-1 }{ \left\{ sin\left( -2\pi t\frac { n }{ N } \right) \times { Re }_{ n }+cos\left( -2\pi t\frac { n }{ N } \right) \times { Im }_{ n } \right\} } $$

となります。

 

あとは周波数\(t\)に関して t = 0 ~ N-1 の時の実部と虚部の値を求めると、各周波数ごとの実部と虚部の成分が求まります。

残すは周波数ごとの実部と虚部の値から、大きさと位相に変換する場合が多いのですが、これは三平方の定理的に考えればいいだけなので簡単です。

 

 

 

 

 

 

 

 

 

 

 

$$大きさ=\sqrt { { 実部 }^{ 2 }+{ 虚部 }^{ 2 } } $$

$$位相={ tan }^{ -1 }\frac { 虚部 }{ 実部 } $$

 

となります。

 

これで以上ですが、実際の離散フーリエ変換では、二次元フーリエ変換でもない限り?変換する元のデータの虚部(Im)は0(ゼロ)の場合が多いので、プログラム的には

実部を

$$\frac { 1 }{ N } \sum _{ n=0 }^{ N-1 }{ { Re }_{ t,n }^{ ‘ } } =\frac { 1 }{ N } \sum _{ n=0 }^{ N-1 }{ \left\{ cos\left( -2\pi t\frac { n }{ N } \right) \times { Re }_{ n } \right\} } $$

虚部を

$$\frac { 1 }{ N } \sum _{ n=0 }^{ N-1 }{ { Im }_{ t,n }^{ ‘ } } =\frac { 1 }{ N } \sum _{ n=0 }^{ N-1 }{ \left\{ sin\left( -2\pi t\frac { n }{ N } \right) \times { Re }_{ n } \right\} }$$

として計算する場合もあります。

まとめ

以上のように、離散フーリエ変換はデータを周波数に合わせてグルグル巻きにして平均を計算します。

1周期の成分は、全データを1回転させて平均する。
2周期の成分は、全データを2回転させて平均する。
3周期の成分は、全データを3回転させて平均する。
・・・
のようにします。

この時、平均した位置と原点までの距離がその周波数の大きさ、X軸と原点ー平均の位置のなす角度が位相になります。フーリエ変換では、何かと大きさばかりに注目されがちですが、位相も有用です。

周波数という言葉を使うと単位はHz(1/秒)?と思ってしまうかもしれませんが、記事中では厳密な使い方をしていません。単位が必要な場合は、全データを何秒でサンプリングしたのか?が重要になってきます。

例えば、全データを2秒でサンプリングした場合は、
1周期と言っていた周波数は、2秒で1回転するので、0.5Hzとなります。
同様に2周期の場合は2秒で2回転するので1Hz
3周期の場合は2秒で3回転するので1.5Hz
となります。

実際のところ、個人的にはFFTや離散フーリエ変換を画像処理に使う機会は、ほとんどありません。
0~N-1周期分の全てのデータが必要ではなく、ある特定の周波数成分はどの程度あるのか?みたいな使い方は良くします。
そのためFFTよりも離散フーリエ変換の方を部分的に使う事が多いので、ぜひ、離散フーリエ変換を理解して応用できるようになって欲しい!という思いで、この記事を書きました。

 

フーリエ変換へ戻る

複素数のイメージ

離散フーリエ変換を勉強すると、突然

 

e = cosθ + i sinθ

 

みたいな式が突然出てきて、これが何だかよく分からないまま、とりあえず公式だけを覚えてみたり...

しかも i  は「実際には存在しない虚数」みたいに教わったので、存在しない物は、なかなかイメージもしにくい。

 

で、いろいろ調べてみて、今ではなんとなく

 

e は複素平面における 長さ1傾きθ のベクトル

 

ぐらいの認識でいます。

 

 

複素平面ではX軸に相当する部分が実部(Real part)、Y軸に相当する部分が虚部(Imaginary part)と呼ばれます。

 

回転方向は反時計まわりが正で時計まわりが負となります。

 

よく e-iθ というのも目にしますが、これはむしろ ei(-θ) と 書いてくれた方が分かりやすいと思いますが、e-iθ は、回転方向がマイナス、つまり時計まわりになります。

 

また、e の頭に係数が付く場合がありますが、その係数分だけベクトルの長さが変化します。

例えば 2e の場合だと

 

のようにベクトルの長さが2となります。

係数の値が負の場合は、ベクトルの向きが逆になり、-2e の場合だと

 

のようになります。

ただ、実際には e の頭に係数はθの値に応じて変化する場合が多く、例えば

 

(sinθ) e

 

のθの値を0~360°で変化させるとどうなるか?というと、このよう↓になります。

 

 

実はフーリエ変換をイメージで覚えるうえで、上図のように思えた事が最大のポイントで、とても嬉しかった~!!!

 

上図をよ~く見てみると、ベクトルの先端は円を描き、この円の中心が(0, 0.5)の座標になっています。

これを少し言い方を変えると、円の中心へのベクトルは長さ0.5、傾きが+90の位置にあります。

この事は、すでに少しだけフーリエ変換になってます。

 

と、ここまで説明しといて、例題としては  (sinθ) e -iθ  の方が良かったかな...ちょっと失敗??

ちなみに(sinθ) e -iθ  だと、円の中心は(0, -0.5)の位置に来ます。

 

複素数、複素関数について、もう少し詳しく知りたい方は、下記ファイルが参考になると思います。

 

複素関数を学ぶ人のために

http://collie.low-temp.sci.yamaguchi-u.ac.jp/~ashida/work/comp.pdf

 

フーリエ変換へ戻る

関連記事

複素数の計算

フーリエ変換の種類

フーリエ変換と言っても、フーリエという言葉が付く変換はいくつかあります。

 

簡単に分類すると、こんな感じ↓でしょうか?

 

フーリエ変換
    フーリエ変換(Fourier Transform)
    離散フーリエ変換(DFT:Discrete Fourier Transform)
            高速フーリエ変換(FFT:Fast Fourier Transform)

 

フーリエとは付きませんが近いのに離散コサイン変換(DCT:Discrete Cosine Transform)というのもあります。

 

フーリエ変換を大きく分類すると、頭に何も付かないフーリエ変換と離散という言葉の付く離散フーリエ変換とがあります。

 

この頭に何も付かないフーリエ変換は変換相手がsinやcosなどの数式で、フーリエ変換の公式にも積分の式()が用いられます。

 

対して離散フーリエ変換は音声や電圧、画像などのように測定データが対象で離散フーリエ変換の公式には積算の式(Σ)が用いられます。

 

高速フーリエ変換(FFT)は離散フーリエ変換の一種で、データの個数が2のn乗(2,4,8,16,32・・・)個の時に特別に高速に処理をできるアルゴリズムで、基本的な処理は離散フーリエ変換となります。

 

私がフーリエ変換を行う場合は、画像や電圧などの測定データを相手にする場合がほとんどなので、単に「フーリエ変換」と言っておきながら「離散フーリエ変換」の事を指して言っている場合もあるかもしれません...

 

フーリエ変換へ戻る

 

フーリエ変換とは?

フーリエ変換は何をしているのか???

 

下図の示すように、任意波形を各周波数の成分に分解し、その大きさ(振幅)および位置(位相)を求める計算がフーリエ変換となります。

 

逆に、振幅と位相から元の波形を求める計算が逆フーリエ変換となります。

 

このフーリエ変換を行うと何が良いのか???

用途は様々なのですが、

 

  • 振幅から、どの周波数成分が強いのか?解析を行う。
  • 位相から、位相のズレで距離(形状)を算出する。
  • 特定周波数の振幅を0(ゼロ)にして、逆フーリエ変換を行い、フィルタ処理(ノイズ除去など)を行う。
  • 振幅および位相から元のデータを算出する事で、元のデータを少ない値から算出する。
    (データ圧縮)

 

などなど。

 

ここで、私がなかなかフーリエ変換を理解できなかったポイントですが、

 

「任意波形はサイン波形を足し合わせると求める事ができる!」

 

みたいな説明が良くあります。さらに位相についても触れられていない場合が多いです。

 

しかし、この説明だとsin(0°)は、いつも0(ゼロ)なので、いくら各周波数成分を足し合わせても0°の値は0になってしまいます。

 

そのため、sin波形で考えるのではなく、cos波形で考えた方が良いと思います。

この時、各周波数の成分は下記の式で表されます。

 

ただし、

N:データ個数

t = 0、1、2、・・・・N-1

n = 0、1、2、・・・・N-1

t:tの時の振幅

φt:tの時の位相

t:データ全体に含まれる1周期分の波の個数。この値で周波数を調整します。

 

また、一度は聞いた事があるであろうパワースペクトルの値は振幅そのものを求めていない場合が多い。

(この事にはついては、別の記事で説明したいと思います。)

 

画像処理をやった事がある人では、フーリエ変換は、各周波数をテンプレートとしたテンプレートマッチングを行い、相関値がフーリエ変換のパワースペクトルに相当していると思えば、感覚的には似ています。

(ただし、数学的には異なります。)

 

と、ざっくりなイメージはこんな感じです。

 

もう少し、まじめにフーリエ変換を理解したい場合は内積正規直交規定を理解しておくと、良いかと思います。

 

フーリエ変換へ戻る

 

【参考書籍】信号処理入門

この本は、知り合いの人から勧められた本なのですが、私も勧めたくなる本でした。

 

画像処理を勉強すると、とりあえずの集大成的なものがフーリエ変換だと思います。

しかし、いきなり公式が出来てきて何だか良く分からないままFFTまでやってみたりと、とても難しいという印象のまま終わってみたりもします。

この本は、そんなもう一度ちゃんとフーリエ変換を理解したいという人にはオススメです。
図が多めで、そんなに難しくは無いと思います。

 

信号処理の概要から始まり、ベクトル→内積→正規直交基底→フーリエ変換→DFT→FFTという流れで説明されています。

この流れがまさにお見事!

フーリエ変換の公式などは本の後半に書いてありますが、いきなり最後の部分を見るのではなく、最初から丁寧に読むと良いと思います。

 

今ではフーリエ変換って正規直交基底との内積でしょ!?

画像処理だと正規化相関も内積、3×3フィルタみたいのも内積、というぐらいの認識でいます。

内積ってサイコ―!!と思えるようになったのも、この本のおかげです。

 

 

目次

1 信号処理とは
1.1 信号処理はどんなとき必要か
1.2 どんな信号があるか
1.3 アナログ信号とディジタル信号
1.4 サンプリング問題
2 信号処理の例
2.1 波形の平滑化
2.2 雑音の圧縮
3 数学の準備体操
3.1 信号処理を学ぶために
3.2 信号の表現
3.3 2次元ベクトルの距離と内積
3.4 正規直交基
3.5 多次元ベクトル空間から関数空間へ
3.6 正規直交関数系
4 相関関数
4.1 関数の類似性を測る
4.2 相互相関関数
4.3 自己相関関数
5 フーリエ級数展開
5.1 フーリエ級数展開とは
5.2 偶関数と奇関数
5.3 周期が2πでない場合
5.4 複素フーリエ級数展開
5.5 パーシバルの定理
5.6 フーリエ級数展開の実例
5.7 フーリエ級数展開の重要な性質
6 DFTとFFT
6.1 ディジタル信号のフーリエ解析
6.2 離散フーリエ変換(DFT)
6.3 DFTの性質
6.4 高速フーリエ変換(FFT)
7 フーリエ変換
7.1 フーリエ級数展開からフーリエ変換へ
7.2 フーリエ変換の性質
7.3 デルタ関数と白色雑音
8 線形システムの解析
8.1 線形システム解析へのアプローチ
8.2 入出力信号の関係
8.3 インパルス応答
8.4 周波数領域でのシステムの表現

 

 

と言っても、フーリエ変換は、あまり使わないな~

 

使える数学へ戻る