【DSP C++】2D_FFT
来源：互联网 发布：前端后端数据库的关系 编辑：程序博客网 时间：2024/05/22 01:58
#include <stdio.h>#include <stdlib.h>#include <math.h>#define intsize sizeof(int)#define complexsize sizeof(complex)#define PI 3.1415926int *a,*b;int nLen,init_nLen,mLen,init_mLen,N,M;//init_nLen,init_mLen：原始图像的大小、N,M：dft变换应该是2的几次方、nLen，mLen：dft的点数FILE *dataFile;typedef struct{ float real; float image;}complex;complex *A,*A_In,*W;//A_In:dft的原矩阵，不够的补0complex Add(complex, complex);complex Sub(complex, complex);complex Mul(complex, complex);int calculate_M(int);//计算dft的点数应为2的N次方、void reverse(int,int);void readData(); //读入原矩阵数据，并补零转化为可以计算dft的矩阵、void fft(int,int);void Ifft();void printResult_fft();void printResult_Ifft();int main(){ int i,j; readData(); A = (complex *)malloc(complexsize*nLen); reverse(nLen,N); for(i=0; i<mLen; i++) { for(j=0; j<nLen; j++) { A[j].real = A_In[i*nLen+b[j]].real; A[j].image = A_In[i*nLen+b[j]].image; } fft(nLen,N); for(j=0; j<nLen; j++) { A_In[i*nLen+j].real = A[j].real; A_In[i*nLen+j].image = A[j].image; } } free(a); free(b); free(A); A = (complex *)malloc(complexsize*mLen); reverse(mLen,M); for(i=0; i<nLen; i++) { for(j=0; j<mLen; j++) { A[j].real = A_In[b[j]*nLen+i].real; A[j].image = A_In[b[j]*nLen+i].image; } fft(mLen,M); for(j=0; j<mLen; j++) { A_In[j*nLen+i].real = A[j].real; A_In[j*nLen+i].image = A[j].image; } } free(A); printResult_fft(); Ifft(); printResult_Ifft(); return 0;}void readData(){ int i,j; dataFile = fopen("dataIn.txt","r"); fscanf(dataFile,"%d %d",&init_mLen,&init_nLen); M = calculate_M(init_mLen); N = calculate_M(init_nLen); nLen = (int)pow(2,N); mLen = (int)pow(2,M); A_In = (complex *)malloc(complexsize*nLen*mLen); for(i=0; i<init_mLen; i++) { for(j=0; j<init_nLen; j++) { fscanf(dataFile,"%f",&A_In[i*nLen+j].real); A_In[i*nLen+j].image = 0.0; } } fclose(dataFile); for(i=0; i<mLen; i++) { for(j=init_nLen; j<nLen; j++) { A_In[i*nLen+j].real = 0.0; A_In[i*nLen+j].image = 0.0; } } for(i=init_mLen; i<mLen; i++) { for(j=0; j<init_nLen; j++) { A_In[i*nLen+j].real = 0.0; A_In[i*nLen+j].image = 0.0; } } printf("Reading initial datas:\n"); for(i=0; i<init_mLen; i++) { for(j=0; j<init_nLen; j++) { if(A_In[i*nLen+j].image < 0) { printf("%f%fi\t",A_In[i*nLen+j].real,A_In[i*nLen+j].image); } else { printf("%f+%fi\t",A_In[i*nLen+j].real,A_In[i*nLen+j].image); } } printf("\n"); } printf("\n"); printf("Reading formal datas:\n"); for(i=0; i<mLen; i++) { for(j=0; j<nLen; j++) { if(A_In[i*nLen+j].image < 0) { printf("%f%fi\t",A_In[i*nLen+j].real,A_In[i*nLen+j].image); } else { printf("%f+%fi\t",A_In[i*nLen+j].real,A_In[i*nLen+j].image); } } printf("\n"); }}void fft(int fft_nLen, int fft_M){ int i; int lev,dist,p,t; complex B; W = (complex *)malloc(complexsize*fft_nLen/2); for(lev=1; lev<=fft_M; lev++) { dist = (int)pow(2,lev-1); for(t=0; t<dist; t++) { p = t*(int)pow(2,fft_M-lev); W[p].real = (float)cos(2*PI*p/fft_nLen); W[p].image = (float)(-1*sin(2*PI*p/fft_nLen)); for(i=t; i<fft_nLen; i=i+(int)pow(2,lev)) { B = Add(A[i],Mul(A[i+dist],W[p])); A[i+dist] = Sub(A[i],Mul(A[i+dist],W[p])); A[i].real = B.real; A[i].image = B.image; } } } free(W);}void printResult_fft(){ int i,j; printf("Output FFT results:\n"); for(i=0; i<mLen; i++) { for(j=0; j<nLen; j++) { if(A_In[i*nLen+j].image < 0) { printf("%f%fi\t",A_In[i*nLen+j].real,A_In[i*nLen+j].image); } else { printf("%f+%fi\t",A_In[i*nLen+j].real,A_In[i*nLen+j].image); } } printf("\n"); }}void printResult_Ifft(){ int i,j; printf("Output IFFT results:\n"); for(i=0; i<mLen; i++) { for(j=0; j<nLen; j++) { if(A_In[i*nLen+j].image < 0) { printf("%f%fi\t",A_In[i*nLen+j].real,A_In[i*nLen+j].image); } else { printf("%f+%fi\t",A_In[i*nLen+j].real,A_In[i*nLen+j].image); } } printf("\n"); } free(A_In);}int calculate_M(int len){ int i; int k; i = 0; k = 1; while(k < len) { k = k*2; i++; } return i;}void reverse(int len, int M){ int i,j; a = (int *)malloc(intsize*M); b = (int *)malloc(intsize*len); for(i=0; i<M; i++) { a[i] = 0; } b[0] = 0; for(i=1; i<len; i++) { j = 0; while(a[j] != 0) { a[j] = 0; j++; } a[j] = 1; b[i] = 0; for(j=0; j<M; j++) { b[i] = b[i]+a[j]*(int)pow(2,M-1-j); } }}complex Add(complex c1, complex c2){ complex c; c.real = c1.real+c2.real; c.image = c1.image+c2.image; return c;}complex Sub(complex c1, complex c2){ complex c; c.real = c1.real-c2.real; c.image = c1.image-c2.image; return c;}complex Mul(complex c1, complex c2){ complex c; c.real = c1.real*c2.real-c1.image*c2.image; c.image = c1.real*c2.image+c2.real*c1.image; return c;}void Ifft(){ int i,j; for(i=0; i<mLen; i++) { for(j=0; j<nLen; j++) { A_In[i*nLen+j].image = -A_In[i*nLen+j].image; } } A = (complex *)malloc(complexsize*nLen); reverse(nLen,N); for(i=0; i<mLen; i++) { for(j=0; j<nLen; j++) { A[j].real = A_In[i*nLen+b[j]].real; A[j].image = A_In[i*nLen+b[j]].image; } fft(nLen,N); for(j=0; j<nLen; j++) { A_In[i*nLen+j].real = A[j].real/nLen; A_In[i*nLen+j].image = A[j].image/nLen; } } free(A); free(a); free(b); A = (complex *)malloc(complexsize*mLen); reverse(mLen,M); for(i=0; i<nLen; i++) { for(j=0; j<mLen; j++) { A[j].real = A_In[b[j]*nLen+i].real; A[j].image = A_In[b[j]*nLen+i].image; } fft(mLen,M); for(j=0; j<mLen; j++) { A_In[j*nLen+i].real = A[j].real/mLen; A_In[j*nLen+i].image = A[j].image/mLen; } } free(A); free(a); free(b);}测试数据及结果如下：数据输入文件data_In.txt的内容如下：3 31 2 5 3 2 5 5 6 4测试结果如下：Reading initial datas:1.000000+0.000000i 2.000000+0.000000i 5.000000+0.000000i3.000000+0.000000i 2.000000+0.000000i 5.000000+0.000000i5.000000+0.000000i 6.000000+0.000000i 4.000000+0.000000iReading formal datas:1.000000+0.000000i 2.000000+0.000000i 5.000000+0.000000i 0.000000+0.000000i3.000000+0.000000i 2.000000+0.000000i 5.000000+0.000000i 0.000000+0.000000i5.000000+0.000000i 6.000000+0.000000i 4.000000+0.000000i 0.000000+0.000000i0.000000+0.000000i 0.000000+0.000000i 0.000000+0.000000i 0.000000+0.000000iOutput FFT results:33.000000+0.000000i -5.000000-10.000000i 13.000000+0.000000i -5.000000+10.000000i-7.000000-10.000000i -7.000000+6.000000i 1.000000-6.000000i -3.000000-2.000000i13.000000+0.000000i -1.000000-6.000000i 1.000000+0.000000i -1.000000+6.000000i-7.000000+10.000000i -3.000000+2.000000i 1.000000+6.000000i -7.000000-6.000000iOutput IFFT results:1.000000+0.000000i 2.000000+0.000000i 5.000000+0.000000i 0.000000-0.000000i3.000000+0.000000i 2.000000+0.000000i 5.000000+0.000000i 0.000000-0.000000i5.000000+0.000000i 6.000000+0.000000i 4.000000+0.000000i -0.000000-0.000000i0.000000-0.000000i 0.000000+0.000000i 0.000000-0.000000i -0.000000+0.000000i 0 0
- 【DSP C++】2D_FFT
- DSP的C语言学习
- DSP C语言基础要点
- DSP C语言基础要点
- DSP C语言基础要点
- DSP系统关于C编程
- DSP开发-C语言环境
- dsp C优化-（四）
- DSP的C/C++开发
- DSP的C/C++语言程序设计
- DSP帮手(2)
- DSP总结2
- dsp调试-----2
- dsp
- dsp
- DSP
- DSP
- dsp
- 计算平方根（牛顿迭代法）
- js基础
- NSUserDefaults 与 数据的归档和反归档
- 详解kernel中watchdog 驱动程序
- Linux性能监测：网络篇
- 【DSP C++】2D_FFT
- js基础（2）
- 再造人类
- ObjectInputStream循环读取对象的方法
- JS基础（3）
- 在iOS应用中跳转到AppStore
- 计算机重装系统后启动出现grub光标解决办法
- Linux性能监测：磁盘IO篇
- 2014Unity亚洲开发者大会会议简录之技术篇
