フーリエ変換をして、周波数成分を取り出したい

1次元フーリエ変換を理解中の者です。

サンプリングした電圧のデータ列を分析して、周波数成分を取り出したいのですが、
フーリエ変換の概念はわかります。
やりたいのは、よく見かける、横軸周波数、縦軸電圧を作りたいです。

探したサンプルプログラムを実行して出力を得ました
実は一般的なソースらしいのですが、意味はまだ理解していません


1Khzでサンプリングした電圧のデータ列
int samp[20] = {0, 0, 0, 0, 0, 0, 0, 100, 100, 100, 100, 100, 100, 0, 0, 0, 0, 0, 0, 0};

結果
0 0 0 0 0 0 0 100 100 100 100 100 100 0 0 0 0 0 0 0
0 12 30 15 -42 -99 -80 60 292 510 600 510 292 60 -80 -100 -42 15 30 12
0 0 0 0 0 0 0 100 100 100 100 100 100 0 0 0 0 0 0 0

ソースは最後に添付します

お聞きしたいのは、この結果は、どう見たら良いのでしょうか
欲しい結果は、
1kHz = 〇v
10kHz = 〇v
．．．


ソース
int test_1_main()
{
double re[20], im[20];
int i;
FILE *fo;

char *outfile="test.txt";

if ((fo = fopen(outfile, "w")) == NULL) {
printf("file open error!!\n");
return (-1);
}

for(i=0; i<20; i++){
re[i] = samp[i];
im[i] = 0;
}

for(i=0; i<20; i++){
fprintf(fo, "%d ", (int)re[i]);
}

dft_swap(re, im, 20);
fprintf(fo, "\n");
dft_idft(re, im, 20, DFT);
dft_swap(re, im, 20);

for(i=0; i<20; i++){
fprintf(fo, "%d ", (int)re[i]);
}
fprintf(fo, "\n");

dft_swap(re, im, 20);
dft_idft(re, im, 20, IDFT);
dft_swap(re, im, 20);

for(i=0; i<20; i++){
fprintf(fo, "%d ", (int)re[i]);
}
fprintf(fo, "\n");

fclose(fo);

return 0;
}

call関数は書けないので抜粋以下
int dft_idft(double *re, double *im, int num, int flag)
{
略
for(i=0; i<num; i++){
temp_re[i] = temp_im[i] = 0.0;
}

if(flag==IDFT)coefficient=num;
else coefficient=1;

for(i=0; i<num; i++){
for(j=0; j<num; j++){
temp_re[i] += re[j]*cos(2*PI*i*j/num) + flag*im[j]*sin(2*PI*i*j/num);
temp_im[i] += -flag*re[j]*sin(2*PI*i*j/num) + im[j]*cos(2*PI*i*j/num);
}
temp_re[i] /= coefficient;
temp_im[i] /= coefficient;
}

for(i=0; i<num; i++){
re[i] = temp_re[i];
im[i] = temp_im[i];
}
略

}

No.3 ベストアンサー 回答者： ICE_FALCON 回答日時： 2011/04/28 01:38

#1です。

> 結果は、（元データ単位がVだとして）
> 0Hz　900v　
> 50Hz　1332.5v
> ・・・
> でしょうか

違います。
0Hzは30V(3列目)です。

pow=(re/20)^2+(im/20)　
は
pow=(re/20)^2+(im/20)^2
の間違いでした。
オイラーの公式知っていますか？
reはreal
imはimaginary
を示しています。
re*cos(θ)+i*im*sin(θ)
のパワーがpowです。
大きさはパワーの平方根でmagです。

”折り返し部分を纏めて”は
08行目のパワーと10行目のパワーは同じ周波数のパワーを示しているので
足しているだけです。
((510/20)^2+(-80/20)^2)*2=1332.5
これはパワーなので大きさに直すと
√1332.5=36.5Vです。

最大100Vの信号なので900Vなんて出るわけありません。

また、分析周波数単位という言葉はわかりませんが、
今回のように
1kHzサンプリングで20サンプルの系列をFFTで20個の周波数の系列に変換した場合、
変換後の周波数分解能は50Hzです。

この回答へのお礼

再度有難うございます
初心者ですみません
内容理解できました

お礼日時：2011/04/28 22:17

