ちゃんと動いた！

ＦＦＴ計算をＥＸＣＥＬでやっちまおうっていうお話しの続き。
ようやく思ったとおりの結果が出るように。（￣ー￣）

fft_dif_excel.xls
（↑このブログサービス、画像以外のファイルも色々貼れる
　ようになったみたいなので、早速ペタリ）

内容的には、8サンプルの周波数間引きＦＦＴを計算するもの。


入力は左側の赤で囲った2箇所。一つはサンプリングレート。
この例では1000sps…ＡＤＣを毎秒1000回行ったと仮定。
もう一つは入力信号としてS0～S7の8つ。基本的には
実数部だけ入力して使うんだけど、虚数も入力可。
そこだけ入力すればあとは自動計算＆グラフも自動描画。

出力的には、G0～G7の8個で、そのうちG0～G4が
各周波数成分として利用する部分。G0が直流成分。
ちなみにビット入れ替えの処理も済ませてあるので、
上から順にG0～G7の順で読めばＯＫ。
まぁ、それよりもグラフに自動反映されるので、グラフを
見るのが一番解りやすいかな…。

例の本読んで理解できたことをそのままワークシート
に落としただけなんだけど、単に計算式をワークシートに
落としただけじゃなくて、リクツを元に計算式を作り直して、
それを実際に動かしてみたら理論通り動いたってところが
ミソ。
しかも、このセルに埋め込んである計算式は、理論上
サンプル数を幾ら増やしても対応できる仕組み。

そう。これを思い切り拡張して、アセンブラでフルクラッチ
する目論見なのだ。これでとりあえず目処は付いたかな…

そう。ただ周波数成分を計算するところでちょっと…。
実数成分と虚数成分それぞれの2乗の和の平方根を
求めないといけないんだけど、アセンブラで平方根がなぁ…

ってことで、とりあえずニュートン法を適用したらどんな
風になるかな…ってことで試算を。

ニュートン法は漸化式で誤差を詰めていく方法だから、
初期条件の与え方が肝になるわけだけど、例によって
表計算のワークシートを旨く利用して試算してみた…。

処理上現実的な方法ってことで初期条件の与え方を
考えてみてから、それを元に誤差およそ２～３％程度
まで近づけるには漸化式を2回計算するだけで十分っぽい。
スバラシイ。
3回も計算すればもはや１％未満になるみたい。うーん、
現実レベルと考えて良さそう。
あとはそれをアセンブラで高速に処理させることができる
かどうか…だな。

だいぶ見えてきたので、手が空いてきたら徐々に
プログラムに落とし込んでいきたいところだな。

32点とするか、64点とするか、128点とするか…
そこが最大の課題になりそうだから、一度ウクレレ弦
を実際に使ってみて、ｐｉｃｏ　ｓｃｏｐｅでＦＦＴを掛けてみて、
それでどんな風かを眺めておきたいところだな…

表計算のワークシートでＦＦＴ…をちょっと進める。

今日も例によってやり直しのための信号数学を元に
表計算ソフトで８点ＦＦＴ計算のワークシートの手を
進めてみる。

とりあえず、落書き帳で最小単位のバタフライ演算の式
（２入力）を「複素数ベクトル表示」→「実数部＋虚数部」
に分解してみる。
ポイントになるのは回転因子のｋ乗の展開。算数苦手な
おいらは、これをワークシート上に落とし込むところで
しばしぼう然。
オイラーの公式を漠然と覚えてるだけじゃダメなんだな…

回転といえばこのオイラーの公式を回転因子に変形して、
それの階乗（ｋ乗）を掛けてあげればいいわけ。
　　e^(jθ)=(cosθ+jsinθ)
で右辺は実数項と虚数項に分かれてて処理上都合いいので、
このθに-2π/Nを代入してから右辺まるごとｋ乗すれば
いいのかな、と思ったんだけど…しばし思考停止。

このままｋ乗しちゃうと、二項定理で愚直に展開しなきゃ
いけなくなっちゃう。それじゃぁ単純なロジックで
処理できないし、できたとしても時間が掛かりすぎ。
ってことはこれはハズレだな。


うーん…、と改めてオイラーの公式を眺めなおす。
おぉ。左辺だ！　左辺をｋ乗すれば…
　　(e^(-j2π/N))^k = e^(-j2πk/N)
ってなって、その結果右辺は…
　　　　= cos(2πk/N)-jsin(2πk/N)
ってことになるのか。これなら二項定理は不要だな。よし。

改めて、２入力のＤＦＴ計算を落書き帳で筆算してみる。
…式はできた。あとはこれを愚直にワークシートに展開。
８入力の周波数間引きＦＦＴ処理として式を設定。

とりあえず完成したので、適当に入力データを設定してみる。


…うーん、出力がなんだか微妙におかしい…。
全セルを見直してみる。設計通りだ。なんだろう？

しばし悩む。うーん、どこかおかしいところがあるはずだ。
もう１回やり直しのための信号数学をパラパラ読み直す。


あ、ベクトル表記の式自体が間違えたまま筆算してた！
しくしく…。（ＴへＴ）設計不良。
下側は、足してから回転因子（ベクトル）を掛けるんだな。
ここが間違えてたよ。

またあとで落書き帳に向かって、筆算からやり直し。
そこだけやり直せば、まぁ今度は動くだろう…。
とりあえず半歩前進…と前向きに考えておこう。