自定義結構體complex（復數）

typedef struct{double re; //實部 double im; //虛部 
}complex;

用于為離散序列倒碼排序的函數

int  reverse(int n,int M){int i,m,dev,temp=0;for(i=0,m=1,dev=1;i<M;i++,m*=2,dev*=2){temp=(temp<<1)+(n&m)/dev;}return temp;
}void reSequence(complex*src, complex*dst,int N){int i,M=(int)log2(N);for(i=0;i<N;i++){dst[i].re=src[reverse(i,M)].re;dst[i].im=src[reverse(i,M)].im;}
}

dft、idft、fft、ifft實現代碼

void dft(complex*src,complex*dst,int N){int i,j; double temp;for(i=0;i<N;i++){dst[i].re=0;dst[i].im=0;for(j=0;j<N;j++){temp=-2*pi*j*i/N;dst[i].re+=src[j].re*cos(temp)-src[j].im*sin(temp);dst[i].im+=src[j].re*sin(temp)+src[j].im*cos(temp);}}
}void idft(complex*src,complex*dst,int N){int i,j;double temp;for(i=0;i<N;i++){dst[i].re=0;dst[i].im=0;for(j=0;j<N;j++){temp=2*pi*i*j/N;dst[i].re+=(src[j].re*cos(temp)-src[j].im*sin(temp))/N;dst[i].im+=(src[j].re*sin(temp)+src[j].im*cos(temp))/N;}}
}void fft(complex*src,complex*dst,int N){int i,j;double temp;int g,n,gap;complex x1,x2;reSequence(src,dst,N);for(g=N/2;g>0;g/=2){	//g為蝴蝶單元群的數量 n=N/g;                 //n為每一蝴蝶單元群中蝴蝶單元的數量 gap=n/2;			//gap為蝴蝶單元上兩個節點的間距 for(i=0;i<N;i+=n){for(j=i;j<i+n/2;j++){temp=-2*pi*g*(j-i)/double(N);x1.re=dst[j].re+dst[j+gap].re*cos(temp)-dst[j+gap].im*sin(temp);x1.im=dst[j].im+dst[j+gap].re*sin(temp)+dst[j+gap].im*cos(temp);x2.re=dst[j].re-dst[j+gap].re*cos(temp)+dst[j+gap].im*sin(temp);x2.im=dst[j].im-dst[j+gap].re*sin(temp)-dst[j+gap].im*cos(temp);dst[j].re=x1.re;dst[j].im=x1.im;dst[j+gap].re=x2.re;dst[j+gap].im=x2.im;}}}
}void ifft(complex*src,complex*dst,int N){int i,j;double temp;int g,n,gap;complex x1,x2;reSequence(src,dst,N);for(g=N/2;g>0;g/=2){	//g為蝴蝶單元群的數量 n=N/g;                 //n為每一蝴蝶單元群中蝴蝶單元的數量 gap=n/2;			//gap為蝴蝶單元上兩個節點的間距 for(i=0;i<N;i+=n){for(j=i;j<i+n/2;j++){temp=2*pi*g*(j-i)/double(N);x1.re=(dst[j].re+dst[j+gap].re*cos(temp)-dst[j+gap].im*sin(temp));x1.im=(dst[j].im+dst[j+gap].re*sin(temp)+dst[j+gap].im*cos(temp));x2.re=(dst[j].re-dst[j+gap].re*cos(temp)+dst[j+gap].im*sin(temp));x2.im=(dst[j].im-dst[j+gap].re*sin(temp)-dst[j+gap].im*cos(temp));dst[j].re=x1.re;dst[j].im=x1.im;dst[j+gap].re=x2.re;dst[j+gap].im=x2.im;}}}for(i=0;i<N;i++){dst[i].re/=N;dst[i].im/=N;}
}

總結

以上是生活随笔為你收集整理的C语言实现傅里叶变换函数dft,idft,fft,ifft的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。