Excelでn=8のFFTを実装してみた

FFTで高速にフーリエ変換できることは、なんとなく知っています。
でも、原理や有効性を肌で感じたことがありません。
そこで、Excelのワークシート関数のみで、n=8のFFTを実装してみました。

FFTの原理
n=8のバタフライ演算
先ずは基本波を準備する
FFTの実装
基本波を解析してみる
混合波を解析してみる
ナイキスト周波数
６倍波を解析してみる
基本波の整数倍でない場合
まとめ

FFTの原理

バタフライ演算などを含め、FFTを分かり易く解説しているサイトがあります。
ところどころタイプミスはあるものの、これ以上の解説記事を書くことは、僕にはできません。
原理を全くご存じない方は、先ずはこちらを一読ください。

うさぎでもわかる信号処理・制御工学　第14羽　高速フーリエ変換(FFT)

今回は、高速フーリエ変換(FFT)について、具体例を用いたアルゴリズムの動き方、バタフライ演算による演算の可視化、およびソースコード例を用いてわかりやすく解説しています。

n=8のバタフライ演算

n=8のバタフライ演算の図を示します。

これをワークシートに実装すれば、FFTができます。
f8_0やf2_7は、ワークシートで使う名前です。

先ずは基本波を準備する

いきなり難しい波形を解析しても何が何だか分からなくなります。
先ずは基本波を用意しましょう。

新たなシートを用意し、「元波形」とでも名前を付けます。

１行目はタイトル行にします。A1、B1、C1にぞれぞれ「n」「time」「f8」を入力します。

A列はインデックスです。A2からA9には0から7の数字を順に入れます。

FFTは１周期分の波形を解析します。
B列に時間を設定します。
分かり易いように、周期は1[s]としましょう。
1を8分割すると、増分は0.125です。
定数を入れても問題ありませんが、ここでは「=(ROW()-2)/8」とします。
さらに、B2:B9を選んで、左上のセル名の欄に「t」と入れてEnterすることで名前を付けます。

