C言語でのFFTについて

質問

http://tsuyu.cocolog-nifty.com/blog/2007/03/publi.html
に掲載されているVBAのFFTプログラムをC言語に書き換えて
実行しているのですがうまくいきません。

どこが、間違っているか教えてください。
＝＝＝＝＝＝以下FFTのサブルーチンソースコード＝＝＝＝＝
void FFT(float Xr[], float Xi[])
{
i=0,j=0,k=0,l=0,m=0,n=0,p=0,q=0;
n=DN;
m = log10(n)/log10(2);
Table(c,s);
l=n,h=1;
for(g=1;g<=m;g++){
l/=2,k=0;
for(q=1;q<=h;q++){
p=0;
for(i=k;i<=l+k-1;i++){
j=i+l;
a=Xr[i]-Xr[j], b=Xi[i]-Xi[j];
Xr[i] = Xr[i] + Xr[j], Xi[i] = Xi[i] + Xi[j];
if(p==0){
Xr[j]=a,Xi[j]=b;
}else{
Xr[j] = a * c[p] + b * s[p], Xi[j] = b * c[p] - a * s[p];
}
p+=h;
}
k+=l+l;
}
h+=h;
}
j=n/2;
for(i=1;i<=n-1;i++){
k=n;
if(j<i){
//ビットリバース
swap(&Xr[i],&Xr[j]);
swap(&Xi[i],&Xi[j]);
}
k=k/2;
while(j>=k){
j=j-k;
k /=2;
}
j = j + k;
}
}

この投稿のマルチメディアは削除されているためご覧いただけません。

回答No.2

pyonmae · ベストアンサー

