newfft.cpp

C++矩阵算法库

   L100:
      REPORT
      if (NO_FOLD) { REPORT goto L600; }
      REPORT
      NO_FOLD=true; K0=M-J;
      T=(C1+S1)*E; S1=(C1-S1)*E; C1=T;
      T=S2; S2=C2; C2=T;
      T=(-C3+S3)*E; S3=(C3+S3)*E; C3=T;
      C4= -C4;
      T= -(C5+S5)*E; S5=(-C5+S5)*E; C5=T;
      T= -S6; S6= -C6; C6=T;
      T=(C7-S7)*E; S7= -(C7+S7)*E; C7=T;
   L200:
      REPORT
      for (K=K0; K<N; K+=M8)
      {
         RS0=X0[K]+X4[K]; IS0=Y0[K]+Y4[K];
         RU0=X0[K]-X4[K]; IU0=Y0[K]-Y4[K];
         RS1=X1[K]+X5[K]; IS1=Y1[K]+Y5[K];
         RU1=X1[K]-X5[K]; IU1=Y1[K]-Y5[K];
         RS2=X2[K]+X6[K]; IS2=Y2[K]+Y6[K];
         RU2=X2[K]-X6[K]; IU2=Y2[K]-Y6[K];
         RS3=X3[K]+X7[K]; IS3=Y3[K]+Y7[K];
         RU3=X3[K]-X7[K]; IU3=Y3[K]-Y7[K];
         RSS0=RS0+RS2; ISS0=IS0+IS2;
         RSU0=RS0-RS2; ISU0=IS0-IS2;
         RSS1=RS1+RS3; ISS1=IS1+IS3;
         RSU1=RS1-RS3; ISU1=IS1-IS3;
         RUS0=RU0-IU2; IUS0=IU0+RU2;
         RUU0=RU0+IU2; IUU0=IU0-RU2;
         RUS1=RU1-IU3; IUS1=IU1+RU3;
         RUU1=RU1+IU3; IUU1=IU1-RU3;
         T=(RUS1+IUS1)*E; IUS1=(IUS1-RUS1)*E; RUS1=T;
         T=(RUU1+IUU1)*E; IUU1=(IUU1-RUU1)*E; RUU1=T;
         X0[K]=RSS0+RSS1; Y0[K]=ISS0+ISS1;
         if (!ZERO)
         {
            REPORT
            R1=RUU0+RUU1; I1=IUU0+IUU1;
            R2=RSU0+ISU1; I2=ISU0-RSU1;
            R3=RUS0+IUS1; I3=IUS0-RUS1;
            R4=RSS0-RSS1; I4=ISS0-ISS1;
            R5=RUU0-RUU1; I5=IUU0-IUU1;
            R6=RSU0-ISU1; I6=ISU0+RSU1;
            R7=RUS0-IUS1; I7=IUS0+RUS1;
            X4[K]=R1*C1+I1*S1; Y4[K]=I1*C1-R1*S1;
            X2[K]=R2*C2+I2*S2; Y2[K]=I2*C2-R2*S2;
            X6[K]=R3*C3+I3*S3; Y6[K]=I3*C3-R3*S3;
            X1[K]=R4*C4+I4*S4; Y1[K]=I4*C4-R4*S4;
            X5[K]=R5*C5+I5*S5; Y5[K]=I5*C5-R5*S5;
            X3[K]=R6*C6+I6*S6; Y3[K]=I6*C6-R6*S6;
            X7[K]=R7*C7+I7*S7; Y7[K]=I7*C7-R7*S7;
         }
         else
         {
            REPORT
            X4[K]=RUU0+RUU1; Y4[K]=IUU0+IUU1;
            X2[K]=RSU0+ISU1; Y2[K]=ISU0-RSU1;
            X6[K]=RUS0+IUS1; Y6[K]=IUS0-RUS1;
            X1[K]=RSS0-RSS1; Y1[K]=ISS0-ISS1;
            X5[K]=RUU0-RUU1; Y5[K]=IUU0-IUU1;
            X3[K]=RSU0-ISU1; Y3[K]=ISU0+RSU1;
            X7[K]=RUS0-IUS1; Y7[K]=IUS0+RUS1;
         }
      }
      goto L100;
   L600: ;
   }

   return;
}

