基2基3定点混合基fft

#include "StdAfx.h" void fixedpfft(long ixr[POINT],long ixi[POINT],double oxr[POINT],double oxi[POINT],int R,int P,int Q) { int n,n1,n0,i,j,k,l,t,par1r,par1i,par2r,par2i,sq2,sq3; // long *mr=new long[POINT],*mi=new long[POINT],(*ar)[POINT]=new long[POINT][POINT],(*ai)[POINT]=new long[POINT][POINT]; long mr[POINT],mi[POINT],(*ar)[POINT]=new long[POINT][POINT],(*ai)[POINT]=new long[POINT][POINT]; // long *wr=new long[POINT],*wi=new long[POINT],*w2r=new long[POINT],*w2i=new long[POINT],*w3r=new long[POINT],*w3i=new long[POINT]; long wr[POINT],wi[POINT],w2r[POINT],w2i[POINT],w3r[POINT],w3i[POINT]; long tmp,upr,upi,midr,midi,downr,downi,tmp1r,tmp1i,tmp2r,tmp2i; sq2 = 1/sqrt(2.00)*MAG+0.5; sq3 = 1/sqrt(3.00)*MAG+0.5; par1r = par2r = -0.5*MAG-0.5; par1i = -sqrt(3.00)/2.00*MAG-0.5; par2i = sqrt(3.00)/2.00*MAG+0.5; for(n=0;n<POINT;n++) { mr[n] = ixr[n]; mi[n] = ixi[n]; } //计算旋转因子； // cout<<"旋转因子"<<endl; for(i=0;i<POINT;i++) { wr[i] = cos(2*PI/POINT*i)*MAG; wi[i] = -1*sin(2*PI/POINT*i)*MAG; // cout<<"w["<<i<<"]="<<wr[i]<<"+j*"<<wi[i]<<endl; } // cout<<"------------------------------------------------"<<endl; //计算基二的旋转因子; // cout<<"基二的旋转因子:"<<endl; for(i=0;i<R;i++) { w2r[i] = cos(2*PI/R*i)*MAG; w2i[i] = -1*sin(2*PI/R*i)*MAG; // cout<<"w2["<<i<<"]="<<w2r[i]<<"+j*"<<w2i[i]<<endl; } // cout<<"------------------------------------------------"<<endl; //计算基三的旋转因子; // cout<<"基三的旋转因子:"<<endl; for(i=0;i<P;i++) { w3r[i] = cos(2*PI/P*i)*MAG; w3i[i] = -1*sin(2*PI/P*i)*MAG; // cout<<"w3["<<i<<"]="<<w3r[i]<<"+j*"<<w3i[i]<<endl; } // cout<<"------------------------------------------------"<<endl; //计算基五的旋转因子 //分解为二维数组; // cout<<"分解为二维数组:"<<endl; for(n1=0;n1<R;n1++) { for(n0=0;n0<P;n0++) { ar[n1][n0] = mr[P*n1+n0]; ai[n1][n0] = mi[P*n1+n0]; // cout<<"ax["<<n1<<"]["<<n0<<"]="<<ar[n1][n0]<<"+j*"<<ai[n1][n0]<<endl; } } // cout<<"------------------------------------------------"<<endl; //对n1部分进行二进制倒位序排序; for(n0=0;n0<P;n0++) { for(i=0;i<R;i++) { k=i; j=0; t=log(R*1.00)/log(2.00); while((t--)>0) { j=j<<1; j=j|(k&1); k=k>>1; } if(j>i) { tmp = ar[i][n0]; ar[i][n0] = ar[j][n0]; ar[j][n0] = tmp; tmp = ai[i][n0]; ai[i][n0] = ai[j][n0]; ai[j][n0] = tmp; } } } /* cout<<"二进制倒位序排序后的序列："<<endl; for(n1=0;n1<R;n1++) { for(n0=0;n0<P;n0++) { cout<<"ax["<<n1<<"]["<<n0<<"]="<<ar[n1][n0]<<"+j*"<<ai[n1][n0]<<endl; } } cout<<"------------------------------------------------"<<endl; */ //基二fft; for(n0=0;n0<P;n0++) { for(i=0;i<log(R*1.00)/log(2.00);i++) { l=pow(2,i); for(j=0;j<R;j+=2*l) { for(k=0;k<l;k++) { tmp1r = ar[j+k+l][n0]*w2r[R*k/2/l]/MAG-ai[j+k+l][n0]*w2i[R*k/2/l]/MAG; tmp1i = ar[j+k+l][n0]*w2i[R*k/2/l]/MAG+ai[j+k+l][n0]*w2r[R*k/2/l]/MAG; upr = ar[j+k][n0] + tmp1r; upi = ai[j+k][n0] + tmp1i; downr = ar[j+k][n0] - tmp1r; downi = ai[j+k][n0] - tmp1i; ar[j+k][n0] = upr*sq2/MAG; ai[j+k][n0] = upi*sq2/MAG; ar[j+k+l][n0] = downr*sq2/MAG; ai[j+k+l][n0] = downi*sq2/MAG; } } } } /* cout<<"基二fft后的序列："<<endl; for(n1=0;n1<R;n1++) { for(n0=0;n0<P;n0++) { cout<<"ax["<<n1<<"]["<<n0<<"]="<<ar[n1][n0]<<"+j*"<<ai[n1][n0]<<endl; } } cout<<"------------------------------------------------"<<endl; */ //乘以旋转因子； for(n1=0;n1<R;n1++) for(n0=0;n0<P;n0++) { tmp = ar[n1][n0]*wr[n1*n0]/MAG-ai[n1][n0]*wi[n1*n0]/MAG; ai[n1][n0] = ar[n1][n0]*wi[n1*n0]/MAG+ai[n1][n0]*wr[n1*n0]/MAG; ar[n1][n0] = tmp; } /* cout<<"乘以旋转因子后序列："<<endl; for(n1=0;n1<R;n1++) { for(n0=0;n0<P;n0++) { cout<<"ax["<<n1<<"]["<<n0<<"]="<<ar[n1][n0]<<"+j*"<<ai[n1][n0]<<endl; } } cout<<"------------------------------------------------"<<endl; */ //对n0部分进行3进制倒位序排序 for(n1=0;n1<R;n1++) { for(i=0;i<P;i++) { k=i; j=0; t=log(P*1.00)/log(3.00); while((t--)>0) { j = j*3; j += k%3; k = k/3; } if(j>i) { tmp = ar[n1][j]; ar[n1][j] = ar[n1][i]; ar[n1][i] = tmp; tmp = ai[n1][j]; ai[n1][j] = ai[n1][i]; ai[n1][i] = tmp; } } } /* cout<<"3进制倒位序排序后的序列："<<endl; for(n1=0;n1<R;n1++) { for(n0=0;n0<P;n0++) { cout<<"ax["<<n1<<"]["<<n0<<"]="<<ar[n1][n0]<<"+j*"<<ai[n1][n0]<<endl; } } cout<<"------------------------------------------------"<<endl; */ //基三fft； for(n1=0;n1<R;n1++) { for(i=0;i<log(P*1.00)/log(3.00);i++) { l=pow(3,i); for(j=0;j<P;j+=3*l) { for(k=0;k<l;k++) { tmp1r = ar[n1][j+k+l]*w3r[P*k/3/l]/MAG - ai[n1][j+k+l]*w3i[P*k/3/l]/MAG; tmp1i = ar[n1][j+k+l]*w3i[P*k/3/l]/MAG + ai[n1][j+k+l]*w3r[P*k/3/l]/MAG; tmp2r = ar[n1][j+k+2*l]*w3r[2*P*k/3/l]/MAG - ai[n1][j+k+2*l]*w3i[2*P*k/3/l]/MAG; tmp2i = ar[n1][j+k+2*l]*w3i[2*P*k/3/l]/MAG + ai[n1][j+k+2*l]*w3r[2*P*k/3/l]/MAG; upr = ar[n1][j+k] + tmp1r + tmp2r; upi = ai[n1][j+k] + tmp1i + tmp2i; midr = ar[n1][j+k] + tmp1r*par1r/MAG - tmp1i*par1i/MAG + tmp2r*par2r/MAG - tmp2i*par2i/MAG; midi = ai[n1][j+k] + tmp1r*par1i/MAG + tmp1i*par1r/MAG + tmp2r*par2i/MAG + tmp2i*par2r/MAG; downr = ar[n1][j+k] + tmp1r*par2r/MAG - tmp1i*par2i/MAG + tmp2r*par1r/MAG - tmp2i*par1i/MAG; downi = ai[n1][j+k] + tmp1r*par2i/MAG + tmp1i*par2r/MAG + tmp2r*par1i/MAG + tmp2i*par1r/MAG ; ar[n1][j+k] = upr*sq3/MAG; ai[n1][j+k] = upi*sq3/MAG; ar[n1][j+k+l] = midr*sq3/MAG; ai[n1][j+k+l] = midi*sq3/MAG; ar[n1][j+k+2*l] = downr*sq3/MAG; ai[n1][j+k+2*l] = downi*sq3/MAG; } } } } /* cout<<"基三fft后的序列"<<endl; for(n1=0;n1<R;n1++) { for(n0=0;n0<P;n0++) { cout<<"ax["<<n1<<"]["<<n0<<"]="<<ar[n1][n0]<<"+j*"<<ai[n1][n0]<<endl; } } cout<<"------------------------------------------------"<<endl; */ //最后整序； for(n1=0;n1<R;n1++) for(n0=0;n0<P;n0++) { mr[R*n0+n1] = ar[n1][n0]; mi[R*n0+n1] = ai[n1][n0]; } /* cout<<"最后的序列："<<endl; for(n=0;n<POINT;n++) cout<<"mx["<<n<<"]="<<setprecision(10)<<mr[n]<<"+j*"<<mi[n]<<endl; cout<<"------------------------------------------------"<<endl; */ // cout<<"fft最后输出序列"<<endl; for(n=0;n<POINT;n++) { ixr[n] = mr[n]; ixi[n] = mi[n]; // cout<<"x["<<n<<"]="<<ixr[n]<<"+j*"<<ixi[n]<<endl; oxr[n] = (double)ixr[n]/MAG; oxi[n] = (double)ixi[n]/MAG; // cout<<"ox["<<n<<"]="<<oxr[n]<<"+j*"<<oxi[n]<<endl; } // cout<<"++++++++++++++++++++++++++++++++++++++++++++++++"<<endl; delete []ar; delete []ai; }