No.2 回答者： ICE_FALCON 回答日時： 2011/04/27 02:18

#1です。

スペースが全部消えてしまった・・・・。
全角スペースで張り直し。

　　　　re　　　im　　　re/20　　　im/20
00　　　12　　　-80　　　0.6　　　-4
01　　　30　　　-95　　　1.5　　　-4.75
02　　　15　　　-30　　　0.75　　　-1.5
03　　　-42　　　58　　　-2.1　　　2.9
04　　　-99　　　100　　-4.95　　　5
05　　　-80　　　58　　　-4　　　　2.9
06　　　60　　　-30　　　3　　　　-1.5
07　　　292　　　-95　　　14.6　　　-4.75
08　　　510　　　-80　　　25.5　　　-4
09　　　600　　　0　　　　30　　　　0
10　　　510　　　-80　　　25.5　　　-4
11　　　292　　　-95　　　14.6　　　-4.75
12　　　60　　　-30　　　　3　　　　-1.5
13　　　-80　　　58　　　　-4　　　　2.9
14　　　-99　　　100　　　-4.95　　　5
15　　　-42　　　58　　　-2.1　　　　2.9
16　　　15　　　-30　　　0.75　　　　-1.5
17　　　30　　　-95　　　1.5　　　　-4.75
18　　　12　　　-80　　　0.6　　　　-4
19　　　0　　　　0　　　　0　　　　　0

折り返し部分のパワーはまとめる。
　　　freq[Hz]　　　pow=(re/20)^2+(im/20)　　　mag=sqrt(pow)
00　　　　　　
01　　　　　　　　　　　　
02　　　　　　　　　　　　
03　　　　　　　　　　　　
04　　　　　　　　　　　　
05　　　　　　　　　　　　
06　　　　　　　　　　　　
07　　　　　　　　　　　　
08　　　　　　　　　　　　
09　　　0　　　　900　　　　30　
10　　　50　　　1332.5　　　36.5
11　　　100　　　471.45　　　21.71
12　　　150　　　22.5　　　　4.74
13　　　200　　　48.82　　　6.99
14　　　250　　　99.01　　　9.95
15　　　300　　　25.64　　　5.06
16　　　350　　　5.63　　　　2.37　
17　　　400　　　49.63　　　7.04
18　　　450　　　32.72　　　5.72
19　　　500　　　0　　　　　0

この回答への補足

回答有難うございます

結果は、（元データ単位がVだとして）
0Hz　900v　
50Hz　1332.5v
・・・

でしょうか
pow=(re/20)^2+(im/20)　から
00　　　12　　　-80　　　0.6　　　-4　のとき
0.6×0.6＋(-4)=4.36
ここから”折り返し部分を纏めて”の何かをすると、900になるのでしょうか

また、サンプリング単位50Hzが分析周波数単位と考えてよろしいのでしょうか

すみません、初心者で、よろしかったら教えて下さい

補足日時：2011/04/27 17:17

No.1 回答者： ICE_FALCON 回答日時： 2011/04/27 02:13

ソースが解らんところがあるから何とも言えんが・・・。

予想で言うと、re,imを全部20で割って、
更にそれを二乗和取ってパワーを求めて、
パワーの折り返し部分を纏めて、
平方根を取ったらＶででると思います。
周波数は1kHzサンプリングで20個のサンプルなので、
50Hz分解能ですね。
下記のような感じ・・見難いかな？

reimre/20im/20
012-800.6-4
130-951.5-4.75
215-300.75-1.5
3-4258-2.12.9
4-99100-4.955
5-8058-42.9
660-303-1.5
7292-9514.6-4.75
8510-8025.5-4
96000300
10510-8025.5-4
11292-9514.6-4.75
1260-303-1.5
13-8058-42.9
14-99100-4.955
15-4258-2.12.9
1615-300.75-1.5
1730-951.5-4.75
1812-800.6-4
190000

折り返し部分のpowはまとめる。
freq[Hz]pow=(re/20)^2+(im/20)mag=sqrt(pow)
0
1
2
3
4
5
6
7
8
9090030
10501332.536.5
11100471.4521.71
1215022.54.74
1320048.826.99
1425099.019.95
1530025.645.06
163505.632.37
1740049.637.04
1845032.725.72
1950000