static void R_16_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,
   Real* X8, Real* Y8, Real* X9, Real* Y9,
   Real* X10, Real* Y10, Real* X11, Real* Y11,
   Real* X12, Real* Y12, Real* X13, Real* Y13,
   Real* X14, Real* Y14, Real* X15, Real* Y15)
//    RADIX SIXTEEN FOURIER TRANSFORM KERNEL
{
   REPORT
   bool NO_FOLD,ZERO;
   int  J,K,K0,M16,M_OVER_2;
   Real ANGLE,EI1,ER1,E2,EI3,ER3,EI5,ER5,T,TWOPI;
   Real RS0,RS1,RS2,RS3,RS4,RS5,RS6,RS7;
   Real IS0,IS1,IS2,IS3,IS4,IS5,IS6,IS7;
   Real RU0,RU1,RU2,RU3,RU4,RU5,RU6,RU7;
   Real IU0,IU1,IU2,IU3,IU4,IU5,IU6,IU7;
   Real RUS0,RUS1,RUS2,RUS3,RUU0,RUU1,RUU2,RUU3;
   Real ISS0,ISS1,ISS2,ISS3,ISU0,ISU1,ISU2,ISU3;
   Real RSS0,RSS1,RSS2,RSS3,RSU0,RSU1,RSU2,RSU3;
   Real IUS0,IUS1,IUS2,IUS3,IUU0,IUU1,IUU2,IUU3;
   Real RSSS0,RSSS1,RSSU0,RSSU1,RSUS0,RSUS1,RSUU0,RSUU1;
   Real ISSS0,ISSS1,ISSU0,ISSU1,ISUS0,ISUS1,ISUU0,ISUU1;
   Real RUSS0,RUSS1,RUSU0,RUSU1,RUUS0,RUUS1,RUUU0,RUUU1;
   Real IUSS0,IUSS1,IUSU0,IUSU1,IUUS0,IUUS1,IUUU0,IUUU1;
   Real R1,R2,R3,R4,R5,R6,R7,R8,R9,R10,R11,R12,R13,R14,R15;
   Real I1,I2,I3,I4,I5,I6,I7,I8,I9,I10,I11,I12,I13,I14,I15;
   Real C1,C2,C3,C4,C5,C6,C7,C8,C9,C10,C11,C12,C13,C14,C15;
   Real S1,S2,S3,S4,S5,S6,S7,S8,S9,S10,S11,S12,S13,S14,S15;

   M16=M*16; M_OVER_2=M/2+1;
   TWOPI=8.0*atan(1.0);
   ER1=cos(TWOPI/16.0); EI1=sin(TWOPI/16.0);
   E2=cos(TWOPI/8.0);
   ER3=cos(3.0*TWOPI/16.0); EI3=sin(3.0*TWOPI/16.0);
   ER5=cos(5.0*TWOPI/16.0); EI5=sin(5.0*TWOPI/16.0);

   for (J=0; J<M_OVER_2; J++)
   {
      NO_FOLD = (J==0 || 2*J==M);
      K0=J;
      ANGLE=TWOPI*Real(J)/Real(M16);
      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;
      C8=C4*C4-S4*S4; S8=C4*S4+S4*C4;
      C9=C8*C1-S8*S1; S9=S8*C1+C8*S1;
      C10=C8*C2-S8*S2; S10=S8*C2+C8*S2;
      C11=C8*C3-S8*S3; S11=S8*C3+C8*S3;
      C12=C8*C4-S8*S4; S12=S8*C4+C8*S4;
      C13=C8*C5-S8*S5; S13=S8*C5+C8*S5;
      C14=C8*C6-S8*S6; S14=S8*C6+C8*S6;
      C15=C8*C7-S8*S7; S15=S8*C7+C8*S7;
      goto L200;
   L100:
      REPORT
      if (NO_FOLD) { REPORT goto L600; }
      REPORT
      NO_FOLD=true; K0=M-J;
      T=C1*ER1+S1*EI1; S1= -S1*ER1+C1*EI1; C1=T;
      T=(C2+S2)*E2; S2=(C2-S2)*E2; C2=T;
      T=C3*ER3+S3*EI3; S3= -S3*ER3+C3*EI3; C3=T;
      T=S4; S4=C4; C4=T;
      T=S5*ER1-C5*EI1; S5=C5*ER1+S5*EI1; C5=T;
      T=(-C6+S6)*E2; S6=(C6+S6)*E2; C6=T;
      T=S7*ER3-C7*EI3; S7=C7*ER3+S7*EI3; C7=T;
      C8= -C8;
      T= -(C9*ER1+S9*EI1); S9=S9*ER1-C9*EI1; C9=T;
      T= -(C10+S10)*E2; S10=(-C10+S10)*E2; C10=T;
      T= -(C11*ER3+S11*EI3); S11=S11*ER3-C11*EI3; C11=T;
      T= -S12; S12= -C12; C12=T;
      T= -S13*ER1+C13*EI1; S13= -(C13*ER1+S13*EI1); C13=T;
      T=(C14-S14)*E2; S14= -(C14+S14)*E2; C14=T;
      T= -S15*ER3+C15*EI3; S15= -(C15*ER3+S15*EI3); C15=T;
   L200:
      REPORT
      for (K=K0; K<N; K+=M16)
      {
         RS0=X0[K]+X8[K]; IS0=Y0[K]+Y8[K];
         RU0=X0[K]-X8[K]; IU0=Y0[K]-Y8[K];
         RS1=X1[K]+X9[K]; IS1=Y1[K]+Y9[K];
         RU1=X1[K]-X9[K]; IU1=Y1[K]-Y9[K];
         RS2=X2[K]+X10[K]; IS2=Y2[K]+Y10[K];
         RU2=X2[K]-X10[K]; IU2=Y2[K]-Y10[K];
         RS3=X3[K]+X11[K]; IS3=Y3[K]+Y11[K];
         RU3=X3[K]-X11[K]; IU3=Y3[K]-Y11[K];
         RS4=X4[K]+X12[K]; IS4=Y4[K]+Y12[K];
         RU4=X4[K]-X12[K]; IU4=Y4[K]-Y12[K];
         RS5=X5[K]+X13[K]; IS5=Y5[K]+Y13[K];
         RU5=X5[K]-X13[K]; IU5=Y5[K]-Y13[K];
         RS6=X6[K]+X14[K]; IS6=Y6[K]+Y14[K];
         RU6=X6[K]-X14[K]; IU6=Y6[K]-Y14[K];
         RS7=X7[K]+X15[K]; IS7=Y7[K]+Y15[K];
         RU7=X7[K]-X15[K]; IU7=Y7[K]-Y15[K];
         RSS0=RS0+RS4; ISS0=IS0+IS4;
         RSS1=RS1+RS5; ISS1=IS1+IS5;
         RSS2=RS2+RS6; ISS2=IS2+IS6;
         RSS3=RS3+RS7; ISS3=IS3+IS7;
         RSU0=RS0-RS4; ISU0=IS0-IS4;
         RSU1=RS1-RS5; ISU1=IS1-IS5;
         RSU2=RS2-RS6; ISU2=IS2-IS6;
         RSU3=RS3-RS7; ISU3=IS3-IS7;
         RUS0=RU0-IU4; IUS0=IU0+RU4;
         RUS1=RU1-IU5; IUS1=IU1+RU5;
         RUS2=RU2-IU6; IUS2=IU2+RU6;
         RUS3=RU3-IU7; IUS3=IU3+RU7;
         RUU0=RU0+IU4; IUU0=IU0-RU4;
         RUU1=RU1+IU5; IUU1=IU1-RU5;
         RUU2=RU2+IU6; IUU2=IU2-RU6;
         RUU3=RU3+IU7; IUU3=IU3-RU7;
         T=(RSU1+ISU1)*E2; ISU1=(ISU1-RSU1)*E2; RSU1=T;
         T=(RSU3+ISU3)*E2; ISU3=(ISU3-RSU3)*E2; RSU3=T;
         T=RUS1*ER3+IUS1*EI3; IUS1=IUS1*ER3-RUS1*EI3; RUS1=T;
         T=(RUS2+IUS2)*E2; IUS2=(IUS2-RUS2)*E2; RUS2=T;
         T=RUS3*ER5+IUS3*EI5; IUS3=IUS3*ER5-RUS3*EI5; RUS3=T;
         T=RUU1*ER1+IUU1*EI1; IUU1=IUU1*ER1-RUU1*EI1; RUU1=T;
         T=(RUU2+IUU2)*E2; IUU2=(IUU2-RUU2)*E2; RUU2=T;
         T=RUU3*ER3+IUU3*EI3; IUU3=IUU3*ER3-RUU3*EI3; RUU3=T;
         RSSS0=RSS0+RSS2; ISSS0=ISS0+ISS2;
         RSSS1=RSS1+RSS3; ISSS1=ISS1+ISS3;
         RSSU0=RSS0-RSS2; ISSU0=ISS0-ISS2;
         RSSU1=RSS1-RSS3; ISSU1=ISS1-ISS3;
         RSUS0=RSU0-ISU2; ISUS0=ISU0+RSU2;
         RSUS1=RSU1-ISU3; ISUS1=ISU1+RSU3;
         RSUU0=RSU0+ISU2; ISUU0=ISU0-RSU2;
         RSUU1=RSU1+ISU3; ISUU1=ISU1-RSU3;
         RUSS0=RUS0-IUS2; IUSS0=IUS0+RUS2;
         RUSS1=RUS1-IUS3; IUSS1=IUS1+RUS3;
         RUSU0=RUS0+IUS2; IUSU0=IUS0-RUS2;
         RUSU1=RUS1+IUS3; IUSU1=IUS1-RUS3;
         RUUS0=RUU0+RUU2; IUUS0=IUU0+IUU2;
         RUUS1=RUU1+RUU3; IUUS1=IUU1+IUU3;
         RUUU0=RUU0-RUU2; IUUU0=IUU0-IUU2;
         RUUU1=RUU1-RUU3; IUUU1=IUU1-IUU3;
         X0[K]=RSSS0+RSSS1; Y0[K]=ISSS0+ISSS1;
         if (!ZERO)
         {
            REPORT
            R1=RUUS0+RUUS1; I1=IUUS0+IUUS1;
            R2=RSUU0+RSUU1; I2=ISUU0+ISUU1;
            R3=RUSU0+RUSU1; I3=IUSU0+IUSU1;
            R4=RSSU0+ISSU1; I4=ISSU0-RSSU1;
            R5=RUUU0+IUUU1; I5=IUUU0-RUUU1;
            R6=RSUS0+ISUS1; I6=ISUS0-RSUS1;
            R7=RUSS0+IUSS1; I7=IUSS0-RUSS1;
            R8=RSSS0-RSSS1; I8=ISSS0-ISSS1;
            R9=RUUS0-RUUS1; I9=IUUS0-IUUS1;
            R10=RSUU0-RSUU1; I10=ISUU0-ISUU1;
            R11=RUSU0-RUSU1; I11=IUSU0-IUSU1;
            R12=RSSU0-ISSU1; I12=ISSU0+RSSU1;
            R13=RUUU0-IUUU1; I13=IUUU0+RUUU1;
            R14=RSUS0-ISUS1; I14=ISUS0+RSUS1;
            R15=RUSS0-IUSS1; I15=IUSS0+RUSS1;
            X8[K]=R1*C1+I1*S1; Y8[K]=I1*C1-R1*S1;
            X4[K]=R2*C2+I2*S2; Y4[K]=I2*C2-R2*S2;
            X12[K]=R3*C3+I3*S3; Y12[K]=I3*C3-R3*S3;
            X2[K]=R4*C4+I4*S4; Y2[K]=I4*C4-R4*S4;
            X10[K]=R5*C5+I5*S5; Y10[K]=I5*C5-R5*S5;
            X6[K]=R6*C6+I6*S6; Y6[K]=I6*C6-R6*S6;
            X14[K]=R7*C7+I7*S7; Y14[K]=I7*C7-R7*S7;
            X1[K]=R8*C8+I8*S8; Y1[K]=I8*C8-R8*S8;
            X9[K]=R9*C9+I9*S9; Y9[K]=I9*C9-R9*S9;
            X5[K]=R10*C10+I10*S10; Y5[K]=I10*C10-R10*S10;
            X13[K]=R11*C11+I11*S11; Y13[K]=I11*C11-R11*S11;
            X3[K]=R12*C12+I12*S12; Y3[K]=I12*C12-R12*S12;
            X11[K]=R13*C13+I13*S13; Y11[K]=I13*C13-R13*S13;
            X7[K]=R14*C14+I14*S14; Y7[K]=I14*C14-R14*S14;
            X15[K]=R15*C15+I15*S15; Y15[K]=I15*C15-R15*S15;
         }
         else
         {
            REPORT
            X8[K]=RUUS0+RUUS1; Y8[K]=IUUS0+IUUS1;
            X4[K]=RSUU0+RSUU1; Y4[K]=ISUU0+ISUU1;
            X12[K]=RUSU0+RUSU1; Y12[K]=IUSU0+IUSU1;
            X2[K]=RSSU0+ISSU1; Y2[K]=ISSU0-RSSU1;
            X10[K]=RUUU0+IUUU1; Y10[K]=IUUU0-RUUU1;
            X6[K]=RSUS0+ISUS1; Y6[K]=ISUS0-RSUS1;
            X14[K]=RUSS0+IUSS1; Y14[K]=IUSS0-RUSS1;
            X1[K]=RSSS0-RSSS1; Y1[K]=ISSS0-ISSS1;
            X9[K]=RUUS0-RUUS1; Y9[K]=IUUS0-IUUS1;
            X5[K]=RSUU0-RSUU1; Y5[K]=ISUU0-ISUU1;
            X13[K]=RUSU0-RUSU1; Y13[K]=IUSU0-IUSU1;
            X3[K]=RSSU0-ISSU1; Y3[K]=ISSU0+RSSU1;
            X11[K]=RUUU0-IUUU1; Y11[K]=IUUU0+RUUU1;
            X7[K]=RSUS0-ISUS1; Y7[K]=ISUS0+RSUS1;
            X15[K]=RUSS0-IUSS1; Y15[K]=IUSS0+RUSS1;
         }
      }
      goto L100;
   L600: ;
   }

   return;
}