C列は解析したい波形を入れます。
Excelの三角関数が使う単位はradです。
周期が1[s]の正弦波は、1秒で2π[rad]となればいいので、「=SIN(2*PI()*@t」とします。
@tとするのは、先に名前を付けたt（B2:B9）の中で、同じ行の値を使うということになります。

実際の値とグラフを示します。

また、A列全体を選んで「n」、C列全体を選んで「_f8」という名前を付けておいてください。（f8という名前は、F8セルを指すことになるので、アンダスコアを付けておく）

FFTの実装

同じブックに新しいシートを追加します。
名前は「FFT」とでもします。

Excelは、複素数計算ができます。
小手試しに、 $w_8^1$ を表現してみます。

$\begin{align*}w_N&=e^{-\frac{2\pi}{N}i}\\w_8&=e^{-\frac{2\pi}{8}i}=e^{-\frac{\pi}{4}}=\cos\left(-\frac{\pi}{4}\right)+i\sin\left(-\frac{\pi}{4}\right)\end{align*}$

B11セルに、タイトルとしてw8と入力します。
C11セルに、ここで求めた複素数を「=COMPLEX(COS(-PI()/4, SIN(-PI()/4))」と入力します。

$w4_4^1$ 、 $w2_2^1$ についても同様にして、B12:C13を埋めます。
埋めたらC11、C12、C13に、それぞれ「_W8」「_W4」「_W2」と名前を付けておきます。

では、いよいよf8(n)～f1(n)を埋めましょう。

先ず、1行目のA～E列に、それぞれ、「n」「f8(n)」「f4(n)」「f2(n)」「f1(n)」と入力します。

A列の2行目から7行目は、0～7で埋めます。

B列のf8(n)は、「元波形」シートのC列をコピー＆ペーストしても良いですが、B2セルに「=COMPLEX(INDEX(_8, MATCH(A2, n)), 0)」とすることで、f8(0)の値を複素数として持ってこれます。f8(1)～f8(7)を持ってくるためには、B2セルをコピーして、B3:B9にペーストします。
B2からB9には、「f8_0」から「f8_7」と名前を付けておきます。面倒ですが、式の分かり易さには代えられません。

C列のf4(n)からいよいよバタフライ演算の出番です。

C2セルのf4_0は、f8_0とf8_4から線が来ているので、単純に複素数の和を取ればいいので「=IMSUM(f8_0, f8_4)」となります。

C6セルのf4_4は、やはり、f8_0とf8_4から線が来ています。でもf8_4の線の脇には-1が書かれているので、f8_4は引き算をする必要があります。さらに、交差した後に $w_8^0$ がかかっています。よって、差に $w_8^0$ がかかるので、C6セルは「=IMPRODUCT(IMSUB(f8_0,f8_4),IMPOWER(_W8,0))」となります。そのままですね。
バタフライ演算の図があれば、機械的に作業ができます。
面倒ですが、C2からC9には、「f4_0」から「f4_7」の名前を付けます。

f2(n)、f1(n)についても同様です。頑張りましょう。

f1(n)まで埋めることが出来たら、Bit reverseの番です。
必要性につきましては、先に紹介したサイトをご覧ください。

ともかく、周波数F(k)として並べるときには、Bit reverseしたときの値で並べます。例えば、f1(3)を考えてみます。
n=8なので、３ビットです。3を3ビットで表すと、「011」です。これをリバースすると「110」であり、10進に戻すと6です。なので、F5セルには6を入れています。
G列は、F列の順番でE列を並べ替えて、ある定数で割った結果です。つまり、F(6)=f1(3)/4となっています。（列名と周波数を表すF、関数のfが出てきて紛らわしいですね）

この割る数ですが、F(0)（G2セル）の場合は分かり易いです。
バタフライ演算をよく見ると、f1(0)は、f8(0)～f8(7)の和になっています。F(0)はDC（直流）成分ですので、平均値を取らなければなりません。よって、8で割ります。すなわち、F(0)=f1(0)/8です。また、F(4)も8で割ります。理由は後述します。

F(0)とF(4)以外の場合は、n/2つまり4で割ります。この理由は、よく分かりません。割ってみると丁度、振幅が良い値になります。ただ、以下の２つの理由のどちらか、又は両方だと思っています。詳しい方は教えてください。

１）フーリエ級数展開で、交流のスペクトルは直流の半分の値で割るから
「フーリエ級数展開の原理」のページの「直交とは」の節で示したように、直流値は積分値を2π、交流値は積分値をπで割ってスペクトルを求める必要があります

２）交流のスペクトルは、正のスペクトルと負のスペクトルに分けられるから
直流のスペクトルはF(0)に出てきます。しかし、交流のスペクトルは正負？に分けられます。例えば、基本波のスペクトルは、F(1)とF(7)で同じ大きさになります。２倍波のスペクトルはF(2)とF(6)、３倍波はF(3)とF(5)です。よって、そもそも半分に分けられているので、更に直流の場合の半分で割れば、丁度良い大きさになるという訳です

G列に対し、H列はその大きさ、I列は実部、J列は虚部となっています。

基本波を解析してみる

まずは、正弦波を解析してみます。 $f=\sin\left(2\pi t\right)$ です。

これがn=8の場合の数値とグラフです。
周期は1[s]、振幅も1で、n=8で１周期になっています。いわゆる基本波です。

FFTの演算過程はこのようになり、

周波数成分はこのように求められます。
基本波なので、n=1のところに振幅１のスペクトルが確認できます。
sin波なので、実数成分は無く、虚数成分が-1になっています。

なんか、虚数成分が負の値になることは、納得いきませんね。
ただ、これは、cosを基本に考えているためのようです。
すなわち、 $\sin\left(x\right)=\cos\left(x-\frac{\pi}{2}\right)$ なので、 $\cos x$ に対する位相差が負であることを示しているのだと思います。

なお、n=7のところにある振幅1のスペクトルは負の基本波成分（-1[Hz]）を意味します。

グラフにするとこのようになりますが、普通はDCから正の周波数領域（n=0～4?）を表示し、n=5～7は非表示にしているはずです。

混合波を解析してみる

次に、 $0.5+2\cos\left(2\pi t\right)+\cos\left(6\pi t+\pi/6\right)$ を解析してみます。

グラフはこんな感じ。見た目では、どんな周波数成分が入っているか、よくわかりません。

でも、FFTをすれば一発です。DC成分が0.5、基本波成分の振幅が２、３倍波の振幅が１であることがはっきりと分かります。
さらに、基本波は実数成分だけなので、cos成分のみで位相差が0であることが分ります。
３倍波の実数と虚数はそれぞれ0.86603と0.5なので、Excelで、=DEGREES(ATAN2(0.86603,0.5))を計算するとほぼ30°になります。30[°]=π/6[rad]なので、あらかじめ与えた位相差と等しくなります。

ナイキスト周波数

ナイキスト周波数は、サンプリング周波数の半分の値です。
今、1 [Hz]で8点をサンプリングしているので、サンプリング周波数は8 [Hz]です。
よって、サンプリング周波数は4[Hz]です。

では、 $f=\cos\left(4\times2\pi t\right)$ を解析してみます。

最大値と最小値を交互にサンプリングします。

先にn=7は、負の基本波を意味していると書きました。同様に、n=6, 5, 4は、それぞれ、負の2倍波、3倍波、4倍波を示します。
n=4は、正の4倍波であって、負の4倍波でもあります。つまり、２倍になっているので、これが $F\left(4\right)=f\left(1\right)/8$ のように、8で割った理由です。

ナイキスト周波数の正弦波、すなわち、 $f=\sin\left(4\times2\pi t\right)$ も見てみましょう。

これのFFT結果は自明だと思います。

丸め誤差はあるものの、恒等的に０です。

何が言いたいかというと、世間的に言われているサンプリング定理「アナログ信号をデジタル信号に変換する際に、元の信号を正確に再現するためには、信号の最大周波数の2倍以上のサンプリング周波数が必要であるという定理です。」という表現は、厳密には誤りということです。周波数が4[Hz]のsin波は、2倍の8[Hz]でサンプリングしても「正確に表現できません」。「最大周波数の２倍を超える周波数でサンプリングしなければいけません」

６倍波を解析してみる

６倍波も解析してみましょう。

滑らかな線でつないでいます。とても6 [Hz]には見えません。2 [Hz]といったところでしょうか。

スペクトラムもn=2に立っています。でも、確かにn=6にも立っています。n=6は、REALもIMAGも約0.70711なので、位相差はπ/4です。この値も合っています。

位相差はともかく、ナイキスト周波数を越えると、スペクトラムは、ナイキスト周波数で折り返したように出現します。この現象をエイリアシングと言います。
エイリアシングの観測を抑えるためには、ナイキスト周波数以下でカットオフするようなLPFを用いることが多いです。

基本波の整数倍でない場合

最後に、周波数が基本波の2.5倍の場合を見ておきます。

随分歪んでいるように見えます。

FFTの結果を示します。
n=2やn=3のスペクトルが大きいことは確かですが、これらが等しいわけではありません。
時間軸のサンプリング結果から高分解能で周波数のピークを得るためには、十分に速いサンプリング周波数で、多くのサンプルを取る必要があることが分ります。

まとめ

学生時代に「いつかはFFTについて再学習しよう」と思いながら、何十年が経ってしまいました。
自分で手を動かしましたので、原理が良く分かりました。
ナイキスト周波数やエイリアシングについても理解が深まりました。
途中経過も面白かったです。ただし、Excelの複素数の表示はもう少しうまく丸めて表示してもらいたいものです。

リンク