00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094 #ifndef WANT_STREAM
00095 #define WANT_STREAM
00096 #endif
00097
00098 #ifndef WANT_MATH
00099 #define WANT_MATH
00100 #endif
00101
00102 #include "newmatap.h"
00103
00104 #ifdef use_namespace
00105 namespace NEWMAT {
00106 #endif
00107
00108 #ifdef DO_REPORT
00109 #define REPORT { static ExeCounter ExeCount(__LINE__,20); ++ExeCount; }
00110 #else
00111 #define REPORT {}
00112 #endif
00113
00114 inline Real square(Real x) { return x*x; }
00115 inline int square(int x) { return x*x; }
00116
00117 static void GR_1D_FS (int PTS, int N_SYM, int N_UN_SYM,
00118 const SimpleIntArray& SYM, int P_SYM, const SimpleIntArray& UN_SYM,
00119 Real* X, Real* Y);
00120 static void GR_1D_FT (int N, int N_FACTOR, const SimpleIntArray& FACTOR,
00121 Real* X, Real* Y);
00122 static void R_P_FTK (int N, int M, int P, Real* X, Real* Y);
00123 static void R_2_FTK (int N, int M, Real* X0, Real* Y0, Real* X1, Real* Y1);
00124 static void R_3_FTK (int N, int M, Real* X0, Real* Y0, Real* X1, Real* Y1,
00125 Real* X2, Real* Y2);
00126 static void R_4_FTK (int N, int M,
00127 Real* X0, Real* Y0, Real* X1, Real* Y1,
00128 Real* X2, Real* Y2, Real* X3, Real* Y3);
00129 static void R_5_FTK (int N, int M,
00130 Real* X0, Real* Y0, Real* X1, Real* Y1, Real* X2, Real* Y2,
00131 Real* X3, Real* Y3, Real* X4, Real* Y4);
00132 static void R_8_FTK (int N, int M,
00133 Real* X0, Real* Y0, Real* X1, Real* Y1,
00134 Real* X2, Real* Y2, Real* X3, Real* Y3,
00135 Real* X4, Real* Y4, Real* X5, Real* Y5,
00136 Real* X6, Real* Y6, Real* X7, Real* Y7);
00137 static void R_16_FTK (int N, int M,
00138 Real* X0, Real* Y0, Real* X1, Real* Y1,
00139 Real* X2, Real* Y2, Real* X3, Real* Y3,
00140 Real* X4, Real* Y4, Real* X5, Real* Y5,
00141 Real* X6, Real* Y6, Real* X7, Real* Y7,
00142 Real* X8, Real* Y8, Real* X9, Real* Y9,
00143 Real* X10, Real* Y10, Real* X11, Real* Y11,
00144 Real* X12, Real* Y12, Real* X13, Real* Y13,
00145 Real* X14, Real* Y14, Real* X15, Real* Y15);
00146 static int BitReverse(int x, int prod, int n, const SimpleIntArray& f);
00147
00148
00149 bool FFT_Controller::ar_1d_ft (int PTS, Real* X, Real *Y)
00150 {
00151
00152
00153 REPORT
00154
00155 int F,J,N,NF,P,PMAX,P_SYM,P_TWO,Q,R,TWO_GRP;
00156
00157
00158
00159 const int NP = 16, NQ = 10;
00160 SimpleIntArray PP(NP), QQ(NQ);
00161
00162 TWO_GRP=16; PMAX=19;
00163
00164
00165
00166
00167
00168
00169 if (PTS<=1) return true;
00170 N=PTS; P_SYM=1; F=2; P=0; Q=0;
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186 while (N > 1)
00187 {
00188 bool fail = true;
00189 for (J=F; J<=PMAX; J++)
00190 if (N % J == 0) { fail = false; F=J; break; }
00191 if (fail || P >= NP || Q >= NQ) return false;
00192 N /= F;
00193 if (N % F != 0) QQ[Q++] = F;
00194 else { N /= F; PP[P++] = F; P_SYM *= F; }
00195 }
00196
00197 R = (Q == 0) ? 0 : 1;
00198
00199 NF = 2*P + Q;
00200 SimpleIntArray FACTOR(NF + 1), SYM(2*P + R);
00201 FACTOR[NF] = 0;
00202
00203
00204 for (J=0; J<P; J++)
00205 { SYM[J]=FACTOR[J]=PP[P-1-J]; FACTOR[P+Q+J]=SYM[P+R+J]=PP[J]; }
00206
00207 if (Q>0)
00208 {
00209 REPORT
00210 for (J=0; J<Q; J++) FACTOR[P+J]=QQ[J];
00211 SYM[P]=PTS/square(P_SYM);
00212 }
00213
00214
00215 P_TWO = 1;
00216 for (J=0; J < NF; J++)
00217 {
00218 if (FACTOR[J]!=2) continue;
00219 P_TWO=P_TWO*2; FACTOR[J]=1;
00220 if (P_TWO<TWO_GRP && FACTOR[J+1]==2) continue;
00221 FACTOR[J]=P_TWO; P_TWO=1;
00222 }
00223
00224 if (P==0) R=0;
00225 if (Q<=1) Q=0;
00226
00227
00228 GR_1D_FT(PTS,NF,FACTOR,X,Y);
00229 GR_1D_FS(PTS,2*P+R,Q,SYM,P_SYM,QQ,X,Y);
00230
00231 return true;
00232
00233 }
00234
00235 static void GR_1D_FS (int PTS, int N_SYM, int N_UN_SYM,
00236 const SimpleIntArray& SYM, int P_SYM, const SimpleIntArray& UN_SYM,
00237 Real* X, Real* Y)
00238 {
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248 REPORT
00249
00250 Real T;
00251 int JJ,KK,P_UN_SYM;
00252
00253
00254
00255
00256 if (N_SYM > 0)
00257 {
00258 REPORT
00259 SimpleIntArray U(N_SYM);
00260 for(MultiRadixCounter MRC(N_SYM, SYM, U); !MRC.Finish(); ++MRC)
00261 {
00262 if (MRC.Swap())
00263 {
00264 int P = MRC.Reverse(); int JJ = MRC.Counter(); Real T;
00265 T=X[JJ]; X[JJ]=X[P]; X[P]=T; T=Y[JJ]; Y[JJ]=Y[P]; Y[P]=T;
00266 }
00267 }
00268 }
00269
00270 int J,JL,K,L,M,MS;
00271
00272
00273
00274
00275
00276
00277
00278 if (N_UN_SYM==0) { REPORT return; }
00279 P_UN_SYM=PTS/square(P_SYM); JL=(P_UN_SYM-3)*P_SYM; MS=P_UN_SYM*P_SYM;
00280
00281 for (J = P_SYM; J<=JL; J+=P_SYM)
00282 {
00283 K=J;
00284 do K = P_SYM * BitReverse(K / P_SYM, P_UN_SYM, N_UN_SYM, UN_SYM);
00285 while (K<J);
00286
00287 if (K!=J)
00288 {
00289 REPORT
00290 for (L=0; L<P_SYM; L++) for (M=L; M<PTS; M+=MS)
00291 {
00292 JJ=M+J; KK=M+K;
00293 T=X[JJ]; X[JJ]=X[KK]; X[KK]=T; T=Y[JJ]; Y[JJ]=Y[KK]; Y[KK]=T;
00294 }
00295 }
00296 }
00297
00298 return;
00299 }
00300
00301 static void GR_1D_FT (int N, int N_FACTOR, const SimpleIntArray& FACTOR,
00302 Real* X, Real* Y)
00303 {
00304
00305
00306 REPORT
00307
00308 int M = N;
00309
00310 for (int i = 0; i < N_FACTOR; i++)
00311 {
00312 int P = FACTOR[i]; M /= P;
00313
00314 switch(P)
00315 {
00316 case 1: REPORT break;
00317 case 2: REPORT R_2_FTK (N,M,X,Y,X+M,Y+M); break;
00318 case 3: REPORT R_3_FTK (N,M,X,Y,X+M,Y+M,X+2*M,Y+2*M); break;
00319 case 4: REPORT R_4_FTK (N,M,X,Y,X+M,Y+M,X+2*M,Y+2*M,X+3*M,Y+3*M); break;
00320 case 5:
00321 REPORT
00322 R_5_FTK (N,M,X,Y,X+M,Y+M,X+2*M,Y+2*M,X+3*M,Y+3*M,X+4*M,Y+4*M);
00323 break;
00324 case 8:
00325 REPORT
00326 R_8_FTK (N,M,X,Y,X+M,Y+M,X+2*M,Y+2*M,
00327 X+3*M,Y+3*M,X+4*M,Y+4*M,X+5*M,Y+5*M,
00328 X+6*M,Y+6*M,X+7*M,Y+7*M);
00329 break;
00330 case 16:
00331 REPORT
00332 R_16_FTK (N,M,X,Y,X+M,Y+M,X+2*M,Y+2*M,
00333 X+3*M,Y+3*M,X+4*M,Y+4*M,X+5*M,Y+5*M,
00334 X+6*M,Y+6*M,X+7*M,Y+7*M,X+8*M,Y+8*M,
00335 X+9*M,Y+9*M,X+10*M,Y+10*M,X+11*M,Y+11*M,
00336 X+12*M,Y+12*M,X+13*M,Y+13*M,X+14*M,Y+14*M,
00337 X+15*M,Y+15*M);
00338 break;
00339 default: REPORT R_P_FTK (N,M,P,X,Y); break;
00340 }
00341 }
00342
00343 }
00344
00345 static void R_P_FTK (int N, int M, int P, Real* X, Real* Y)
00346
00347
00348 {
00349 REPORT
00350 bool NO_FOLD,ZERO;
00351 Real ANGLE,IS,IU,RS,RU,T,TWOPI,XT,YT;
00352 int J,JJ,K0,K,M_OVER_2,MP,PM,PP,U,V;
00353
00354 Real AA [9][9], BB [9][9];
00355 Real A [18], B [18], C [18], S [18];
00356 Real IA [9], IB [9], RA [9], RB [9];
00357
00358 TWOPI=8.0*atan(1.0);
00359 M_OVER_2=M/2+1; MP=M*P; PP=P/2; PM=P-1;
00360
00361 for (U=0; U<PP; U++)
00362 {
00363 ANGLE=TWOPI*Real(U+1)/Real(P);
00364 JJ=P-U-2;
00365 A[U]=cos(ANGLE); B[U]=sin(ANGLE);
00366 A[JJ]=A[U]; B[JJ]= -B[U];
00367 }
00368
00369 for (U=1; U<=PP; U++)
00370 {
00371 for (V=1; V<=PP; V++)
00372 { JJ=U*V-U*V/P*P; AA[V-1][U-1]=A[JJ-1]; BB[V-1][U-1]=B[JJ-1]; }
00373 }
00374
00375 for (J=0; J<M_OVER_2; J++)
00376 {
00377 NO_FOLD = (J==0 || 2*J==M);
00378 K0=J;
00379 ANGLE=TWOPI*Real(J)/Real(MP); ZERO=ANGLE==0.0;
00380 C[0]=cos(ANGLE); S[0]=sin(ANGLE);
00381 for (U=1; U<PM; U++)
00382 {
00383 C[U]=C[U-1]*C[0]-S[U-1]*S[0];
00384 S[U]=S[U-1]*C[0]+C[U-1]*S[0];
00385 }
00386 goto L700;
00387 L500:
00388 REPORT
00389 if (NO_FOLD) { REPORT goto L1500; }
00390 REPORT
00391 NO_FOLD=true; K0=M-J;
00392 for (U=0; U<PM; U++)
00393 { T=C[U]*A[U]+S[U]*B[U]; S[U]= -S[U]*A[U]+C[U]*B[U]; C[U]=T; }
00394 L700:
00395 REPORT
00396 for (K=K0; K<N; K+=MP)
00397 {
00398 XT=X[K]; YT=Y[K];
00399 for (U=1; U<=PP; U++)
00400 {
00401 RA[U-1]=XT; IA[U-1]=YT;
00402 RB[U-1]=0.0; IB[U-1]=0.0;
00403 }
00404 for (U=1; U<=PP; U++)
00405 {
00406 JJ=P-U;
00407 RS=X[K+M*U]+X[K+M*JJ]; IS=Y[K+M*U]+Y[K+M*JJ];
00408 RU=X[K+M*U]-X[K+M*JJ]; IU=Y[K+M*U]-Y[K+M*JJ];
00409 XT=XT+RS; YT=YT+IS;
00410 for (V=0; V<PP; V++)
00411 {
00412 RA[V]=RA[V]+RS*AA[V][U-1]; IA[V]=IA[V]+IS*AA[V][U-1];
00413 RB[V]=RB[V]+RU*BB[V][U-1]; IB[V]=IB[V]+IU*BB[V][U-1];
00414 }
00415 }
00416 X[K]=XT; Y[K]=YT;
00417 for (U=1; U<=PP; U++)
00418 {
00419 if (!ZERO)
00420 {
00421 REPORT
00422 XT=RA[U-1]+IB[U-1]; YT=IA[U-1]-RB[U-1];
00423 X[K+M*U]=XT*C[U-1]+YT*S[U-1]; Y[K+M*U]=YT*C[U-1]-XT*S[U-1];
00424 JJ=P-U;
00425 XT=RA[U-1]-IB[U-1]; YT=IA[U-1]+RB[U-1];
00426 X[K+M*JJ]=XT*C[JJ-1]+YT*S[JJ-1];
00427 Y[K+M*JJ]=YT*C[JJ-1]-XT*S[JJ-1];
00428 }
00429 else
00430 {
00431 REPORT
00432 X[K+M*U]=RA[U-1]+IB[U-1]; Y[K+M*U]=IA[U-1]-RB[U-1];
00433 JJ=P-U;
00434 X[K+M*JJ]=RA[U-1]-IB[U-1]; Y[K+M*JJ]=IA[U-1]+RB[U-1];
00435 }
00436 }
00437 }
00438 goto L500;
00439 L1500: ;
00440 }
00441 return;
00442 }
00443
00444 static void R_2_FTK (int N, int M, Real* X0, Real* Y0, Real* X1, Real* Y1)
00445
00446 {
00447 REPORT
00448 bool NO_FOLD,ZERO;
00449 int J,K,K0,M2,M_OVER_2;
00450 Real ANGLE,C,IS,IU,RS,RU,S,TWOPI;
00451
00452 M2=M*2; M_OVER_2=M/2+1;
00453 TWOPI=8.0*atan(1.0);
00454
00455 for (J=0; J<M_OVER_2; J++)
00456 {
00457 NO_FOLD = (J==0 || 2*J==M);
00458 K0=J;
00459 ANGLE=TWOPI*Real(J)/Real(M2); ZERO=ANGLE==0.0;
00460 C=cos(ANGLE); S=sin(ANGLE);
00461 goto L200;
00462 L100:
00463 REPORT
00464 if (NO_FOLD) { REPORT goto L600; }
00465 REPORT
00466 NO_FOLD=true; K0=M-J; C= -C;
00467 L200:
00468 REPORT
00469 for (K=K0; K<N; K+=M2)
00470 {
00471 RS=X0[K]+X1[K]; IS=Y0[K]+Y1[K];
00472 RU=X0[K]-X1[K]; IU=Y0[K]-Y1[K];
00473 X0[K]=RS; Y0[K]=IS;
00474 if (!ZERO) { X1[K]=RU*C+IU*S; Y1[K]=IU*C-RU*S; }
00475 else { X1[K]=RU; Y1[K]=IU; }
00476 }
00477 goto L100;
00478 L600: ;
00479 }
00480
00481 return;
00482 }
00483
00484 static void R_3_FTK (int N, int M, Real* X0, Real* Y0, Real* X1, Real* Y1,
00485 Real* X2, Real* Y2)
00486
00487 {
00488 REPORT
00489 bool NO_FOLD,ZERO;
00490 int J,K,K0,M3,M_OVER_2;
00491 Real ANGLE,A,B,C1,C2,S1,S2,T,TWOPI;
00492 Real I0,I1,I2,IA,IB,IS,R0,R1,R2,RA,RB,RS;
00493
00494 M3=M*3; M_OVER_2=M/2+1; TWOPI=8.0*atan(1.0);
00495 A=cos(TWOPI/3.0); B=sin(TWOPI/3.0);
00496
00497 for (J=0; J<M_OVER_2; J++)
00498 {
00499 NO_FOLD = (J==0 || 2*J==M);
00500 K0=J;
00501 ANGLE=TWOPI*Real(J)/Real(M3); ZERO=ANGLE==0.0;
00502 C1=cos(ANGLE); S1=sin(ANGLE);
00503 C2=C1*C1-S1*S1; S2=S1*C1+C1*S1;
00504 goto L200;
00505 L100:
00506 REPORT
00507 if (NO_FOLD) { REPORT goto L600; }
00508 REPORT
00509 NO_FOLD=true; K0=M-J;
00510 T=C1*A+S1*B; S1=C1*B-S1*A; C1=T;
00511 T=C2*A-S2*B; S2= -C2*B-S2*A; C2=T;
00512 L200:
00513 REPORT
00514 for (K=K0; K<N; K+=M3)
00515 {
00516 R0 = X0[K]; I0 = Y0[K];
00517 RS=X1[K]+X2[K]; IS=Y1[K]+Y2[K];
00518 X0[K]=R0+RS; Y0[K]=I0+IS;
00519 RA=R0+RS*A; IA=I0+IS*A;
00520 RB=(X1[K]-X2[K])*B; IB=(Y1[K]-Y2[K])*B;
00521 if (!ZERO)
00522 {
00523 REPORT
00524 R1=RA+IB; I1=IA-RB; R2=RA-IB; I2=IA+RB;
00525 X1[K]=R1*C1+I1*S1; Y1[K]=I1*C1-R1*S1;
00526 X2[K]=R2*C2+I2*S2; Y2[K]=I2*C2-R2*S2;
00527 }
00528 else { REPORT X1[K]=RA+IB; Y1[K]=IA-RB; X2[K]=RA-IB; Y2[K]=IA+RB; }
00529 }
00530 goto L100;
00531 L600: ;
00532 }
00533
00534 return;
00535 }
00536
00537 static void R_4_FTK (int N, int M,
00538 Real* X0, Real* Y0, Real* X1, Real* Y1,
00539 Real* X2, Real* Y2, Real* X3, Real* Y3)
00540
00541 {
00542 REPORT
00543 bool NO_FOLD,ZERO;
00544 int J,K,K0,M4,M_OVER_2;
00545 Real ANGLE,C1,C2,C3,S1,S2,S3,T,TWOPI;
00546 Real I1,I2,I3,IS0,IS1,IU0,IU1,R1,R2,R3,RS0,RS1,RU0,RU1;
00547
00548 M4=M*4; M_OVER_2=M/2+1;
00549 TWOPI=8.0*atan(1.0);
00550
00551 for (J=0; J<M_OVER_2; J++)
00552 {
00553 NO_FOLD = (J==0 || 2*J==M);
00554 K0=J;
00555 ANGLE=TWOPI*Real(J)/Real(M4); ZERO=ANGLE==0.0;
00556 C1=cos(ANGLE); S1=sin(ANGLE);
00557 C2=C1*C1-S1*S1; S2=S1*C1+C1*S1;
00558 C3=C2*C1-S2*S1; S3=S2*C1+C2*S1;
00559 goto L200;
00560 L100:
00561 REPORT
00562 if (NO_FOLD) { REPORT goto L600; }
00563 REPORT
00564 NO_FOLD=true; K0=M-J;
00565 T=C1; C1=S1; S1=T;
00566 C2= -C2;
00567 T=C3; C3= -S3; S3= -T;
00568 L200:
00569 REPORT
00570 for (K=K0; K<N; K+=M4)
00571 {
00572 RS0=X0[K]+X2[K]; IS0=Y0[K]+Y2[K];
00573 RU0=X0[K]-X2[K]; IU0=Y0[K]-Y2[K];
00574 RS1=X1[K]+X3[K]; IS1=Y1[K]+Y3[K];
00575 RU1=X1[K]-X3[K]; IU1=Y1[K]-Y3[K];
00576 X0[K]=RS0+RS1; Y0[K]=IS0+IS1;
00577 if (!ZERO)
00578 {
00579 REPORT
00580 R1=RU0+IU1; I1=IU0-RU1;
00581 R2=RS0-RS1; I2=IS0-IS1;
00582 R3=RU0-IU1; I3=IU0+RU1;
00583 X2[K]=R1*C1+I1*S1; Y2[K]=I1*C1-R1*S1;
00584 X1[K]=R2*C2+I2*S2; Y1[K]=I2*C2-R2*S2;
00585 X3[K]=R3*C3+I3*S3; Y3[K]=I3*C3-R3*S3;
00586 }
00587 else
00588 {
00589 REPORT
00590 X2[K]=RU0+IU1; Y2[K]=IU0-RU1;
00591 X1[K]=RS0-RS1; Y1[K]=IS0-IS1;
00592 X3[K]=RU0-IU1; Y3[K]=IU0+RU1;
00593 }
00594 }
00595 goto L100;
00596 L600: ;
00597 }
00598
00599 return;
00600 }
00601
00602 static void R_5_FTK (int N, int M,
00603 Real* X0, Real* Y0, Real* X1, Real* Y1, Real* X2, Real* Y2,
00604 Real* X3, Real* Y3, Real* X4, Real* Y4)
00605
00606
00607 {
00608 REPORT
00609 bool NO_FOLD,ZERO;
00610 int J,K,K0,M5,M_OVER_2;
00611 Real ANGLE,A1,A2,B1,B2,C1,C2,C3,C4,S1,S2,S3,S4,T,TWOPI;
00612 Real R0,R1,R2,R3,R4,RA1,RA2,RB1,RB2,RS1,RS2,RU1,RU2;
00613 Real I0,I1,I2,I3,I4,IA1,IA2,IB1,IB2,IS1,IS2,IU1,IU2;
00614
00615 M5=M*5; M_OVER_2=M/2+1;
00616 TWOPI=8.0*atan(1.0);
00617 A1=cos(TWOPI/5.0); B1=sin(TWOPI/5.0);
00618 A2=cos(2.0*TWOPI/5.0); B2=sin(2.0*TWOPI/5.0);
00619
00620 for (J=0; J<M_OVER_2; J++)
00621 {
00622 NO_FOLD = (J==0 || 2*J==M);
00623 K0=J;
00624 ANGLE=TWOPI*Real(J)/Real(M5); ZERO=ANGLE==0.0;
00625 C1=cos(ANGLE); S1=sin(ANGLE);
00626 C2=C1*C1-S1*S1; S2=S1*C1+C1*S1;
00627 C3=C2*C1-S2*S1; S3=S2*C1+C2*S1;
00628 C4=C2*C2-S2*S2; S4=S2*C2+C2*S2;
00629 goto L200;
00630 L100:
00631 REPORT
00632 if (NO_FOLD) { REPORT goto L600; }
00633 REPORT
00634 NO_FOLD=true; K0=M-J;
00635 T=C1*A1+S1*B1; S1=C1*B1-S1*A1; C1=T;
00636 T=C2*A2+S2*B2; S2=C2*B2-S2*A2; C2=T;
00637 T=C3*A2-S3*B2; S3= -C3*B2-S3*A2; C3=T;
00638 T=C4*A1-S4*B1; S4= -C4*B1-S4*A1; C4=T;
00639 L200:
00640 REPORT
00641 for (K=K0; K<N; K+=M5)
00642 {
00643 R0=X0[K]; I0=Y0[K];
00644 RS1=X1[K]+X4[K]; IS1=Y1[K]+Y4[K];
00645 RU1=X1[K]-X4[K]; IU1=Y1[K]-Y4[K];
00646 RS2=X2[K]+X3[K]; IS2=Y2[K]+Y3[K];
00647 RU2=X2[K]-X3[K]; IU2=Y2[K]-Y3[K];
00648 X0[K]=R0+RS1+RS2; Y0[K]=I0+IS1+IS2;
00649 RA1=R0+RS1*A1+RS2*A2; IA1=I0+IS1*A1+IS2*A2;
00650 RA2=R0+RS1*A2+RS2*A1; IA2=I0+IS1*A2+IS2*A1;
00651 RB1=RU1*B1+RU2*B2; IB1=IU1*B1+IU2*B2;
00652 RB2=RU1*B2-RU2*B1; IB2=IU1*B2-IU2*B1;
00653 if (!ZERO)
00654 {
00655 REPORT
00656 R1=RA1+IB1; I1=IA1-RB1;
00657 R2=RA2+IB2; I2=IA2-RB2;
00658 R3=RA2-IB2; I3=IA2+RB2;
00659 R4=RA1-IB1; I4=IA1+RB1;
00660 X1[K]=R1*C1+I1*S1; Y1[K]=I1*C1-R1*S1;
00661 X2[K]=R2*C2+I2*S2; Y2[K]=I2*C2-R2*S2;
00662 X3[K]=R3*C3+I3*S3; Y3[K]=I3*C3-R3*S3;
00663 X4[K]=R4*C4+I4*S4; Y4[K]=I4*C4-R4*S4;
00664 }
00665 else
00666 {
00667 REPORT
00668 X1[K]=RA1+IB1; Y1[K]=IA1-RB1;
00669 X2[K]=RA2+IB2; Y2[K]=IA2-RB2;
00670 X3[K]=RA2-IB2; Y3[K]=IA2+RB2;
00671 X4[K]=RA1-IB1; Y4[K]=IA1+RB1;
00672 }
00673 }
00674 goto L100;
00675 L600: ;
00676 }
00677
00678 return;
00679 }
00680
00681 static void R_8_FTK (int N, int M,
00682 Real* X0, Real* Y0, Real* X1, Real* Y1,
00683 Real* X2, Real* Y2, Real* X3, Real* Y3,
00684 Real* X4, Real* Y4, Real* X5, Real* Y5,
00685 Real* X6, Real* Y6, Real* X7, Real* Y7)
00686
00687 {
00688 REPORT
00689 bool NO_FOLD,ZERO;
00690 int J,K,K0,M8,M_OVER_2;
00691 Real ANGLE,C1,C2,C3,C4,C5,C6,C7,E,S1,S2,S3,S4,S5,S6,S7,T,TWOPI;
00692 Real R1,R2,R3,R4,R5,R6,R7,RS0,RS1,RS2,RS3,RU0,RU1,RU2,RU3;
00693 Real I1,I2,I3,I4,I5,I6,I7,IS0,IS1,IS2,IS3,IU0,IU1,IU2,IU3;
00694 Real RSS0,RSS1,RSU0,RSU1,RUS0,RUS1,RUU0,RUU1;
00695 Real ISS0,ISS1,ISU0,ISU1,IUS0,IUS1,IUU0,IUU1;
00696
00697 M8=M*8; M_OVER_2=M/2+1;
00698 TWOPI=8.0*atan(1.0); E=cos(TWOPI/8.0);
00699
00700 for (J=0;J<M_OVER_2;J++)
00701 {
00702 NO_FOLD= (J==0 || 2*J==M);
00703 K0=J;
00704 ANGLE=TWOPI*Real(J)/Real(M8); ZERO=ANGLE==0.0;
00705 C1=cos(ANGLE); S1=sin(ANGLE);
00706 C2=C1*C1-S1*S1; S2=C1*S1+S1*C1;
00707 C3=C2*C1-S2*S1; S3=S2*C1+C2*S1;
00708 C4=C2*C2-S2*S2; S4=S2*C2+C2*S2;
00709 C5=C4*C1-S4*S1; S5=S4*C1+C4*S1;
00710 C6=C4*C2-S4*S2; S6=S4*C2+C4*S2;
00711 C7=C4*C3-S4*S3; S7=S4*C3+C4*S3;
00712 goto L200;
00713 L100:
00714 REPORT
00715 if (NO_FOLD) { REPORT goto L600; }
00716 REPORT
00717 NO_FOLD=true; K0=M-J;
00718 T=(C1+S1)*E; S1=(C1-S1)*E; C1=T;
00719 T=S2; S2=C2; C2=T;
00720 T=(-C3+S3)*E; S3=(C3+S3)*E; C3=T;
00721 C4= -C4;
00722 T= -(C5+S5)*E; S5=(-C5+S5)*E; C5=T;
00723 T= -S6; S6= -C6; C6=T;
00724 T=(C7-S7)*E; S7= -(C7+S7)*E; C7=T;
00725 L200:
00726 REPORT
00727 for (K=K0; K<N; K+=M8)
00728 {
00729 RS0=X0[K]+X4[K]; IS0=Y0[K]+Y4[K];
00730 RU0=X0[K]-X4[K]; IU0=Y0[K]-Y4[K];
00731 RS1=X1[K]+X5[K]; IS1=Y1[K]+Y5[K];
00732 RU1=X1[K]-X5[K]; IU1=Y1[K]-Y5[K];
00733 RS2=X2[K]+X6[K]; IS2=Y2[K]+Y6[K];
00734 RU2=X2[K]-X6[K]; IU2=Y2[K]-Y6[K];
00735 RS3=X3[K]+X7[K]; IS3=Y3[K]+Y7[K];
00736 RU3=X3[K]-X7[K]; IU3=Y3[K]-Y7[K];
00737 RSS0=RS0+RS2; ISS0=IS0+IS2;
00738 RSU0=RS0-RS2; ISU0=IS0-IS2;
00739 RSS1=RS1+RS3; ISS1=IS1+IS3;
00740 RSU1=RS1-RS3; ISU1=IS1-IS3;
00741 RUS0=RU0-IU2; IUS0=IU0+RU2;
00742 RUU0=RU0+IU2; IUU0=IU0-RU2;
00743 RUS1=RU1-IU3; IUS1=IU1+RU3;
00744 RUU1=RU1+IU3; IUU1=IU1-RU3;
00745 T=(RUS1+IUS1)*E; IUS1=(IUS1-RUS1)*E; RUS1=T;
00746 T=(RUU1+IUU1)*E; IUU1=(IUU1-RUU1)*E; RUU1=T;
00747 X0[K]=RSS0+RSS1; Y0[K]=ISS0+ISS1;
00748 if (!ZERO)
00749 {
00750 REPORT
00751 R1=RUU0+RUU1; I1=IUU0+IUU1;
00752 R2=RSU0+ISU1; I2=ISU0-RSU1;
00753 R3=RUS0+IUS1; I3=IUS0-RUS1;
00754 R4=RSS0-RSS1; I4=ISS0-ISS1;
00755 R5=RUU0-RUU1; I5=IUU0-IUU1;
00756 R6=RSU0-ISU1; I6=ISU0+RSU1;
00757 R7=RUS0-IUS1; I7=IUS0+RUS1;
00758 X4[K]=R1*C1+I1*S1; Y4[K]=I1*C1-R1*S1;
00759 X2[K]=R2*C2+I2*S2; Y2[K]=I2*C2-R2*S2;
00760 X6[K]=R3*C3+I3*S3; Y6[K]=I3*C3-R3*S3;
00761 X1[K]=R4*C4+I4*S4; Y1[K]=I4*C4-R4*S4;
00762 X5[K]=R5*C5+I5*S5; Y5[K]=I5*C5-R5*S5;
00763 X3[K]=R6*C6+I6*S6; Y3[K]=I6*C6-R6*S6;
00764 X7[K]=R7*C7+I7*S7; Y7[K]=I7*C7-R7*S7;
00765 }
00766 else
00767 {
00768 REPORT
00769 X4[K]=RUU0+RUU1; Y4[K]=IUU0+IUU1;
00770 X2[K]=RSU0+ISU1; Y2[K]=ISU0-RSU1;
00771 X6[K]=RUS0+IUS1; Y6[K]=IUS0-RUS1;
00772 X1[K]=RSS0-RSS1; Y1[K]=ISS0-ISS1;
00773 X5[K]=RUU0-RUU1; Y5[K]=IUU0-IUU1;
00774 X3[K]=RSU0-ISU1; Y3[K]=ISU0+RSU1;
00775 X7[K]=RUS0-IUS1; Y7[K]=IUS0+RUS1;
00776 }
00777 }
00778 goto L100;
00779 L600: ;
00780 }
00781
00782 return;
00783 }
00784
00785 static void R_16_FTK (int N, int M,
00786 Real* X0, Real* Y0, Real* X1, Real* Y1,
00787 Real* X2, Real* Y2, Real* X3, Real* Y3,
00788 Real* X4, Real* Y4, Real* X5, Real* Y5,
00789 Real* X6, Real* Y6, Real* X7, Real* Y7,
00790 Real* X8, Real* Y8, Real* X9, Real* Y9,
00791 Real* X10, Real* Y10, Real* X11, Real* Y11,
00792 Real* X12, Real* Y12, Real* X13, Real* Y13,
00793 Real* X14, Real* Y14, Real* X15, Real* Y15)
00794
00795 {
00796 REPORT
00797 bool NO_FOLD,ZERO;
00798 int J,K,K0,M16,M_OVER_2;
00799 Real ANGLE,EI1,ER1,E2,EI3,ER3,EI5,ER5,T,TWOPI;
00800 Real RS0,RS1,RS2,RS3,RS4,RS5,RS6,RS7;
00801 Real IS0,IS1,IS2,IS3,IS4,IS5,IS6,IS7;
00802 Real RU0,RU1,RU2,RU3,RU4,RU5,RU6,RU7;
00803 Real IU0,IU1,IU2,IU3,IU4,IU5,IU6,IU7;
00804 Real RUS0,RUS1,RUS2,RUS3,RUU0,RUU1,RUU2,RUU3;
00805 Real ISS0,ISS1,ISS2,ISS3,ISU0,ISU1,ISU2,ISU3;
00806 Real RSS0,RSS1,RSS2,RSS3,RSU0,RSU1,RSU2,RSU3;
00807 Real IUS0,IUS1,IUS2,IUS3,IUU0,IUU1,IUU2,IUU3;
00808 Real RSSS0,RSSS1,RSSU0,RSSU1,RSUS0,RSUS1,RSUU0,RSUU1;
00809 Real ISSS0,ISSS1,ISSU0,ISSU1,ISUS0,ISUS1,ISUU0,ISUU1;
00810 Real RUSS0,RUSS1,RUSU0,RUSU1,RUUS0,RUUS1,RUUU0,RUUU1;
00811 Real IUSS0,IUSS1,IUSU0,IUSU1,IUUS0,IUUS1,IUUU0,IUUU1;
00812 Real R1,R2,R3,R4,R5,R6,R7,R8,R9,R10,R11,R12,R13,R14,R15;
00813 Real I1,I2,I3,I4,I5,I6,I7,I8,I9,I10,I11,I12,I13,I14,I15;
00814 Real C1,C2,C3,C4,C5,C6,C7,C8,C9,C10,C11,C12,C13,C14,C15;
00815 Real S1,S2,S3,S4,S5,S6,S7,S8,S9,S10,S11,S12,S13,S14,S15;
00816
00817 M16=M*16; M_OVER_2=M/2+1;
00818 TWOPI=8.0*atan(1.0);
00819 ER1=cos(TWOPI/16.0); EI1=sin(TWOPI/16.0);
00820 E2=cos(TWOPI/8.0);
00821 ER3=cos(3.0*TWOPI/16.0); EI3=sin(3.0*TWOPI/16.0);
00822 ER5=cos(5.0*TWOPI/16.0); EI5=sin(5.0*TWOPI/16.0);
00823
00824 for (J=0; J<M_OVER_2; J++)
00825 {
00826 NO_FOLD = (J==0 || 2*J==M);
00827 K0=J;
00828 ANGLE=TWOPI*Real(J)/Real(M16);
00829 ZERO=ANGLE==0.0;
00830 C1=cos(ANGLE); S1=sin(ANGLE);
00831 C2=C1*C1-S1*S1; S2=C1*S1+S1*C1;
00832 C3=C2*C1-S2*S1; S3=S2*C1+C2*S1;
00833 C4=C2*C2-S2*S2; S4=S2*C2+C2*S2;
00834 C5=C4*C1-S4*S1; S5=S4*C1+C4*S1;
00835 C6=C4*C2-S4*S2; S6=S4*C2+C4*S2;
00836 C7=C4*C3-S4*S3; S7=S4*C3+C4*S3;
00837 C8=C4*C4-S4*S4; S8=C4*S4+S4*C4;
00838 C9=C8*C1-S8*S1; S9=S8*C1+C8*S1;
00839 C10=C8*C2-S8*S2; S10=S8*C2+C8*S2;
00840 C11=C8*C3-S8*S3; S11=S8*C3+C8*S3;
00841 C12=C8*C4-S8*S4; S12=S8*C4+C8*S4;
00842 C13=C8*C5-S8*S5; S13=S8*C5+C8*S5;
00843 C14=C8*C6-S8*S6; S14=S8*C6+C8*S6;
00844 C15=C8*C7-S8*S7; S15=S8*C7+C8*S7;
00845 goto L200;
00846 L100:
00847 REPORT
00848 if (NO_FOLD) { REPORT goto L600; }
00849 REPORT
00850 NO_FOLD=true; K0=M-J;
00851 T=C1*ER1+S1*EI1; S1= -S1*ER1+C1*EI1; C1=T;
00852 T=(C2+S2)*E2; S2=(C2-S2)*E2; C2=T;
00853 T=C3*ER3+S3*EI3; S3= -S3*ER3+C3*EI3; C3=T;
00854 T=S4; S4=C4; C4=T;
00855 T=S5*ER1-C5*EI1; S5=C5*ER1+S5*EI1; C5=T;
00856 T=(-C6+S6)*E2; S6=(C6+S6)*E2; C6=T;
00857 T=S7*ER3-C7*EI3; S7=C7*ER3+S7*EI3; C7=T;
00858 C8= -C8;
00859 T= -(C9*ER1+S9*EI1); S9=S9*ER1-C9*EI1; C9=T;
00860 T= -(C10+S10)*E2; S10=(-C10+S10)*E2; C10=T;
00861 T= -(C11*ER3+S11*EI3); S11=S11*ER3-C11*EI3; C11=T;
00862 T= -S12; S12= -C12; C12=T;
00863 T= -S13*ER1+C13*EI1; S13= -(C13*ER1+S13*EI1); C13=T;
00864 T=(C14-S14)*E2; S14= -(C14+S14)*E2; C14=T;
00865 T= -S15*ER3+C15*EI3; S15= -(C15*ER3+S15*EI3); C15=T;
00866 L200:
00867 REPORT
00868 for (K=K0; K<N; K+=M16)
00869 {
00870 RS0=X0[K]+X8[K]; IS0=Y0[K]+Y8[K];
00871 RU0=X0[K]-X8[K]; IU0=Y0[K]-Y8[K];
00872 RS1=X1[K]+X9[K]; IS1=Y1[K]+Y9[K];
00873 RU1=X1[K]-X9[K]; IU1=Y1[K]-Y9[K];
00874 RS2=X2[K]+X10[K]; IS2=Y2[K]+Y10[K];
00875 RU2=X2[K]-X10[K]; IU2=Y2[K]-Y10[K];
00876 RS3=X3[K]+X11[K]; IS3=Y3[K]+Y11[K];
00877 RU3=X3[K]-X11[K]; IU3=Y3[K]-Y11[K];
00878 RS4=X4[K]+X12[K]; IS4=Y4[K]+Y12[K];
00879 RU4=X4[K]-X12[K]; IU4=Y4[K]-Y12[K];
00880 RS5=X5[K]+X13[K]; IS5=Y5[K]+Y13[K];
00881 RU5=X5[K]-X13[K]; IU5=Y5[K]-Y13[K];
00882 RS6=X6[K]+X14[K]; IS6=Y6[K]+Y14[K];
00883 RU6=X6[K]-X14[K]; IU6=Y6[K]-Y14[K];
00884 RS7=X7[K]+X15[K]; IS7=Y7[K]+Y15[K];
00885 RU7=X7[K]-X15[K]; IU7=Y7[K]-Y15[K];
00886 RSS0=RS0+RS4; ISS0=IS0+IS4;
00887 RSS1=RS1+RS5; ISS1=IS1+IS5;
00888 RSS2=RS2+RS6; ISS2=IS2+IS6;
00889 RSS3=RS3+RS7; ISS3=IS3+IS7;
00890 RSU0=RS0-RS4; ISU0=IS0-IS4;
00891 RSU1=RS1-RS5; ISU1=IS1-IS5;
00892 RSU2=RS2-RS6; ISU2=IS2-IS6;
00893 RSU3=RS3-RS7; ISU3=IS3-IS7;
00894 RUS0=RU0-IU4; IUS0=IU0+RU4;
00895 RUS1=RU1-IU5; IUS1=IU1+RU5;
00896 RUS2=RU2-IU6; IUS2=IU2+RU6;
00897 RUS3=RU3-IU7; IUS3=IU3+RU7;
00898 RUU0=RU0+IU4; IUU0=IU0-RU4;
00899 RUU1=RU1+IU5; IUU1=IU1-RU5;
00900 RUU2=RU2+IU6; IUU2=IU2-RU6;
00901 RUU3=RU3+IU7; IUU3=IU3-RU7;
00902 T=(RSU1+ISU1)*E2; ISU1=(ISU1-RSU1)*E2; RSU1=T;
00903 T=(RSU3+ISU3)*E2; ISU3=(ISU3-RSU3)*E2; RSU3=T;
00904 T=RUS1*ER3+IUS1*EI3; IUS1=IUS1*ER3-RUS1*EI3; RUS1=T;
00905 T=(RUS2+IUS2)*E2; IUS2=(IUS2-RUS2)*E2; RUS2=T;
00906 T=RUS3*ER5+IUS3*EI5; IUS3=IUS3*ER5-RUS3*EI5; RUS3=T;
00907 T=RUU1*ER1+IUU1*EI1; IUU1=IUU1*ER1-RUU1*EI1; RUU1=T;
00908 T=(RUU2+IUU2)*E2; IUU2=(IUU2-RUU2)*E2; RUU2=T;
00909 T=RUU3*ER3+IUU3*EI3; IUU3=IUU3*ER3-RUU3*EI3; RUU3=T;
00910 RSSS0=RSS0+RSS2; ISSS0=ISS0+ISS2;
00911 RSSS1=RSS1+RSS3; ISSS1=ISS1+ISS3;
00912 RSSU0=RSS0-RSS2; ISSU0=ISS0-ISS2;
00913 RSSU1=RSS1-RSS3; ISSU1=ISS1-ISS3;
00914 RSUS0=RSU0-ISU2; ISUS0=ISU0+RSU2;
00915 RSUS1=RSU1-ISU3; ISUS1=ISU1+RSU3;
00916 RSUU0=RSU0+ISU2; ISUU0=ISU0-RSU2;
00917 RSUU1=RSU1+ISU3; ISUU1=ISU1-RSU3;
00918 RUSS0=RUS0-IUS2; IUSS0=IUS0+RUS2;
00919 RUSS1=RUS1-IUS3; IUSS1=IUS1+RUS3;
00920 RUSU0=RUS0+IUS2; IUSU0=IUS0-RUS2;
00921 RUSU1=RUS1+IUS3; IUSU1=IUS1-RUS3;
00922 RUUS0=RUU0+RUU2; IUUS0=IUU0+IUU2;
00923 RUUS1=RUU1+RUU3; IUUS1=IUU1+IUU3;
00924 RUUU0=RUU0-RUU2; IUUU0=IUU0-IUU2;
00925 RUUU1=RUU1-RUU3; IUUU1=IUU1-IUU3;
00926 X0[K]=RSSS0+RSSS1; Y0[K]=ISSS0+ISSS1;
00927 if (!ZERO)
00928 {
00929 REPORT
00930 R1=RUUS0+RUUS1; I1=IUUS0+IUUS1;
00931 R2=RSUU0+RSUU1; I2=ISUU0+ISUU1;
00932 R3=RUSU0+RUSU1; I3=IUSU0+IUSU1;
00933 R4=RSSU0+ISSU1; I4=ISSU0-RSSU1;
00934 R5=RUUU0+IUUU1; I5=IUUU0-RUUU1;
00935 R6=RSUS0+ISUS1; I6=ISUS0-RSUS1;
00936 R7=RUSS0+IUSS1; I7=IUSS0-RUSS1;
00937 R8=RSSS0-RSSS1; I8=ISSS0-ISSS1;
00938 R9=RUUS0-RUUS1; I9=IUUS0-IUUS1;
00939 R10=RSUU0-RSUU1; I10=ISUU0-ISUU1;
00940 R11=RUSU0-RUSU1; I11=IUSU0-IUSU1;
00941 R12=RSSU0-ISSU1; I12=ISSU0+RSSU1;
00942 R13=RUUU0-IUUU1; I13=IUUU0+RUUU1;
00943 R14=RSUS0-ISUS1; I14=ISUS0+RSUS1;
00944 R15=RUSS0-IUSS1; I15=IUSS0+RUSS1;
00945 X8[K]=R1*C1+I1*S1; Y8[K]=I1*C1-R1*S1;
00946 X4[K]=R2*C2+I2*S2; Y4[K]=I2*C2-R2*S2;
00947 X12[K]=R3*C3+I3*S3; Y12[K]=I3*C3-R3*S3;
00948 X2[K]=R4*C4+I4*S4; Y2[K]=I4*C4-R4*S4;
00949 X10[K]=R5*C5+I5*S5; Y10[K]=I5*C5-R5*S5;
00950 X6[K]=R6*C6+I6*S6; Y6[K]=I6*C6-R6*S6;
00951 X14[K]=R7*C7+I7*S7; Y14[K]=I7*C7-R7*S7;
00952 X1[K]=R8*C8+I8*S8; Y1[K]=I8*C8-R8*S8;
00953 X9[K]=R9*C9+I9*S9; Y9[K]=I9*C9-R9*S9;
00954 X5[K]=R10*C10+I10*S10; Y5[K]=I10*C10-R10*S10;
00955 X13[K]=R11*C11+I11*S11; Y13[K]=I11*C11-R11*S11;
00956 X3[K]=R12*C12+I12*S12; Y3[K]=I12*C12-R12*S12;
00957 X11[K]=R13*C13+I13*S13; Y11[K]=I13*C13-R13*S13;
00958 X7[K]=R14*C14+I14*S14; Y7[K]=I14*C14-R14*S14;
00959 X15[K]=R15*C15+I15*S15; Y15[K]=I15*C15-R15*S15;
00960 }
00961 else
00962 {
00963 REPORT
00964 X8[K]=RUUS0+RUUS1; Y8[K]=IUUS0+IUUS1;
00965 X4[K]=RSUU0+RSUU1; Y4[K]=ISUU0+ISUU1;
00966 X12[K]=RUSU0+RUSU1; Y12[K]=IUSU0+IUSU1;
00967 X2[K]=RSSU0+ISSU1; Y2[K]=ISSU0-RSSU1;
00968 X10[K]=RUUU0+IUUU1; Y10[K]=IUUU0-RUUU1;
00969 X6[K]=RSUS0+ISUS1; Y6[K]=ISUS0-RSUS1;
00970 X14[K]=RUSS0+IUSS1; Y14[K]=IUSS0-RUSS1;
00971 X1[K]=RSSS0-RSSS1; Y1[K]=ISSS0-ISSS1;
00972 X9[K]=RUUS0-RUUS1; Y9[K]=IUUS0-IUUS1;
00973 X5[K]=RSUU0-RSUU1; Y5[K]=ISUU0-ISUU1;
00974 X13[K]=RUSU0-RUSU1; Y13[K]=IUSU0-IUSU1;
00975 X3[K]=RSSU0-ISSU1; Y3[K]=ISSU0+RSSU1;
00976 X11[K]=RUUU0-IUUU1; Y11[K]=IUUU0+RUUU1;
00977 X7[K]=RSUS0-ISUS1; Y7[K]=ISUS0+RSUS1;
00978 X15[K]=RUSS0-IUSS1; Y15[K]=IUSS0+RUSS1;
00979 }
00980 }
00981 goto L100;
00982 L600: ;
00983 }
00984
00985 return;
00986 }
00987
00988
00989
00990
00991 bool FFT_Controller::CanFactor(int PTS)
00992 {
00993 REPORT
00994 const int NP = 16, NQ = 10, PMAX=19;
00995
00996 if (PTS<=1) { REPORT return true; }
00997
00998 int N = PTS, F = 2, P = 0, Q = 0;
00999
01000 while (N > 1)
01001 {
01002 bool fail = true;
01003 for (int J = F; J <= PMAX; J++)
01004 if (N % J == 0) { fail = false; F=J; break; }
01005 if (fail || P >= NP || Q >= NQ) { REPORT return false; }
01006 N /= F;
01007 if (N % F != 0) Q++; else { N /= F; P++; }
01008 }
01009
01010 return true;
01011
01012 }
01013
01014 bool FFT_Controller::OnlyOldFFT;
01015
01016
01017
01018 MultiRadixCounter::MultiRadixCounter(int nx, const SimpleIntArray& rx,
01019 SimpleIntArray& vx)
01020 : Radix(rx), Value(vx), n(nx), reverse(0),
01021 product(1), counter(0), finish(false)
01022 {
01023 REPORT for (int k = 0; k < n; k++) { Value[k] = 0; product *= Radix[k]; }
01024 }
01025
01026 void MultiRadixCounter::operator++()
01027 {
01028 REPORT
01029 counter++; int p = product;
01030 for (int k = 0; k < n; k++)
01031 {
01032 Value[k]++; int p1 = p / Radix[k]; reverse += p1;
01033 if (Value[k] == Radix[k]) { REPORT Value[k] = 0; reverse -= p; p = p1; }
01034 else { REPORT return; }
01035 }
01036 finish = true;
01037 }
01038
01039
01040 static int BitReverse(int x, int prod, int n, const SimpleIntArray& f)
01041 {
01042
01043
01044
01045
01046
01047 REPORT
01048 const int* d = f.Data() + n; int sum = 0; int q = 1;
01049 while (n--)
01050 {
01051 prod /= *(--d);
01052 int c = x / prod; x-= c * prod;
01053 sum += q * c; q *= *d;
01054 }
01055 return sum;
01056 }
01057
01058
01059 #ifdef use_namespace
01060 }
01061 #endif
01062
01063