WindowsでFFTWを使う
by K.I
2010/03/24
Index
- 周波数解析を行うためのFFT（高速フーリエ変換）ライブラリ FFTWを
- とりあえず、VC5を良く使うので、VC5で試してみた。
- VC5で動くんなら、VisualStudio等で使うのも可能だろう。
- インストールというより、fftwのDLLを使うためのLIBファイルの作成。
- 解凍して、fftw3というディレクトリ名に変更しておく。
- Windowsでコマンドプロンプトを起動、VCVARS32.BATを実行
lib /machine:i386 /def:libfftw3-3.def
lib /machine:i386 /def:libfftw3f-3.def
lib /machine:i386 /def:libfftw3l-3.def
以下のファイルが生成される。
libfftw3-3.lib （単精度用）
libfftw3f-3.lib （倍精度用）
libfftw3l-3.lib （ロング精度用）
- インストールで作成したlibファイルと、dllファイル、それにfftw3.hを、プロジェクトファイルにコピー1。
- あとは、fftw3.hをプロジェクトに追加しておく。
- 最終的にWindowsプログラムにしたいので、ちょっとだけ改造。WinMainに変更してみた。
#include <windows.h>
#include <stdio.h>
#include <math.h>
#include "fftw3.h"
#pragma comment(lib, "libfftw3-3.lib")
#define SIZE 128
#define RANGE 10
int WINAPI WinMain(HINSTANCE hInst, HINSTANCE hPrevInstance, LPSTR lpCmdLine, int nCmdShow)
{
char line[256];
// コンソールを生成
AllocConsole();
freopen("CONOUT$", "w", stdout);
freopen("CONIN$", "r", stdin);
int i;
fftw_plan plan;
fftw_complex *in_buf, *out_buf;
in_buf = (fftw_complex *)fftw_malloc(sizeof(fftw_complex)*SIZE);
out_buf = (fftw_complex *)fftw_malloc(sizeof(fftw_complex)*SIZE);
for(i = 0; i < SIZE; i++){
in_buf[i][0]=0.0;
in_buf[i][1]=0.0;
}
for(i = 0; i < RANGE; i++)
in_buf[i][0] = in_buf[SIZE - i - 1][0] = 1.0;
printf("#Input Function\n");
for (i = 0; i< SIZE; i++)
printf("%d %f %f\n", i, in_buf[i][0], in_buf[i][1]);
printf("\n\n");
plan = fftw_plan_dft_1d(SIZE, in_buf, out_buf, FFTW_FORWARD, FFTW_ESTIMATE);
fftw_execute(plan);
printf("#FFT by fftw\n");
for (i = 0; i< SIZE; i++)
printf("%d %f %f\n", i, out_buf[i][0] / sqrt(SIZE), out_buf[i][1] / sqrt(SIZE));
gets(line); // Returnキー入力待ち
return 0;
}
とりあえず、サンプルのお陰で、ここまで簡単に出来ました。感謝です！
- データだと、ちょっと分かりづらいので、CSVファイルを作成して、
FILE *fp;
fp = fopen("test.csv","w");
for (i = 0; i< SIZE; i++)
fprintf(fp,"%d,%f,%f\n", i, out_buf[i][0] / sqrt(SIZE), out_buf[i][1] / sqrt(SIZE));
fclose(fp);
- Excelはとにかく手軽にグラフに出来るので便利。これで普通は十分だろう。
- まぁ、如何にもExcelで描きましたって感じのグラフになるのが欠点かな。
- Graph出力ツールとしては、定番のgnuplotをWindowsで使ってみる。
- gnuplotのダウンロードページから、gp440win32.zipをダウンロード。
- これで、コマンドプロンプトから、以下のように起動出来る。
- wgnuplotは、ウィンドウを表示して実行する。
wgnuplot
- デフォルトでは、文字が小さくて見えないので、
- 右クリック→ChooseFont...で、フォントを設定
- 右クリック→Update xxxx.iniで、設定を更新しておく。
- プログラムからは、単純なgnuplotコマンドの方が使いやすい。
- pipe用にpgnuplotというのもあるらしい。
- 以下のように、plotコマンドをファイルにして、systemでシェル実行すればプロットされる。
fp = fopen ("test.gp","w");
fprintf(fp,"plot \"test.csv\" using 1:2 w l,\"test.csv\" using 1:3 w l\n");
fclose(fp);
system("gnuplot -persist test.gp");
本当は、もうちょっとちゃんと設定した方が良いけど、これでもそれなりに出力される。
- とりあえず、これで試しにグラフ出力するのに使えそうだ。
plot [:][:] 'test.dat' using ($1/1000):4 with impulses linewidth 1
- これは、'test.dat'ファイルのデータをプロットする。
- [:][:]で、X,Yの表示範囲を指定する。この場合は、指定していないので、データの範囲によって自動的に設定される。
- usingで、表示カラムを指定する。
- カラムnを$nと指定すると、計算値を表示することも出来る。この場合、計算式は（）で括る必要がある。
- FFTWによる結果は、カラム１がHz単位なので、($1/1000)とすることで、kHz単位とする。
1本当は、サンプルのページみたいにちゃんとした場所に入れるべきなんだけど、とりあえず。
[top]
- 参考にさせて頂いたサンプルプログラムはちゃんと理解していなかったので、改めて使い方を考える。
- 元々のデータは時間の情報を含まないものなので、サンプリング周波数が基準になるはずだ。
- でも、サンプリング周波数と、FFTWの出力の関係がイマイチ分からないので実験してみる。
- 単純なサイン波をFFTWに掛けて実際にどうなるのか見てみよう。
- １秒分のサイン波を生成するプログラムがこれ
void singen(WORD data[], int sample, double f)
{
double pi = 3.1415926536;
int max = 65535;
double d;
double out;
int t;
d = 2*pi*f/sample;
for(t=0; t<sample ;t++) {
out = sin(t*d);
data[t] = (WORD)((out+1)*max/2);
}
}
- 1kHzのサイン波、サンプリング周波数1MHzで１秒分、つまり1000000個のデータを生成して、
- gnuplotで、0〜1000の範囲だけ表示してみる。
#define N 1000000
WORD data[N];
int WINAPI WinMain(HINSTANCE hInst, HINSTANCE hPrevInstance, LPSTR lpCmdLine, int nCmdShow)
{
int i;
FILE *fp;
singen(data,N,1000);
fp = fopen("test.csv","w");
for (i = 0; i< N; i++)
fprintf(fp,"%d\n",data[i]);
fclose(fp);
fp = fopen ("test.gp","w");
fprintf(fp, "plot [0:1000] \"test.csv\"\n");
fclose(fp);
system("gnuplot -persist test.gp");
return 0;
}
- FFTWで、１秒分のデータのDFTを掛けて、結果をファイルに入れるプログラムを書いてみた。
void fft_execute(WORD data[], int sample, char *fname) {
FILE *fp;
int i;
double re,im,mag,ang;
fftw_complex *in, *out;
fftw_plan p;
//sample個のFFTW用の複素数配列を確保する
in = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*sample);
out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*sample);
// 複素数配列に入力データをセット
for(i=0; i<sample; i++) { // t=0〜1secのデータ
in[i][0]=data[i]; // 実数部
in[i][1]=0.0; // 虚数部
}
// １次元のDFTを実行
p = fftw_plan_dft_1d(sample,in,out,FFTW_FORWARD,FFTW_ESTIMATE);
fftw_execute(p);
// 結果を出力する
fp = fopen(fname,"w");
for(i=0; i<sample; i++) { // f=0〜1MHz（サンプリング周波数）のデータ
re = out[i][0]; // 複素数の実数部
im = out[i][1]; // 複素数の虚数部
mag = sqrt(re*re + im*im); // 大きさを計算
ang = atan2(im,re); // 位相角を計算
fprintf(fp,"%d,%.5lf,%.5lf,%.5lf,%.5lf\n",i,re,im,mag,ang);
}
// 後始末
fftw_destroy_plan(p);
fftw_free(in);
fftw_free(out);
}
- Toshi1960さんの指摘で、位相角の計算式がimとreが逆になっていたのを修正しました。（160920）
- 結局のところ、FFTWの出力は、
- サンプリング周波数に比例した周波数分解能で、複素数出力される様だ。
- 普通は、周波数成分の大きさが必要なので、以下のように求める。
- で、DFTを掛けてから、gnuplotでグラフ出力するプログラムを書いてみる。
#define N 1000000
WORD data[N];
int WINAPI WinMain(HINSTANCE hInst, HINSTANCE hPrevInstance, LPSTR lpCmdLine, int nCmdShow)
{
FILE *fp;
singen(data,N,1000);
fft_execute(data,N,"test.csv");
fp = fopen ("test.gp","w");
fprintf(fp, "set datafile separator ','\n"
"set logscale x\n"
"set logscale y\n"
"plot [:][:] 'test.csv' using ($1/1000):4 with impulses linewidth 1\n"
);
fclose(fp);
system("gnuplot -persist test.gp");
return 0;
}
X軸は、kHz単位になるように1/1000スケールで出力した。
- また、誤差等も見たいので、X,Y軸ともLOGスケールにした。
- 当然、1kHzに強い信号があるが、それ以外にも矩形波なので、1kHz刻みに高調波が含まれている。
- そして、スペクトルが漏れたように、高調波が裾野を持って広がっている。
- 1kHzの周期の整数倍（1秒）の窓なので、サイン波に対する窓の影響は無いはずだが。。。。
- 雑音を含む原因を、いろいろ考えてみる。
- データが、16bitで量子化されている
- 1秒で切り取った（方形窓）データである
- 計算誤差
- 量子化の誤差は当然あるだろう。
- しかし量子化の誤差も1kHzの周期を持つはずなので、窓の影響は無い。（と思う）
- でも1kHzの周期を持たない雑音があれば、方形窓の影響を受けるかもしれない。
- 計算誤差は小さいが、1kHzの周期ではなく、蓄積する可能性がある。
- もしかすると、そのために量子化誤差も1kHz周期で無くなっているのかも。
- 試しに、サイン波を１周期の繰返しになるようにプログラムを変えてみた。
- 但し、これはサンプル数が周波数の整数倍の場合のみしか有効ではない。
void singen(WORD data[], int sample, double f)
{
double pi = 3.1415926536;
int max = 65535;
double d;
double out;
int t=0,i,j;
int m,n;
d = 2*pi*f/sample;
n = sample/f; // １周期当たりのサンプル数
m = f; // 繰返し数（つまり周波数）
for (i=0; i<m ;i++) {
for(j=0; j<n ;j++) {
out = sin(j*d);
data[t++] = (WORD)((out+1)*max/2);
}
}
}
→見た目はずいぶん変わった
- 小さな誤差であっても、窓の周期と異なっていれば、窓の影響を受けることが良く分かった。
- この例では計算誤差なので、非常に小さい（10^-8）無視できる程度の値だったが、実際のデータでDFTを掛ける時は注意が必要だ。
- サンプル数、つまり方形窓の大きさMを変えてみる。
- サンプリング周波数Nとの関係で、FFT出力のスケールをkHz単位になるようにした。
#define N 1000000 // サンプリング周波数
#define M 10000 // サンプル数（方形窓の大きさ）
WORD data[N];
int WINAPI WinMain(HINSTANCE hInst, HINSTANCE hPrevInstance, LPSTR lpCmdLine, int nCmdShow)
{
int i;
FILE *fp;
double f = 1000;
singen(data,N,f);
fft_execute(data,M,"test.csv");
fp = fopen ("test.gp","w");
fprintf(fp, "set datafile separator ','\nset logscale x\n"
"set logscale y\n"
"plot [:][:] 'test.csv' using ($1*%f):4 with impulses linewidth 1\n"
,(double)N/M/1000.
);
fclose(fp);
system("gnuplot -persist test.gp");
return 0;
}
・サンプル数、つまり窓の大きさMを変えることで、
- 主成分の1kHzは別として、高調波とかには変化があると思ったが、
- 予想に反して、高調波の分布には変化が見られない。
- 計算で生成したデータなので、窓の大きさが1kHz周期の倍数である限り影響が無いらしい。
- でも、1kHzの大きさは変化がある。→サンプル数が多い程、大きな値になっている。
- 今度は、サンプル数は固定で、サンプリング周波数を変えてみる。
→1kHzと3kHzにピークがある
→1kHzと9kHzにピークがある
→1kHzのピークのみ。99kHzにピークが無い？
- サンプリング周波数が変わっても、主成分の1kHzの大きさは変わらないが、高調波の分布はけっこう変わる。
- サンプリング周波数Nが高い方が多いように見えるが、結局はNが高い方が高調波は小さい。
- それから、サンプリング周波数の半分で折り返した分布になるので、
- N=4kHzでは3kHz、N=10kHzでは9kHz、N=100kHzでは99kHzに1kHzの折り返しが出るはず。
- サンプリング周波数4kHz,10kHzの時は折り返しが見られるが、100kHzでは何故か折り返しが見られない。
→1kHzと99kHzにピーク
- サンプル数を10000個と大きくすると、99kHzの折り返しが見えるようになった。
- 窓が小さいと、誤差で範囲ギリギリの所は見えなかったりするのかな？
- これまでのデータは、１周期分のデータのみFFTWに掛ける場合は、周波数分だけ同じデータを並べて１秒分のデータにして掛けていた。
- しかし１周期分のデータがあれば、同様に求められるはずだ。
- 短いデータをFFTWに掛けた場合、結果のスケールが変わってしまう。
- 例えば、0.2us毎の500個のデータが丁度、１周期分のデータだったとする。
- この場合、全部で100us分のデータとなる。1/10000秒のデータになるが、
- これをFFTWに掛けると、周波数レンジが1/10000になってしまう。
- つまり、以下のように設定すれば、１秒分のデータを入れたのと同じスケールになる。
plot [:][:] 'test.dat' using ($1/1000*10000):4 with impulses linewidth 1
- もう少し一般的に言えば、FFTを掛けるデータがn秒分とした時、以下のようにする。
- FFTを掛けるデータが１周期（或は周期の整数倍）分のデータではない場合は、方形窓の影響が出る。
- この場合は、窓関数を掛けて、窓の影響を緩和するようにした方が良い。
- 実際の波形データにFFTを掛ける場合、波形データの長さで繰返した波形になる。
- 波形データの長さが、含まれる周波数の周期に一致していれば問題ないが、異なる場合はそれが余計な周波数成分に見える。
- 窓関数は、スパッとデータを切り取るのではなく、切り取ったエッジを丸めて、この影響を少なくするような感じ。
- 単純にデータを取り込んだ場合、窓関数は単純な方形波になるので、
- 理想的には、１つの周波数が含まれる場合、１本の幅の無いスペクトルになって欲しいが、
- データを切り取ることで、スペクトルが漏れたような裾野を持った周波数分布になってしまう。
- 方形窓は、周波数分解能は高いが、裾野の部分（サイドローブ）が大きい。
- いずれにせよ、サイドローブを持ったスペクトルになっている場合は、窓の大きさと違う周期のデータを含んでいるということになる。
- 窓関数に関しては、測定器玉手箱の 窓関数の解説が分かりやすい。
[top]
- FFTWは、ちゃんと使えることが分かったが、
- GPLライセンスなので、配布時に著作権表示を残すだけでなく、作ったプログラムを複製・改変できるように公開しなければならない2。
- 実験的なものまで公開する必要はないと思うが、ちょっと配布したい場合に、
- 大したことないプログラムをわざわざ公開するのは、ちょっと面倒だ。
- FFTWと使用方法が互換性のある FFTSSというライブラリがあるらしい。
- ソースコードのfftss-3.0-20071031.zipをダウンロードして、win32フォルダにVisualStudioのslnファイルがあったので、
- コンパイルしたら、VS2008Expressでは問題なくコンパイル出来た。→実際には動かしてないけど。
- でもVC5で、試しにプロジェクトを作ってやってみたけど、ちょっとコンパイル出来なかった。
- VisualStudioで出来れば、良いんだろうか。（未確認）
- FFTSSの場合、MITライセンスなので、著作権表示は必要だが、それで配布可能なので、こちらの方が良いかもしれない。
2良くわからないが、そういうことなんだと思う。
[top]
[プログラムの部屋に戻る]
→1kHzと3kHzにピークがある
→1kHzと9kHzにピークがある
→1kHzのピークのみ。99kHzにピークが無い？
→1kHzと99kHzにピーク