取样光纤光栅传输函数源程序.docxVIP

源程序：#include complex#include iostream.h#include math.h#define pi 3define N1 8#define N2 8typedef std::complexdouble complex;/*----二维基2FFT算法----*/void fft2(complex f[N1][N2]){int i,j,k,m,n,l,r,M;int la,lb,lc;complex t;/*----逐行进行一维FFT变换----*/for(i=0;iN1;i++){/*----时域抽取的FFT算法,对原始图像矩阵第i行数据进行FFT变换----*//*----计算分解的级数M=log2(N2)----*/for(m=N2,M=1;(m=m/2)!=1;M++); /*----按照倒位序重新排列原信号----*/for(j=1,k=N2/2;j=N2-2;j++){if(jk){t=f[i][k];f[i][k]=f[i][j];f[i][j]=t;}l=N2/2;while(l=k){k=k-l;l=l/2;}k=k+l;}/*----FFT算法----*/for(m=1;m=M;m++){la=pow(2,m); //la=2^m代表第m级每个分组所含节点数lb=la/2; //lb代表第m级每个分组所含蝶形单元数 //同时它也表示每个蝶形单元上下节点之间的距离/*----蝶形运算----*/for(l=1;l=lb;l++){r=(l-1)*pow(2,M-m);for(n=l-1;nN2-1;n=n+la) //遍历每个分组，分组总数为N2/la{lc=n+lb; //n,lc分别代表一个蝶形单元的上、下节点编号t=f[i][lc]*complex(cos(2*pi*r/N2),-sin(2*pi*r/N2));f[i][lc]=f[i][n]-t;f[i][n]=f[i][n]+t;}}}}/*----逐列进行一维FFT变换----*/for(j=0;jN2;j++){/*----时域抽取的FFT算法,对原始图像矩阵第j列数据进行FFT变换----*//*----计算分解的级数M=log2(N1)----*/for(m=N1,M=1;(m=m/2)!=1;M++); /*----按照倒位序重新排列原信号----*/for(i=1,k=N1/2;i=N1-2;i++){if(ik){t=f[k][j];f[k][j]=f[i][j];f[i][j]=t;}l=N1/2;while(l=k){k=k-l;l=l/2;}k=k+l;}/*----FFT算法----*/for(m=1;m=M;m++){la=pow(2,m); //la=2^m代表第m级每个分组所含节点数lb=la/2; //lb代表第m级每个分组所含蝶形单元数 //同时它也表示每个蝶形单元上下节点之间的距离/*----蝶形运算----*/for(l=1;l=lb;l++){r=(l-1)*pow(2,M-m);for(n=l-1;nN1-1;n=n+la) //遍历每个分组，分组总数为N1/la{lc=n+lb; //n,lc分别代表一个蝶形单元的上、下节点编号t=f[lc][j]*complex(cos(2*pi*r/N1),-sin(2*pi*r/N1));f[lc][j]=f[n][j]-t;f[n][j]=f[n][j]+t;}}}}}/*----显示信号数据----*/void display(complex f[N1][N2]){int i,j;for(i=0;iN1;i++){for(j=0;jN2;j++){if((j)%3==0) coutendl;cout.width(9);cout.setf(ios::fixed);cout.precision(4);coutf[i][j].real();cout.flags(ios::showpos);cout.width(9);cout.setf(ios::fixed);cout.precision(4);coutf[i][j].imag()i\t;cout.unsetf(ios::showpos);}}}/*----主函数----*/void main(){int i,j;/*----原始图像矩阵数据x[N1][N2]----*/double x[N1][N2]={{0,0,0,0,0,0,0,0},{0,0,0,1,1,0,0,0},{0,0,2,4,4,2,0,0},{0,1,4,8,8,4,1,0},{0,1,4,8,8,4,1,0},{0,0,2,4,4,