FFTをやりたいのですが、


	BBS

FFTをやりたいのですが、

Nih - 2002年04月03日 22時15分39秒

WINDOWS95下でVC++6で使っています。
0.1秒間に4096サンプルデータある波形データ
（XR[0]からXR[4095]まで）をFFTしたいのですが、、、

①　結果が、どうもしっくり行きません。
　　バタフライ演算がうまく行っていないせいか、
　　新たに算出された
　　XR[1]+iXI[1]の複素スペクトルの絶対値と
　　XR[4]+iXI[4]の複素スペクトルの絶対値が
　　イコールになったりして。。。）
　　コードは下記です。
　　4096データがあるので、2＾12=4096になるので、
　　バタフライをやる時の処理列は12列にしているのですが、
　　どこがおかしいのでしょうか？

②　また、下記コードで、
　　フーリエ変換した後にはXR[FI]でフーリエスペクトルの
　　実部、XI[FI]に同じく虚部を入れているのですが、
　　周波数との対応はどのようになるのでしょうか？
　　（要はXR[0]+iXI[0]がどの周波数の複素スペクトルで、
　　　XR[1]+iXI[1]ががどの周波数の複素スペクトルで、
　　　.............
　　　.............

　　　XR[4095]+iXI[4095]ががどの周波数の複素スペクトル
　　　に対応するのか？）

②　また、逆フーリエ変換をする場合においては、
　　下記コードで、何所がどう変わって来て、
　　波形のサンプリング値はどの変数に戻って来るのでしょうか？

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
/* XR[0]～XR[4095]に前もって波形のサンプリングデータが　*/
/* 入っていて、これをフーリエ変換したい。　*/

double XR[4096]; /* For FFT */
double XI[4096]; /* For FFT */

int FN; /* For FFT */
int FL; /* For FFT */
int FR; /* For FFT */
int FIR; /* For FFT */
int FJ2; /* For FFT */
double FDEG; /* For FFT */
int FK; /* For FFT */
int FL1; /* For FFT */
double FP; /* For FFT */
int FB1; /* For FFT */
int FBIT; /* For FFT */
int FB2; /* For FFT */
double FARG; /* For FFT */
double FC; /* For FFT */
double FS; /* For FFT */
int FK1; /* For FFT */
double FA; /* For FFT */
double FB; /* For FFT */
int FKK; /* For FFT */
int FI; /* For FFT */
int FII; /* For FFT */
double FELR; /* For FFT */
double FELI; /* For FFT */
double PWR[4096]; /* For FFT */　

FR = 12; /* FFT Routine */
　　 FIR = 1; /* FFT Routine */
FN = 4096; 　　　　　　　/* 下まで */
FJ2 = FN;
FDEG = 2*3.14159/FN;
FK = 0;
FL1 = FR - 1;
for (FL=1;FL<=FR;FL++)
{
FJ2 = FJ2/2;
for (FI=1;FI<=FJ2;FI++)
{
FP = FK/pow(2,FL1);
FB1 = FP;
FBIT = 0;
for (FII=1;FII<=FR;FII++)
{
FB2 = FB1/2;
FBIT = FBIT*2+FB1-2*FB2;
FB1 = FB2;
}
FARG = FDEG*FBIT;
FC = cos(FARG);
FS = FIR*sin(FARG);
FK1 = FK + FJ2;
FA = XR[FK1]*FC+XI[FK1]*FS;
FB = XI[FK1]*FC-XR[FK1]*FS;
XR[FK1]=XR[FK]-FA;
XI[FK1]=XI[FK]-FB;
XR[FK]=XR[FK]+FA;
XI[FK]=XI[FK]+FB;
FK = FK + 1;
}
FK = FK + FJ2;
if (FK >= FN)
{
FK = 0;
FL1 = FL1 - 1;
}
}

for (FKK = 0;FKK<=FN-1;FKK++)
{
FP = FKK;
FB1 = FP;
FBIT = 0;
for (FII=1;FII<=FR;FII++)
{
FB2 = FB1/2;
FBIT = FBIT*2+FB1-2*FB2;
FB1 = FB2;
}
if (FBIT > FKK)
{
FELR = XR[FKK];
XR[FKK] = XR[FBIT];
XR[FBIT] = FELR;
FELI = XI[FKK];
XI[FKK] = XI[FBIT];
XI[FBIT] = FELI;
}
}

for (FI=0;FI<FN-1;FI++)
{
PWR[FI]=sqrt(pow(XR[FI],2)+pow(XR[FI],2));
printf("FN FI XR XI PWR = %d %d %d %d \n",
FN,FI,XR[FI],XI[FI],PWR[FI]);
}

この投稿へのコメント

追加質問です。

| インデックス | ホームページ |