ExcelでDSPをしよう(フーリエ変換編)

フーリエ変換を用いると信号を時間領域から周波数領域へ変換し、スペクトルを見る事ができます。
離散信号では特に「離散フーリエ変換(DFT)」と呼ばれます。
実際DFTしようとすると面倒な事になるのですが、なんとExcelにはDFT機能が搭載されています。
これは使うしかないですね。

ExcelのDFTを使うためにはアドインを設定する必要があります。
オプション→アドイン を選択し、その中に「分析ツール」がアクティブになっている事を確かめましょう。
無ければ、Excelアドイン→設定 から「分析ツール」にチェックを打ち、OKボタンを押せばインストールされます。

余談:
DFTにはいくつか手法があるのですが、Excelでは一般的に使われる「FFT(Fast Fourier Transform)」が使われていると予想されます。
FFTとは「バタフライ演算」という演算を用いた手法で、並列演算に特化した構造をなしています。
勘違いしがちですが、FFTとは「DFTの手法の一つ」で、両者は全く同じではありません。

さて、アドインの設定が終わったところで何か試してみたいと思います。
手始めにデルタ関数のスペクトルを見てみましょう。
デルタ関数をフーリエ変換した結果を知っているなら、スペクトルがどうなるか見当がつくでしょう。

※ここから複素数を取り扱いますが、虚数単位はExcelに合わせて「i^2=-1」とします。
電気系の慣例である「j」を使わないので注意してください。
また、単位と離散領域関数(配列のインデックス)の括弧を、どちらも角括弧として書くので間違えないでください。

デルタ関数δ(t)の最大値は∞ですが、Excelにおいて振幅は有限なのでそんな値はとることができません。
ここでは最大値を1に制限したデルタ関数d[t]を使います。
定義は
・d[0]=1
・(それ以外)=0
です。

ここで、ある関数g(t)をフーリエ変換したとき、その操作、変換後のg(t)をそれぞれ
F[g(t)]=g~(f)
と書くことにします。
Zoom パラメータを用意します。
・fs:サンプリング周波数
・n:サンプル番号(0からスタートを推奨)
・t:時間 [s]
・d[t]:制限デルタ関数 [V]
サンプリング周波数は0以外の適当な数値を入れてください。
tのアキュムレータの構成方法は自由です。
d[t]は上に書いた通りです。
それぞれ、128個以上用意しましょう。
Zoom データが用意できたところで、DFTを行います。
DFT機能は
データ→解析→フーリエ解析
から選択できます。
Zoom パラメータを入力するウインドウが出ます。
入力範囲はd[t]の0サンプル〜127サンプルを選択します。
FFTの構造上、2^n個の範囲しか選ぶことができません。
この選択範囲(変換範囲)を「窓、ウインドウ」と言い、窓に収まるデータ個数(窓の大きさ)の単位は普通「ポイント、pts」と表されます。
したがって、今回は128 ptsのFFTです。

出力は「出力先」ボタンを選択し、時間領域信号部分から少し離れた場所を指定します。
設定が完了したら、「OK」をクリックしましょう。
Zoom エラーの付いたセルが128個できるはずです。
エラーは「数字が文字列として書かれている」というものなので、無視してもかまいません。
結果を見ると、すべて同じ値になっているはずです。
ですが、各成分に対応する周波数がありません。
出力の一番上を0番目としたn番目の周波数をF[n] [Hz]、窓の大きさをL[pts]として
F[n]=n*fs/L
と書けます。
よく見ると結果はすべて実数です。
連続時間領域でデルタ関数をフーリエ変換すると
δ~(f)=exp(-i*2πf*0)=cos(-2πf*0)+i*sin(-2πf*0)=1
すなわち、周波数にかかわらず実数で一定です。

DFTの結果はこれの定数倍であり、離散信号においてもフーリエ変換が成立することが分かります。
フーリエ変換の出力は必ず複素数なので、それに対応した出力となっています。
Excelでの複素数は「(実数)±(虚数)i」で表記します。
もちろん、Excelには複素数に対応した関数が用意されています。
Zoom 直交座標だと分かりにくいので、周波数成分の電力を示す「パワースペクトル」に直しましょう。
電力は電圧の2乗に比例するので、複素数の絶対値を2乗します。
Excelの関数で複素数の絶対値はIMABS()で求められます。

変換後のd[t]をd~[f] [Hz]、パワースペクトルを|d~[f]|^2 [W]と置くと
|d~[f]|^2=IMABS(d~[f])^2
と書けます。
Zoom 結果は定数なので面白くありませんが、グラフ化してみましょう。

ナイキストのサンプリング定理より、ナイキスト周波数以上のスペクトルは利用できません。
今回は関数の性質上確認できませんが、DFTのパワースペクトルはナイキスト周波数を軸に左右対称になります。
正確に言うと、複素スペクトルにおいて複素共役になっており、逆相、すなわち負周波数のスペクトルを示しています。
余談:
ちなみに極座標の位相を求めるにはIMARGUMENT()です。
もちろん角度の単位はラジアンです。
余談2:
フーリエ変換が可逆変換であることは周知のとおりです。
先ほどのパラメータ指定ウインドウで気づいたかもしれませんが、ExcelのDFTも可逆変換です。
気になる方は線スペクトルなど、時間領域に変換して確かめてみてください。

これで信号解析が一気に楽になります。
今後、離散時間領域から離散周波数領域へ変換するにはこの方法を用います。

もう少し調べてみましょう。
「パーセバルの定理」というものがあります。
電気的な解釈すると「フーリエ変換の前後で電力の総和は変化しない」となります。
式で書くと連続領域では
Σ(g(t))^2=Σ(g~(f))^2
DFTを用いた離散領域では
Σ(g[t])^2=1/L*Σ(g~[f])^2
となります。

Zoom まず、時間領域の電力の総和を求めます。
式は先ほどのとおりですが、Excelには二乗総和を求めるSUMSQUARE()関数があります。
窓の大きさと同じ分だけ選択しましょう。
値は見やすい適当な場所へ配置します。
Zoom 同様に、周波数領域にも適用します。
こちらはパワースペクトルで電力を求めたので、SUM()関数を使い総和を求めましょう。
Zoom 値には1か0しかないので、計算せずとも想像はつくでしょう。
時間領域は1 W、周波数領域は128 Wです。
今度は、窓の大きさを変えて今度は256 ptsでフーリエ変換してみましょう。
その前に、デルタ関数も256サンプル用意し、総和の範囲も256サンプルに変えます。
Zoom すると周波数領域の総和は256 Wになります。
引き続き、窓の大きさを変えてもいいのですが、もうお分かりでしょう。
周波数領域の電力の総和は、窓の大きさとの積になります。
したがって、パーセバルの定理は成立することが分かります。
戻る
inserted by FC2 system