newfft.cpp

C++矩阵算法库

   M_OVER_2=M/2+1; MP=M*P; PP=P/2; PM=P-1;

   for (U=0; U<PP; U++)
   {
      ANGLE=TWOPI*Real(U+1)/Real(P);
      JJ=P-U-2;
      A[U]=cos(ANGLE); B[U]=sin(ANGLE);
      A[JJ]=A[U]; B[JJ]= -B[U];
   }

   for (U=1; U<=PP; U++)
   {
      for (V=1; V<=PP; V++)
         { JJ=U*V-U*V/P*P; AA[V-1][U-1]=A[JJ-1]; BB[V-1][U-1]=B[JJ-1]; }
   }

   for (J=0; J<M_OVER_2; J++)
   {
      NO_FOLD = (J==0 || 2*J==M);
      K0=J;
      ANGLE=TWOPI*Real(J)/Real(MP); ZERO=ANGLE==0.0;
      C[0]=cos(ANGLE); S[0]=sin(ANGLE);
      for (U=1; U<PM; U++)
      {
         C[U]=C[U-1]*C[0]-S[U-1]*S[0];
         S[U]=S[U-1]*C[0]+C[U-1]*S[0];
      }
      goto L700;
   L500:
      REPORT
      if (NO_FOLD) { REPORT goto L1500; }
      REPORT
      NO_FOLD=true; K0=M-J;
      for (U=0; U<PM; U++)
         { T=C[U]*A[U]+S[U]*B[U]; S[U]= -S[U]*A[U]+C[U]*B[U]; C[U]=T; }
   L700:
      REPORT
      for (K=K0; K<N; K+=MP)
      {
         XT=X[K]; YT=Y[K];
         for (U=1; U<=PP; U++)
         {
            RA[U-1]=XT; IA[U-1]=YT;
            RB[U-1]=0.0; IB[U-1]=0.0;
         }
         for (U=1; U<=PP; U++)
         {
            JJ=P-U;
            RS=X[K+M*U]+X[K+M*JJ]; IS=Y[K+M*U]+Y[K+M*JJ];
            RU=X[K+M*U]-X[K+M*JJ]; IU=Y[K+M*U]-Y[K+M*JJ];
            XT=XT+RS; YT=YT+IS;
            for (V=0; V<PP; V++)
            {
               RA[V]=RA[V]+RS*AA[V][U-1]; IA[V]=IA[V]+IS*AA[V][U-1];
               RB[V]=RB[V]+RU*BB[V][U-1]; IB[V]=IB[V]+IU*BB[V][U-1];
            }
         }
         X[K]=XT; Y[K]=YT;
         for (U=1; U<=PP; U++)
         {
            if (!ZERO)
            {
               REPORT
               XT=RA[U-1]+IB[U-1]; YT=IA[U-1]-RB[U-1];
               X[K+M*U]=XT*C[U-1]+YT*S[U-1]; Y[K+M*U]=YT*C[U-1]-XT*S[U-1];
               JJ=P-U;
               XT=RA[U-1]-IB[U-1]; YT=IA[U-1]+RB[U-1];
               X[K+M*JJ]=XT*C[JJ-1]+YT*S[JJ-1];
               Y[K+M*JJ]=YT*C[JJ-1]-XT*S[JJ-1];
            }
            else
            {
               REPORT
               X[K+M*U]=RA[U-1]+IB[U-1]; Y[K+M*U]=IA[U-1]-RB[U-1];
               JJ=P-U;
               X[K+M*JJ]=RA[U-1]-IB[U-1]; Y[K+M*JJ]=IA[U-1]+RB[U-1];
            }
         }
      }
      goto L500;
L1500: ;
   }
   return;
}

