高速フーリエ変換プログラムのC言語ソース高速フーリエ変換プログラムのC言語ソースを以下に示します。
なお、データ入力とデータ出力の部分は省略しています。
#include
#include
#include
int L,n;
double X[4100],Y[4100],f;
void main()
{
int i,j,ns,L1,arg,arg2,i0,i1,m[4100],n1;
double s,c,sc,X1,Y1,t,PI;
IN_FM();//データ入力
PI=3.14159265358979;//円周率を設定
n=(int)pow(2,L);//データ数
n1=n/2;//データ数の1/2
ns=n/2;//2進数の最上位の値
sc=2.0*PI/(double)n;
for(i=0;i < n;i++){m[i]=0;} //m[i]の初期化
while(ns >= 1) //2進数の最上位の値が1になるまで繰り返し
{
for(L1=0;L1 < n;L1=L1+2*ns) //L1の値を0からnまで、2×最上位の値のstepで繰り返し
{
arg=m[L1]/2; //回転子=前の回転子の値/2
arg2=arg+n1; //回転子2＝回転子+データ数の1/2
c=cos(sc*(double)arg);
s=sin(f*sc*(double)arg);
for(i0=L1;i0 < L1+ns;i0++) //i0の値をL1の値からL1+2進数の最上位の値まで繰り返し
{
i1=i0+ns; //i1＝i0+2進数の最上位の値
X1=X[i1]*c-Y[i1]*s;
Y1=Y[i1]*c+X[i1]*s;
X[i1]=X[i0]-X1;
Y[i1]=Y[i0]-Y1;
X[i0]=X[i0]+X1;
Y[i0]=Y[i0]+Y1;
m[i0]=arg; //m[i0]に回転子を保存
m[i1]=arg2; //m[i1]に回転子2を保存
}
}
ns=ns/2;//2進数の最上位の値を繰り下げ
}
if(f < 0)//正規化処理
{
for(i=0;i < n; i++)
{
X[i]=X[i]/(double)n;
Y[i]=Y[i]/(double)n;
}
}
for(i=0;i < n; i++)//ビット反転処理
{
if(i < m[i])
{
j=m[i];
t=X[i];X[i]=X[j];X[j]=t;
t=Y[i];Y[i]=Y[j];Y[j]=t;
}
}
OUT_FM();//データ出力
}
高速フーリエ変換の基礎理論と比較してプログラムのソースが短いと感じるかもしれません？
しかし、短くても難解なプログラムのソースなのです。
４章：高速フーリエ変換プログラムのVBAソースに行く。