C语言数据结构算法之实现快速傅立叶变换
C语言数据结构算法之实现快速傅立叶变换
本实例将实现二维快速傅立叶变换,同时也将借此实例学习用c语言实现矩阵的基本操作、复数的基本掾作,复习所学过的动态内存分配、文件操作、结构指针的函数调用等内容。
很久以来,傅立叶变换一直是许多领域,如线性系统、光学、概率论、量子物理、天线、数字图像处理和信号分析等的一个基本分析工具,但是即便使用计算速度惊人的计算机计算离散傅立叶变换所花费的时间也常常是难以接受的,因此导致了快速傅立叶变换(FFT)的产生。
本实例将对一个二维数组进行正、反快速傅立叶变换。正傅立叶变换时dfft()函数先调用fft()按行对数组进行变换,再对结果调用fft()按列进行变换,此时完成了快速傅立叶变换,再调用rdfft()函数进行傅立叶逆变换。如果程序设计正确的话,变换的结果应与原数组相同。
实例代码:
#include#include #include #definePI3.14159265358979323846 structCOMPLEX { floatre; floatim; }cplx,*Hfield,*S,*R,*w; intn,m; intln,lm; voidinitiate(); voiddfft(); voidrdfft(); voidshowresult(); voidfft(intl,intk); intreverse(intt,intk); voidW(intl); intloop(intl); voidconjugate(); voidadd(structCOMPLEX*x,structCOMPLEX*y,structCOMPLEX*z); voidsub(structCOMPLEX*x,structCOMPLEX*y,structCOMPLEX*z); voidmul(structCOMPLEX*x,structCOMPLEX*y,structCOMPLEX*z); structCOMPLEX*Hread(inti,intj); voidHwrite(inti,intj,structCOMPLEXx); voidmain() { initiate(); printf("\n原始数据:\n"); showresult(); getchar(); dfft(); printf("\n快速复利叶变换后的结果:\n"); showresult(); getchar(); rdfft(); printf("\n快速复利叶逆变换后的结果:\n"); showresult(); getchar(); free(Hfield); } voidinitiate() {//程序初始化操作,包括分配内存、读入要处理的数据、进行显示等 FILE*df; df=fopen("data.txt","r"); fscanf(df,"%5d",&n); fscanf(df,"%5d",&m); if((ln=loop(n))==-1) { printf("列数不是2的整数次幂"); exit(1); } if((lm=loop(m))==-1) { printf("行数不是2的整数次幂"); exit(1); } Hfield=(structCOMPLEX*)malloc(n*m*sizeof(cplx)); if(fread(Hfield,sizeof(cplx),m*n,df)!=(unsigned)(m*n)) { if(feof(df))printf("Prematureendoffile"); elseprintf("Filereaderror"); } fclose(df); } voiddfft() {//进行二维快速复利叶变换 inti,j; intl,k; l=n; k=ln; w=(structCOMPLEX*)calloc(l,sizeof(cplx)); R=(structCOMPLEX*)calloc(l,sizeof(cplx)); S=(structCOMPLEX*)calloc(l,sizeof(cplx)); W(l); for(i=0;i re; S[j].im=Hread(i,j)->im; } fft(l,k); for(j=0;j re; S[j].im=Hread(j,i)->im; } fft(l,k); for(j=0;j re,Hread(i,j)->im); } } } voidfft(intl,intk) { inti,j,s,nv,t; floatc; structCOMPLEXmp,r; nv=l; c=(float)l; c=pow(c,0.5); for(i=0;i >(k-i-1); s=reverse(s,k); r.re=S[t+j].re; r.im=S[t+j].im; mul(&w[s],&S[t+j+nv/2],&mp);/////////讲解传递结构指针和结构本身的区别 add(&r,&mp,&S[t+j]); sub(&r,&mp,&S[t+j+nv/2]); } } nv=nv>>1; } for(i=0;i >1; y=(y<<1)+x; } returny; } voidW(intl) { inti; floatc,a; c=(float)l; c=2*PI/c; for(i=0;i >i; if(m==0) break; } if(l==(1<<(i-1))) returni-1; } return-1; } voidconjugate() {//求复数矩阵的共轭矩阵 inti,j; for(i=0;i im*=-1; } } } structCOMPLEX*Hread(inti,intj) {//按读矩阵方式返回Hfield中指定位置的指针 return(Hfield+i*n+j); } voidHwrite(inti,intj,structCOMPLEXx) {//按写矩阵方式将复数结构x写到指定的Hfield位置上 (Hfield+i*n+j)->re=x.re; (Hfield+i*n+j)->im=x.im; } voidadd(structCOMPLEX*x,structCOMPLEX*y,structCOMPLEX*z) {//定义复数加法 z->re=x->re+y->re; z->im=x->im+y->im; } voidsub(structCOMPLEX*x,structCOMPLEX*y,structCOMPLEX*z) {//定义复数减法 z->re=x->re-y->re; z->im=x->im-y->im; } voidmul(structCOMPLEX*x,structCOMPLEX*y,structCOMPLEX*z) {//定义复数乘法 z->re=(x->re)*(y->re)-(x->im)*(y->im); z->im=(x->im)*(y->re)+(x->re)*(y->im); }
感谢阅读,希望能帮助到大家,谢谢大家对本站的支持!