static void R_2_FTK (int N, int M, Real* X0, Real* Y0, Real* X1, Real* Y1)
//    RADIX TWO FOURIER TRANSFORM KERNEL;
{
   REPORT
   bool NO_FOLD,ZERO;
   int  J,K,K0,M2,M_OVER_2;
   Real ANGLE,C,IS,IU,RS,RU,S,TWOPI;

   M2=M*2; M_OVER_2=M/2+1;
   TWOPI=8.0*atan(1.0);

   for (J=0; J<M_OVER_2; J++)
   {
      NO_FOLD = (J==0 || 2*J==M);
      K0=J;
      ANGLE=TWOPI*Real(J)/Real(M2); ZERO=ANGLE==0.0;
      C=cos(ANGLE); S=sin(ANGLE);
      goto L200;
   L100:
      REPORT
      if (NO_FOLD) { REPORT goto L600; }
      REPORT
      NO_FOLD=true; K0=M-J; C= -C;
   L200:
      REPORT
      for (K=K0; K<N; K+=M2)
      {
         RS=X0[K]+X1[K]; IS=Y0[K]+Y1[K];
         RU=X0[K]-X1[K]; IU=Y0[K]-Y1[K];
         X0[K]=RS; Y0[K]=IS;
         if (!ZERO) { X1[K]=RU*C+IU*S; Y1[K]=IU*C-RU*S; }
         else { X1[K]=RU; Y1[K]=IU; }
      }
      goto L100;
   L600: ;
   }

   return;
}

static void R_3_FTK (int N, int M, Real* X0, Real* Y0, Real* X1, Real* Y1,
   Real* X2, Real* Y2)
//    RADIX THREE FOURIER TRANSFORM KERNEL
{
   REPORT
   bool NO_FOLD,ZERO;
   int  J,K,K0,M3,M_OVER_2;
   Real ANGLE,A,B,C1,C2,S1,S2,T,TWOPI;
   Real I0,I1,I2,IA,IB,IS,R0,R1,R2,RA,RB,RS;

   M3=M*3; M_OVER_2=M/2+1; TWOPI=8.0*atan(1.0);
   A=cos(TWOPI/3.0); B=sin(TWOPI/3.0);

   for (J=0; J<M_OVER_2; J++)
   {
      NO_FOLD = (J==0 || 2*J==M);
      K0=J;
      ANGLE=TWOPI*Real(J)/Real(M3); ZERO=ANGLE==0.0;
      C1=cos(ANGLE); S1=sin(ANGLE);
      C2=C1*C1-S1*S1; S2=S1*C1+C1*S1;
      goto L200;
   L100:
      REPORT
      if (NO_FOLD) { REPORT goto L600; }
      REPORT
      NO_FOLD=true; K0=M-J;
      T=C1*A+S1*B; S1=C1*B-S1*A; C1=T;
      T=C2*A-S2*B; S2= -C2*B-S2*A; C2=T;
   L200:
      REPORT
      for (K=K0; K<N; K+=M3)
      {
         R0 = X0[K]; I0 = Y0[K];
         RS=X1[K]+X2[K]; IS=Y1[K]+Y2[K];
         X0[K]=R0+RS; Y0[K]=I0+IS;
         RA=R0+RS*A; IA=I0+IS*A;
         RB=(X1[K]-X2[K])*B; IB=(Y1[K]-Y2[K])*B;
         if (!ZERO)
         {
            REPORT
            R1=RA+IB; I1=IA-RB; R2=RA-IB; I2=IA+RB;
            X1[K]=R1*C1+I1*S1; Y1[K]=I1*C1-R1*S1;
            X2[K]=R2*C2+I2*S2; Y2[K]=I2*C2-R2*S2;
         }
         else { REPORT X1[K]=RA+IB; Y1[K]=IA-RB; X2[K]=RA-IB; Y2[K]=IA+RB; }
      }
      goto L100;
   L600: ;
   }

   return;
}

static void R_4_FTK (int N, int M,
   Real* X0, Real* Y0, Real* X1, Real* Y1,
   Real* X2, Real* Y2, Real* X3, Real* Y3)
