C言語でのFFT(構造体とポインタの使用)
今、VC++2008 Express Editionを使用して先日、作成したFFTのプログラムを構造体とポインタを使用して書き換えました。
ビルド＆コンパイルは通ったのですが、実行するとエラーが起きます。何がおかしいのでしょうか。
ソースコードを掲載しますので、教えてください。
//main.c
#include <stdio.h>
#include "stdlib.h"
#include "fft.h"
#include <math.h>
#include "data.h"
#define PI 3.1415926
FILE *ifp, *ofp;
struct DATA f;
double dT, dF;
double *amp;
int main (void)
{
//波形の生成と配列入力
for(i=0;i<1024;i++){
dumx[i]=(double)i/1024;
dumy[i]=sin(2.0*PI*50*(double)i/1024);
}
DN=1024;
//FFT用データ配列に配列のコピー
time = (double*)calloc(DN,sizeof(double));
f.Re =(double*)calloc(DN,sizeof(double));
f.Im =(double*)calloc(DN,sizeof(double));
for(i=0;i<DN;i++){
time[i]=dumx[i];
f.Re[i]=dumy[i];
f.Im[i]=0.00;
}
//サンプリング周期
dT = time[1];
//サンプリング周波数
dF = 1.0/dT;
//FFT
FFT(f);
//大きさを出力するための配列の確保
amp=(double*)calloc(DN,sizeof(double));
for(i=0;i<DN;i++){
amp[i]=sqrt(Sqr(f.Re[i])*Sqr(f.Im[i]));
}
//FFT用データ配列の初期化
free(f.Re),f.Re=NULL;
free(f.Im),f.Im=NULL;
//出力データの書き込み
ofp=fopen("c:\\make\\data\\output.csv","w");
for(i=0;i<DN/2;i++){
fprintf(ofp,"%lf,%lf\n",i*dF, amp[i]);
}
fclose(ofp);
return 0;
}
//fft.h
extern void FFT(struct DATA);
#define Sqr(x) ((x)*(x))
data.h
int g, h, i, j, k, l, m, n, p, q, DN;
double a, b, tmp;
double dumx[1024];
double dumy[1024];
double *time;
struct DATA {
double *Re;
double *Im;
};
struct ROT {
double *c;
double *s;
};
//FFT解析用ソースコード
//analysis.c
#include "data.h"
#include <math.h>
#include <stdlib.h>
#define PI 3.1415926
void Table( struct ROT);
void swap( double*, double*);
void FFT(struct DATA in)
{
struct ROT r;
r.c=NULL;
r.s=NULL;
in.Re=(double*)calloc(DN,sizeof(double));
in.Im=(double*)calloc(DN,sizeof(double));
n=DN;
m = (int)(log10(n)/log10(2));
//回転因子テーブルの作成
Table(r);
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=in.Re[i]-in.Re[j], b=in.Im[i] - in.Im[j];
in.Re[i] = in.Re[i] + in.Re[j], in.Im[i] = in.Im[i] + in.Im[j];
if(p==0){
in.Re[j]=a,in.Im[j]=b;
}else{
in.Re[j] = a * (r.c[p]) + b * (r.s[p]), in.Im[j] = b * (r.c[p]) - a * (r.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(&(in.Re[i]),&(in.Re[j]));
swap(&(in.Im[i]),&(in.Im[j]));
}
k/=2;
while(j>=k){
j=j-k,k /=2;
//無限ループ回避(ここを消すと無限ループに陥ります。)ここから
if((j==0)&&(k==0)){
break;
}
//ここまで
}
j += k;
}
}
void Table( struct ROT w)
{
w.c = (double*)calloc(DN/2,sizeof(double));
w.s = (double*)calloc(DN/2,sizeof(double));
a=0.00, b = 2.0 * PI /DN;
for(i=0;i<DN/2;i++)
{
(w.c[i])=cos(a+i*b), (w.s[i])=sin(a+i*b);
}
}
void swap(double *var1, double *var2)
{
tmp = *var1, *var1=*var2, *var2=tmp;
}
補足
すみません。もう少し具体的に書きます。 FFTは N-1 Ak = Σ A(n)・W^kn, W^kn = e^-2πn/N にて示されますが、 n=0 サイトでよく出ているFFT算出Cプログラムを見ると、 Re0,Im0,Re1,Im1という実数＋虚数入力に対して下記のように演算されています。 +----------------------------------------------+ | Wre = cos, Wre = sinとしてTableを持ち | | | | XRe(0) = Re(0) + Wre x Re(1) - Wim x Im(1) | | XIm(0) = Im(0) + Wim x Re(1) + Wre x Im(1) | | XRe(1) = Re(0) - Wre x Re(1) + Wim x Im(1) | | XIm(1) = Im(0) - Wim x Re(1) - Wre x Im(1) | +----------------------------------------------+ 単純にDITのバタフライ図を書いてみると(見づらくてすみません) Re(0) + jIm(0) -+-+- XRe(0) + jXIm(0) X Re(1) + jIm(1) -+-+- XRe(1) + jXIm(1) W - 上記のバタフライを単純計算すると、 XRe(0) + jXIm(0) = Re(0) + jIm(0) + W(Re(1) + jIm(1)) = Re(0) + jIm(0) + (cos - jsin)(Re(1) + jIm(1)) = Re(0) + jIm(0) + cos・Re(1) - jsin・Re(1) + jcos・Im(1) - j^2sin・Im(1) = Re(0) + cos・Re(1) + sin・Im(1) + j(Im(0) - sin・Re(1) + cos・Im) ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 実数項:XRe(0) 虚数項：Xim(0) XRe(1) + jXIm(1) = Re(0) + jIm(0) - W(Re(1) + jIm(1)) = Re(0) + jIm(0) - (cos - jsin)(Re(1) + jIm(1)) = Re(0) + jIm(0) - cos・Re(1) + jsin・Re(1) - jcos・Im(1) + j^2sin・Im(1) = Re(0) - cos・Re(1) - sin・Im(1) + j(Im(0) + sin・Re(1) - cos・Im) ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 実数項:XRe(1) 虚数項：Xim(1) ここで、単純計算で求めた、各項とCプログラムの式の各項を比較してみると 微妙に符号が合いません。どこで勘違いしているかご指摘頂ければ助かります。