// can the number of points be factorised sufficiently
// for the fft to run

bool FFT_Controller::CanFactor(int PTS)
{
   REPORT
   const int NP = 16, NQ = 10, PMAX=19;

   if (PTS<=1) { REPORT return true; }

   int N = PTS, F = 2, P = 0, Q = 0;

   while (N > 1)
   {
      bool fail = true;
      for (int J = F; J <= PMAX; J++)
         if (N % J == 0) { fail = false; F=J; break; }
      if (fail || P >= NP || Q >= NQ) { REPORT return false; }
      N /= F;
      if (N % F != 0) Q++; else { N /= F; P++; }
   }

   return true;    // can factorise

}

bool FFT_Controller::OnlyOldFFT;         // static variable

// **************************** multi radix counter **********************

MultiRadixCounter::MultiRadixCounter(int nx, const SimpleIntArray& rx,
   SimpleIntArray& vx)
   : Radix(rx), Value(vx), n(nx), reverse(0),
      product(1), counter(0), finish(false)
{
   REPORT for (int k = 0; k < n; k++) { Value[k] = 0; product *= Radix[k]; }
}

void MultiRadixCounter::operator++()
{
   REPORT
   counter++; int p = product;
   for (int k = 0; k < n; k++)
   {
      Value[k]++; int p1 = p / Radix[k]; reverse += p1;
      if (Value[k] == Radix[k]) { REPORT Value[k] = 0; reverse -= p; p = p1; }
      else { REPORT return; }
   }
   finish = true;
}


static int BitReverse(int x, int prod, int n, const SimpleIntArray& f)
{
   // x = c[0]+f[0]*(c[1]+f[1]*(c[2]+...
   // return c[n-1]+f[n-1]*(c[n-2]+f[n-2]*(c[n-3]+...
   // prod is the product of the f[i]
   // n is the number of f[i] (don't assume f has the correct length)

   REPORT
   const int* d = f.Data() + n; int sum = 0; int q = 1;
   while (n--)
   {
      prod /= *(--d);
      int c = x / prod; x-= c * prod;
      sum += q * c; q *= *d;
   }
   return sum;
}


#ifdef use_namespace
}
#endif