こんにちは。
変数の宣言や、未定義の関数、そもそも何がうまくいかないのか、など補足して頂きたい事がいっぱいあるのですが、私自身、FFTに強い関心がありますので、勝手に補完して検証してみました。
まず、未定義の部分を以下のように定義しました。
　変数 → 全てグローバル(Integer → int、Single → float)
　DN → 1024
　Table関数 → SinとCosのテーブルを作成
　swap関数 → 2つのfloat変数の内容を入れ替える
　FFT対象データ → 適当にSIN波を入力
これでプログラムを組んでみますと、コンパイルは通りました。
しかし実行させてみると途中で無限ループにはまったらしく、関数から帰ってきません。
調べてみますと、上記ソースの下から7行目の
＞while(j>=k){
で、jとkが共に0となり、ループから抜け出せなくなっていました。

では、元になったVBAはどうなのかと見てみますと、上記に相当する部分
＞Do While j >= k
において、jが0、kが0.5となり、比較結果がFalseとなりループを抜けていました。
さて、ここで大きな疑問。Integerで宣言されているはずのkがなぜ0.5になるのか？
VBAにおける宣言文は、以下です。
＞Dim g, h, i, j, k, l, m, n, o, p, q As Integer
一見、全て整数になりそうですが、いろいろ試してみると、整数はqだけで、それ以外g～pは実数になっていました。
つまり、「As Integer」文が有効なのは直前のqだけであり、それ以外はデフォルトの型(バリアントかな？)になってしまうようです。
VBAの作者の方がこれを狙ってコーディングをされたのか、たまたま動作しているのかは知る由もありませんが。。。
ですので、最後のループを実数で比較させるように変更してみたところ、とりあえずそれっぽく動作いたしました。
ただ、整数と実数が入り乱れた非常に見づらいプログラムになるので、出来れば内容を理解した上でもう少しスマートにコーディングしたいものです。

今回、FFTとVBAの宣言のクセという2つの収穫があり、質問者様に感謝を申し上げたく思います。
しかし、やはり無駄な苦労がありましたので、ご質問の際には全てのソースを載せていただきたく思います。

お礼

早速の回答ありがとうございます。
//FFTルーチンについて
お察しのとおりです。
入れ替えるところのwhile文で、jとｋが0になり無限ループにはまってしまうので、私はwhile文の中にjとｋがともに０になったときに抜けるようにしました。
そして実行結果を見たらVBAバージョンとC言語バージョンとでは大きく結果が違っていました。（もちろんVBAバージョンの結果は正しく表示されます。）なぜこうなるのか、どうすれば正しく表示できるのか、トライ＆エラーを繰り返して探っている最中です。
解決策がわかったら教えてください。
よろしくお願いいたします。
//ユーザー関数について
すいません。　ユーザー関数の部分を掲載するのを忘れていたので、
追加掲載します。参考にしてください。
Table関数
void Table(float Rx[],float Ry[]){
a = 0.00, b = 2.00 * PI / DN;
for(i=0;i<DN/2;i++){
Rx[i]=cos( a + i * b ),Ry[i]=sin( a + i * b );
}

}
swap関数
void swap(float *L,float *R){
float tmp;
tmp = *L, *L= *R, *R=tmp;

}

回答No.3

pyonmae

こんにちは。
私のところでは、それっぽく値が取れていますが・・・？
補足で挙げていただいた関数はVBAにも記述されているものであり、当方でも実装済みです。
問題なのは、対象データや、変数の宣言、最終的な演算（VBAで言う、"データ出力"以降）だと思います。
何故、小出しなのでしょうか。
とは言え、大きくは間違ってはいないと思いますので、これ以降の回答は不要な気がします。
がんばって下さい。

あと、先の回答で私が得意げにVBAの宣言について申し上げましたが、VBA版のブログをよく見ると、コメントでちゃんと言及されていました。
こちらはお恥ずかしいです。

お礼

どうやら、出力時に直流分以外のところで2倍していなかったこととサンプル数で割っていなかったことが原因だったみたいです。（お騒がせしました。）
上記の処理を施した後の実行結果を振幅値で表示したところ、VBAバージョンとC言語バージョンの結果が一致しました。しかし、ｄB表示すると双方で、約267dBの違いがありました。
この原因（振幅値表示が一致しているのにｄB値表示が異なること）は、いまもって謎です。

どうもありがとうございました。

回答No.1

Tacosan

「うまくいかない」とはどういう現象を指すのでしょうか?
コンパイルできる完全なプログラムとデータ並びに「期待した結果」と「実際に得られた結果」を出してください.

お礼

すいません。全体をのせることができないので、FFTルーチンとそれに使用しているユーザー関数を載せます。
void FFT(float Xr[], float Xi[])
{
i=0,j=0,k=0,l=0,m=0,n=0,p=0,q=0;
n=DN;
m = log10(n)/log10(2);
Table(c,s);
l=n,h=1;
printf("FFT Processing computing.......\n");
for(g=1;g<=m;g++){
l/=2,k=0;
for(q=1;q<=h;q++){
p=0;
for(i=k;i<=l+k-1;i++){
j=i+l;
a=Xr[i]-Xr[j], b=Xi[i]-Xi[j];
Xr[i] = Xr[i] + Xr[j], Xi[i] = Xi[i] + Xi[j];
if(p==0){
Xr[j]=a,Xi[j]=b;
}else{
Xr[j] = a * c[p] + b * s[p], Xi[j] = b * c[p] - a * s[p];
}
p+=h;
}
k+=l+l;
}
h+=h;
}
j=n/2;
for(i=1;i<=n-1;i++){
k=n;
if(j<i){
swap(&Xr[i],&Xr[j]);
swap(&Xi[i],&Xi[j]);
}
k/=2;
while(j>=k){
j=j-k,k /=2;
if((j==0)&&(k==0)){
break;
}
}
j += k;
}
printf("FFT Processing Complete!\n");
}
void Table(float Rx[],float Ry[]){
a = 0.00, b = 2.00 * PI / DN;
for(i=0;i<DN/2;i++){
Rx[i]=cos( a + i * b ),Ry[i]=sin( a + i * b );
}
}
void swap(float *L,float *R){
float tmp;
tmp = *L, *L= *R, *R=tmp;
}

補足

入力データは、０～1秒、1024サンプル、振幅１、５０Hzの正弦波です。
出力結果はここでは、表示できないので、画像で掲載しますので、ご検討ください。