//    RADIX FOUR FOURIER TRANSFORM KERNEL
{
   REPORT
   bool NO_FOLD,ZERO;
   int  J,K,K0,M4,M_OVER_2;
   Real ANGLE,C1,C2,C3,S1,S2,S3,T,TWOPI;
   Real I1,I2,I3,IS0,IS1,IU0,IU1,R1,R2,R3,RS0,RS1,RU0,RU1;

   M4=M*4; M_OVER_2=M/2+1;
   TWOPI=8.0*atan(1.0);

   for (J=0; J<M_OVER_2; J++)
   {
      NO_FOLD = (J==0 || 2*J==M);
      K0=J;
      ANGLE=TWOPI*Real(J)/Real(M4); ZERO=ANGLE==0.0;
      C1=cos(ANGLE); S1=sin(ANGLE);
      C2=C1*C1-S1*S1; S2=S1*C1+C1*S1;
      C3=C2*C1-S2*S1; S3=S2*C1+C2*S1;
      goto L200;
   L100:
      REPORT
      if (NO_FOLD) { REPORT goto L600; }
      REPORT
      NO_FOLD=true; K0=M-J;
      T=C1; C1=S1; S1=T;
      C2= -C2;
      T=C3; C3= -S3; S3= -T;
   L200:
      REPORT
      for (K=K0; K<N; K+=M4)
      {
         RS0=X0[K]+X2[K]; IS0=Y0[K]+Y2[K];
         RU0=X0[K]-X2[K]; IU0=Y0[K]-Y2[K];
         RS1=X1[K]+X3[K]; IS1=Y1[K]+Y3[K];
         RU1=X1[K]-X3[K]; IU1=Y1[K]-Y3[K];
         X0[K]=RS0+RS1; Y0[K]=IS0+IS1;
         if (!ZERO)
         {
            REPORT
            R1=RU0+IU1; I1=IU0-RU1;
            R2=RS0-RS1; I2=IS0-IS1;
            R3=RU0-IU1; I3=IU0+RU1;
            X2[K]=R1*C1+I1*S1; Y2[K]=I1*C1-R1*S1;
            X1[K]=R2*C2+I2*S2; Y1[K]=I2*C2-R2*S2;
            X3[K]=R3*C3+I3*S3; Y3[K]=I3*C3-R3*S3;
         }
         else
         {
            REPORT
            X2[K]=RU0+IU1; Y2[K]=IU0-RU1;
            X1[K]=RS0-RS1; Y1[K]=IS0-IS1;
            X3[K]=RU0-IU1; Y3[K]=IU0+RU1;
         }
      }
      goto L100;
   L600: ;
   }

   return;
}

static void R_5_FTK (int N, int M,
   Real* X0, Real* Y0, Real* X1, Real* Y1, Real* X2, Real* Y2,
   Real* X3, Real* Y3, Real* X4, Real* Y4)
//    RADIX FIVE FOURIER TRANSFORM KERNEL

