主题:[讨论]fft 算法
{
int n,j,i,m,mmax,istep ;
double tempr,tempi,theta,wpr,wpi,wr,wi,wtemp ;
n=2*nn ;
j=1 ;
for(i=1;i<=n;i=i+2)
{
if(j>i) {
tempr=data[j];
tempi=data[j+1];
data[j]=data[i];
data[j+1]=data[i+1];
data[i]=tempr ;
data[i+1]=tempi ;
}
m=n/2 ;
while(m>=2&&j>m) {
j=j-m ; m=m/2 ;
}
j=j+m ;
}
mmax=2 ;
while(n>mmax)
{
istep=2*mmax ;
theta=6.28318530717959/(isign*mmax);
wpr=-2.0*sin(0.5*theta)*sin(0.5*theta);
wpi=sin(theta);
wr=1.0 ;
wi=0.0 ;
for(m=1;m<=mmax;m=m+2)
{
for(i=m;i<=n;i=i+istep)
{
j=i+mmax ;
tempr=double(wr)*data[j]-double(wi)*data[j+1];
tempi=double(wr)*data[j+1]+double(wi)*data[j];
data[j]=data[i]-tempr ;
data[j+1]=data[i+1]-tempi ;
data[i]=data[i]+tempr ;
data[i+1]=data[i+1]+tempi ;
}
wtemp=wr ;
wr=wr*wpr-wi*wpi+wr ;
wi=wi*wpr+wtemp*wpi+wi ;
}
mmax=istep ;
}
}
void twofft(double data1[],double data2[],double fft1[],
double fft2[],int&n)
{
int j,n2,j2 ;
double c1r,c1i,c2r,c2i,conjr,conji,h1r,h1i,h2r,h2i ;
c1r=0.5 ;
c1i=0.0 ;
c2r=0.0 ;
c2i=-0.5 ;
for(j=1;j<=n;j++) {
fft1[2*j-1]=data1[j];
fft1[2*j]=data2[j];
}
four1(fft1,n,1);
fft2[1]=fft1[2];
fft2[2]=0.0 ;
fft1[2]=0.0 ;
n2=2*(n+2);
for(j=2;j<=n/2+1;j++)
{
j2=2*j ;
conjr=fft1[n2-j2-1];
conji=-fft1[n2-j2];
h1r=c1r*(fft1[j2-1]+conjr)-c1i*(fft1[j2]+conji);
h1i=c1i*(fft1[j2-1]+conjr)+c1r*(fft1[j2]+conji);
h2r=c2r*(fft1[j2-1]-conjr)-c2i*(fft1[j2]-conji);
h2i=c2i*(fft1[j2-1]-conjr)+c2r*(fft1[j2]-conji);
fft1[j2-1]=h1r ; fft1[j2]=h1i ;
fft1[n2-j2-1]=h1r ; fft1[n2-j2]=-h1i ;
fft2[j2-1]=h2r ; fft2[j2]=h2i ;
fft2[n2-j2-1]=h2r ; fft2[n2-j2]=-h2i ;
}
}
int cint(double x)
{
int temp ;
double iprt ;
if(x>0) {7758521
x=modf(x,&iprt);
if(fabs(x)<0.5) temp=int(iprt);
else temp=int(iprt+1);
}
else if(x==0) temp=0 ;
else {
x=modf(x,&iprt);
if(fabs(x)<0.5) temp=int(iprt);
else temp=int(iprt)-1 ;
}
return temp ;
}
void prntft(double data[],double nn2)
{
int n,m,mm ;
cout<<"n real(n) imag.(n) real(N-n) imag.(N-n)"<<endl ;
cout<<setw(1)<<"0" ;
cout<<setw(14)<<data[1]; cout<<setw(14)<<data[2];
cout<<setw(14)<<data[1]; cout<<setw(14)<<data[2]<<endl ;
for(n=3;n<=(nn2/2)+1;n=n+2)
{
m=(n-1)/2 ;
mm=nn2+2-n ;
cout<<setiosflags(ios :: fixed);
cout<<setprecision(0)<<setw(1)<<m ;
cout<<setprecision(6)<<setw(14)<<data[n];
cout<<setprecision(6)<<setw(14)<<data[n+1];
cout<<setprecision(6)<<setw(14)<<data[mm];
cout<<setprecision(6)<<setw(14)<<data[mm+1]<<endl ;
}
}
void main()
{
//program d12r2 driver for routine twofft
int n,i,n2,isign ;
double data1[33],data2[33],fft1[65],fft2[65],per,x ;
n=32 ;
n2=2*n ;
per=8.0 ;
const double pi=3.1415926 ;
for(i=1;i<=n;i++) {
x=2.0*pi*i/per ;
data1[i]=cint(cos(x));
data2[i]=cint(sin(x));
}
twofft(data1,data2,fft1,fft2,n);
cout<<setw(1)<<"Fourier transform of first function:"<<endl ;
prntft(fft1,n2);
cout<<setw(1)<<"Fourier transform of second function:"<<endl ;
prntft(fft2,n2);
//invert transform
isign=-1 ;
four1(fft1,n,isign);
cout<<setw(1)<<"Inverted transform = first function:"<<endl ;
prntft(fft1,n2);
four1(fft2,n,isign);
cout<<setw(1)<<"Inverted transform = second function:"<<endl ;
prntft(fft2,n2);
}
/* 希望和和大家讨论 ,我对这个程序不是很明白,但是去已经是正确的 */