fourl.cpp

本课程提供的示范程序包括三个方面

#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr;


void fourl(data,nn,isign)
float data[];
int nn,isign;
{
 int n,mmax,m,j,istep,i;
 float wtemp,wr,wpr,wpi,wi,theta;
 float tempr,tempi;
 n=nn<<1;
 j=1;
 for(i=1;i<n;i+=2) 
 {
  if(j>i)
    {
     SWAP(data[j],data[i]);
     SWAP(data[j+1],data[i+1]);
    }
  m=n>>1;
   while(m>=2&&j>m)
   {
    j-=m;
    m>>=1;
    }
   j+=m;
 }

  mmax=2;
  while(n>mmax)
  {
   istep=2*mmax;
   theta=6.28318530717959/(isign*mmax);
   wtemp=sin(0.5*theta);
   wpr=-2.0*wtemp*wtemp;
   wpi=sin(theta);
   wr=1.0;
   wi=0.0;
   for(m=1;m<mmax;m+=2)
    {
     for(i=m;i<=n;i+=istep)
      {
       j=i+mmax;
       tempr=wr*data[j]-wi*data[j+1];
       tempi=wr*data[j+1]+wi*data[j];
       data[j]=data[i]-tempr;
       data[j+1]=data[i+1]-tempi;
       data[i]+=tempr;
       data[i+1]+=tempi;
      }
      wr=(wtemp=wr)*wpr-wi*wpi+wr;
      wi=wi*wpr+wtemp*wpi+wi;
     }
     mmax=istep;
    }
   if(isign==-1)
   for(i=1;i<=2*nn;i++)
   data[i]/=nn;
  }