{
   REPORT
   bool NO_FOLD,ZERO;
   int  J,K,K0,M5,M_OVER_2;
   Real ANGLE,A1,A2,B1,B2,C1,C2,C3,C4,S1,S2,S3,S4,T,TWOPI;
   Real R0,R1,R2,R3,R4,RA1,RA2,RB1,RB2,RS1,RS2,RU1,RU2;
   Real I0,I1,I2,I3,I4,IA1,IA2,IB1,IB2,IS1,IS2,IU1,IU2;

   M5=M*5; M_OVER_2=M/2+1;
   TWOPI=8.0*atan(1.0);
   A1=cos(TWOPI/5.0); B1=sin(TWOPI/5.0);
   A2=cos(2.0*TWOPI/5.0); B2=sin(2.0*TWOPI/5.0);

   for (J=0; J<M_OVER_2; J++)
   {
      NO_FOLD = (J==0 || 2*J==M);
      K0=J;
      ANGLE=TWOPI*Real(J)/Real(M5); ZERO=ANGLE==0.0;
      C1=cos(ANGLE); S1=sin(ANGLE);
      C2=C1*C1-S1*S1; S2=S1*C1+C1*S1;
      C3=C2*C1-S2*S1; S3=S2*C1+C2*S1;
      C4=C2*C2-S2*S2; S4=S2*C2+C2*S2;
      goto L200;
   L100:
      REPORT
      if (NO_FOLD) { REPORT goto L600; }
      REPORT
      NO_FOLD=true; K0=M-J;
      T=C1*A1+S1*B1; S1=C1*B1-S1*A1; C1=T;
      T=C2*A2+S2*B2; S2=C2*B2-S2*A2; C2=T;
      T=C3*A2-S3*B2; S3= -C3*B2-S3*A2; C3=T;
      T=C4*A1-S4*B1; S4= -C4*B1-S4*A1; C4=T;
   L200:
      REPORT
      for (K=K0; K<N; K+=M5)
      {
         R0=X0[K]; I0=Y0[K];
         RS1=X1[K]+X4[K]; IS1=Y1[K]+Y4[K];
         RU1=X1[K]-X4[K]; IU1=Y1[K]-Y4[K];
         RS2=X2[K]+X3[K]; IS2=Y2[K]+Y3[K];
         RU2=X2[K]-X3[K]; IU2=Y2[K]-Y3[K];
         X0[K]=R0+RS1+RS2; Y0[K]=I0+IS1+IS2;
         RA1=R0+RS1*A1+RS2*A2; IA1=I0+IS1*A1+IS2*A2;
         RA2=R0+RS1*A2+RS2*A1; IA2=I0+IS1*A2+IS2*A1;
         RB1=RU1*B1+RU2*B2; IB1=IU1*B1+IU2*B2;
         RB2=RU1*B2-RU2*B1; IB2=IU1*B2-IU2*B1;
         if (!ZERO)
         {
            REPORT
            R1=RA1+IB1; I1=IA1-RB1;
            R2=RA2+IB2; I2=IA2-RB2;
            R3=RA2-IB2; I3=IA2+RB2;
            R4=RA1-IB1; I4=IA1+RB1;
            X1[K]=R1*C1+I1*S1; Y1[K]=I1*C1-R1*S1;
            X2[K]=R2*C2+I2*S2; Y2[K]=I2*C2-R2*S2;
            X3[K]=R3*C3+I3*S3; Y3[K]=I3*C3-R3*S3;
            X4[K]=R4*C4+I4*S4; Y4[K]=I4*C4-R4*S4;
         }
         else
         {
            REPORT
            X1[K]=RA1+IB1; Y1[K]=IA1-RB1;
            X2[K]=RA2+IB2; Y2[K]=IA2-RB2;
            X3[K]=RA2-IB2; Y3[K]=IA2+RB2;
            X4[K]=RA1-IB1; Y4[K]=IA1+RB1;
         }
      }
      goto L100;
   L600: ;
   }

   return;
}

static void R_8_FTK (int N, int M,
   Real* X0, Real* Y0, Real* X1, Real* Y1,
   Real* X2, Real* Y2, Real* X3, Real* Y3,
   Real* X4, Real* Y4, Real* X5, Real* Y5,
   Real* X6, Real* Y6, Real* X7, Real* Y7)
//    RADIX EIGHT FOURIER TRANSFORM KERNEL
{
   REPORT
   bool NO_FOLD,ZERO;
   int  J,K,K0,M8,M_OVER_2;
   Real ANGLE,C1,C2,C3,C4,C5,C6,C7,E,S1,S2,S3,S4,S5,S6,S7,T,TWOPI;
   Real R1,R2,R3,R4,R5,R6,R7,RS0,RS1,RS2,RS3,RU0,RU1,RU2,RU3;
   Real I1,I2,I3,I4,I5,I6,I7,IS0,IS1,IS2,IS3,IU0,IU1,IU2,IU3;
   Real RSS0,RSS1,RSU0,RSU1,RUS0,RUS1,RUU0,RUU1;
   Real ISS0,ISS1,ISU0,ISU1,IUS0,IUS1,IUU0,IUU1;

   M8=M*8; M_OVER_2=M/2+1;
   TWOPI=8.0*atan(1.0); E=cos(TWOPI/8.0);

   for (J=0;J<M_OVER_2;J++)
   {
      NO_FOLD= (J==0 || 2*J==M);
      K0=J;
      ANGLE=TWOPI*Real(J)/Real(M8); ZERO=ANGLE==0.0;
      C1=cos(ANGLE); S1=sin(ANGLE);
      C2=C1*C1-S1*S1; S2=C1*S1+S1*C1;
      C3=C2*C1-S2*S1; S3=S2*C1+C2*S1;
      C4=C2*C2-S2*S2; S4=S2*C2+C2*S2;
      C5=C4*C1-S4*S1; S5=S4*C1+C4*S1;
      C6=C4*C2-S4*S2; S6=S4*C2+C4*S2;
      C7=C4*C3-S4*S3; S7=S4*C3+C4*S3;
      goto L200;