てつ缶

2003-11　

　FFT で振幅スペクトラムとかパワースペクトラムくらいは求めたい。なぜか特定用途向けのスペルトラム分析ソフトはあるのに、各種方面に応用できて、しかも信頼できるソフトウェアは少ないように思います。

ずぼらな「てつ缶」ですから、Web で探しました。浮動小数点 FFT です。

FFT (高速フーリエ・コサイン・サイン変換)の概略と設計法
大浦拓哉先生、どうもありがとうございます。

これを使います。まずは cdft 複素 FFT です。

とりあえずは cdft() の動作テスト (ソースプログラム)

数の数え方でつまづきました。n というのはデータ数？ double a[] の数？ともかく動くようになって、1024 点 (doubleデータの数は 2048) の fast DFT ができました。

元データとして、cos(3x) / 3 + cos(5x) / 5 + cos(7x) / 7 を与えると、

n 1024 Re         Im       sqrt(Re^2 + Im^2) (正規化あり)
 0:   -0.000000   0.000000   0.000000
 1:   -0.000000   0.000000   0.000000
 2:   -0.000000   0.000000   0.000000
 3:  170.666667   0.000000   0.333333　＜ |X[3]| は 1/3
 4:   -0.000000   0.000000   0.000000
 5:  102.400000   0.000000   0.200000　＜ |X[5]| は 1/5
 6:   -0.000000  -0.000000   0.000000
 7:   73.142857   0.000000   0.142857　＜ |X[7]| は 1/7
 8:    0.000000  -0.000000   0.000000
　(後略)

スペクトラムアナライザー風の表示画面を作ってみました。

ソースプログラム　環境は Cygwin です。

(fft/sample1/fft4g.o は fft/sample1/Makefile で作ります。)

FFT のデータポイント数 1024 以外はちょっと手抜きしています。 7 種類のテストデータを FFT して、グラフィックに表示させるところまで進みました。

数字キーを押すと、それに対応したテストデータを作成し、元データと FFT 後のデータを表示します。処理に関連した情報がリストボックスに表示されます。
```
    0 : 0 にクリアー
    1 : Re[0] に振幅 51.2 のインパルス
    2 : 実数部は 1.0 * cos(2x)、虚数部は 0
    3 : 実数部 cos(3x)/3 + sin(254x)/10 + sin(510x)/10
        虚数部 -cos(5x)/5 + sin(7x)/7
    4 : 3 と同様、虚数部に振幅 0.5 のランダムノイズを加算
    5 : 矩形波。 幅は実行ごとに増加します。
    6 : 実数部に cos(1x) を入れて、振幅の目盛りに合うことを確認します。
    7 : ランダムノイズ (1)。 #define した RND() を使います。
    8 : ランダムノイズ (2)。 drand48() を使います。
```
目で見た限りでは RND() と drand48() のノイズ品質の区別はつきませんでした。きっと周期のランダム性に違いがあるのでしょう。
はて、応用は？

ちょっと変わった応用を考えました。工業計測の分野などでは、比較的ゆっくり変動するノイズにまみれた信号を測定する必要があります。グランドの絶縁も必須ですし、計測用の A/D 変換器などを使用したかなり大掛かりなものになりがちです。ちょっと手軽にやってみる方法はないものか。そこで、V/F 変換器 → 絶縁用低周波トランス → PC オーディオ入力 → 周波数分析をして入力電圧に換算する、という流れを思いつきました。 (もし 10 年前だったら、実にあほらしく思えた方法です。)

V/F 変換器とは、入力されたアナログ電圧を振幅が一定の周波数信号 (多くは矩形波) に変換するもので、安価なモノリシック IC としても販売されています。出力周波数は入力アナログ電圧の一次関数で表されます。この部分は電子工作で可能です。 PC オーディオ入力と周波数分析以降がソフトウェアの仕事です。

まず V/F 変換器の出力周波数範囲を決めます。応答速度の面から最低 1kHz くらいは必要でしょう。上限は PC オーディオの性能によって制限されます。周波数範囲が広いほうが有利ですから、きりのよいところで 1kHz ～ 15kHz (中心 8kHz) とか 2kHz ～ 16kHz (中心 9kHz) あたりが妥当だと思います。
```
サンプリング周波数 44.1kHz として、FFT データ点数と入力計測時間の関係は次のようになります。 

　　FFT 点数　　計測時間 (fs = 44.1kHz)
　───────────────────
　　　1024　　　　23.2ms
　　　4096　　　　92.9ms
　　 16384　　　　372ms
　　 65536　　　　1.49s
```
最終的なデータ出力周期と計測時間には、データ出力周期 ≧ 計測時間、という関係があります。例えば FFT 点数を 4096 とした場合、データ出力周期は約 1/10 秒になり、入力換算電圧分解能はフルスケールの約 1/2000 になります。

ちょっとおそまつですかね～。スペクトラムの広がりを計算して、入力電圧換算の分解能を高めることはできるでしょう。さらには、時間軸上の補間をして見かけ上の出力データレートを上げることはできます。 V/F 変換器の出力周波数の変化は実質的に連続 (*1) であって、PC オーディオ入力のサンプリング周波数とは非同期ですから、補間する意味は充分あります。ここはデジタルフィルターの出番でしょう。 (*1) 出力信号が内部クロックに同期しているタイプの V/F 変換器もありますが、このクロック周波数は PC オーディオのサンプリング周波数よりずっと高く、同期もしていません。

そうでした。 FFT の前に窓関数を乗じる演算を忘れてはいけません。これを怠るとスペクトラムが不必要に広がってしまいます。

計測周期は、例えば 1 秒に何回とか、正確に決める必要があります。多くの場合、測定時刻との関連付けも必要でしょう。

V/F 変換器の入力部分は純然たるアナログ回路ですから、温度ドリフトや経時変化があります。これを補正する方策も考えてなくてはなりません。
この続きは「や缶」にあります。
　や缶　FFT を利用した電圧計測

(つづく)

　や缶　　　　FFT を利用した電圧計測を試しました。

　ぶりき缶　「や缶」のフロントエンドになる V/F 変換器の部分などを考えました。

　あるみ缶　 GNU MP を使って多桁の演算をしてみました。

　river_r@d1.dion.ne.jp