newfft.cpp
Go to the documentation of this file.
1 
6 
7 
8 // This is originally by Sande and Gentleman in 1967! I have translated from
9 // Fortran into C and a little bit of C++.
10 
11 // It takes about twice as long as fftw
12 // (http://theory.lcs.mit.edu/~fftw/homepage.html)
13 // but is much shorter than fftw and so despite its age
14 // might represent a reasonable
15 // compromise between speed and complexity.
16 // If you really need the speed get fftw.
17 
18 
19 // THIS SUBROUTINE WAS WRITTEN BY G.SANDE OF PRINCETON UNIVERSITY AND
20 // W.M.GENTLMAN OF THE BELL TELEPHONE LAB. IT WAS BROUGHT TO LONDON
21 // BY DR. M.D. GODFREY AT THE IMPERIAL COLLEGE AND WAS ADAPTED FOR
22 // BURROUGHS 6700 BY D. R. BRILLINGER AND J. PEMBERTON
23 // IT REPRESENTS THE STATE OF THE ART OF COMPUTING COMPLETE FINITE
24 // DISCRETE FOURIER TRANSFORMS AS OF NOV.1967.
25 // OTHER PROGRAMS REQUIRED.
26 // ONLY THOSE SUBROUTINES INCLUDED HERE.
27 // USAGE.
28 // CALL AR1DFT(N,X,Y)
29 // WHERE N IS THE NUMBER OF POINTS IN THE SEQUENCE .
30 // X - IS A ONE-DIMENSIONAL ARRAY CONTAINING THE REAL
31 // PART OF THE SEQUENCE.
32 // Y - IS A ONE-DIMENSIONAL ARRAY CONTAINING THE
33 // IMAGINARY PART OF THE SEQUENCE.
34 // THE TRANSFORM IS RETURNED IN X AND Y.
35 // METHOD
36 // FOR A GENERAL DISCUSSION OF THESE TRANSFORMS AND OF
37 // THE FAST METHOD FOR COMPUTING THEM, SEE GENTLEMAN AND SANDE,
38 // @FAST FOURIER TRANSFORMS - FOR FUN AND PROFIT,@ 1966 FALL JOINT
39 // COMPUTER CONFERENCE.
40 // THIS PROGRAM COMPUTES THIS FOR A COMPLEX SEQUENCE Z(T) OF LENGTH
41 // N WHOSE ELEMENTS ARE STORED AT(X(I) , Y(I)) AND RETURNS THE
42 // TRANSFORM COEFFICIENTS AT (X(I), Y(I)).
43 // DESCRIPTION
44 // AR1DFT IS A HIGHLY MODULAR ROUTINE CAPABLE OF COMPUTING IN PLACE
45 // THE COMPLETE FINITE DISCRETE FOURIER TRANSFORM OF A ONE-
46 // DIMENSIONAL SEQUENCE OF RATHER GENERAL LENGTH N.
47 // THE MAIN ROUTINE , AR1DFT ITSELF, FACTORS N. IT THEN CALLS ON
48 // ON GR 1D FT TO COMPUTE THE ACTUAL TRANSFORMS, USING THESE FACTORS.
49 // THIS GR 1D FT DOES, CALLING AT EACH STAGE ON THE APPROPRIATE KERN
50 // EL R2FTK, R4FTK, R8FTK, R16FTK, R3FTK, R5FTK, OR RPFTK TO PERFORM
51 // THE COMPUTATIONS FOR THIS PASS OVER THE SEQUENCE, DEPENDING ON
52 // WHETHER THE CORRESPONDING FACTOR IS 2, 4, 8, 16, 3, 5, OR SOME
53 // MORE GENERAL PRIME P. WHEN GR1DFT IS FINISHED THE TRANSFORM IS
54 // COMPUTED, HOWEVER, THE RESULTS ARE STORED IN "DIGITS REVERSED"
55 // ORDER. AR1DFT THEREFORE, CALLS UPON GR 1S FS TO SORT THEM OUT.
56 // TO RETURN TO THE FACTORIZATION, SINGLETON HAS POINTED OUT THAT
57 // THE TRANSFORMS ARE MORE EFFICIENT IF THE SAMPLE SIZE N, IS OF THE
58 // FORM B*A**2 AND B CONSISTS OF A SINGLE FACTOR. IN SUCH A CASE
59 // IF WE PROCESS THE FACTORS IN THE ORDER ABA THEN
60 // THE REORDERING CAN BE DONE AS FAST IN PLACE, AS WITH SCRATCH
61 // STORAGE. BUT AS B BECOMES MORE COMPLICATED, THE COST OF THE DIGIT
62 // REVERSING DUE TO B PART BECOMES VERY EXPENSIVE IF WE TRY TO DO IT
63 // IN PLACE. IN SUCH A CASE IT MIGHT BE BETTER TO USE EXTRA STORAGE
64 // A ROUTINE TO DO THIS IS, HOWEVER, NOT INCLUDED HERE.
65 // ANOTHER FEATURE INFLUENCING THE FACTORIZATION IS THAT FOR ANY FIXED
66 // FACTOR N WE CAN PREPARE A SPECIAL KERNEL WHICH WILL COMPUTE
67 // THAT STAGE OF THE TRANSFORM MORE EFFICIENTLY THAN WOULD A KERNEL
68 // FOR GENERAL FACTORS, ESPECIALLY IF THE GENERAL KERNEL HAD TO BE
69 // APPLIED SEVERAL TIMES. FOR EXAMPLE, FACTORS OF 4 ARE MORE
70 // EFFICIENT THAN FACTORS OF 2, FACTORS OF 8 MORE EFFICIENT THAN 4,ETC
71 // ON THE OTHER HAND DIMINISHING RETURNS RAPIDLY SET IN, ESPECIALLY
72 // SINCE THE LENGTH OF THE KERNEL FOR A SPECIAL CASE IS ROUGHLY
73 // PROPORTIONAL TO THE FACTOR IT DEALS WITH. HENCE THESE PROBABLY ARE
74 // ALL THE KERNELS WE WISH TO HAVE.
75 // RESTRICTIONS.
76 // AN UNFORTUNATE FEATURE OF THE SORTING PROBLEM IS THAT THE MOST
77 // EFFICIENT WAY TO DO IT IS WITH NESTED DO LOOPS, ONE FOR EACH
78 // FACTOR. THIS PUTS A RESTRICTION ON N AS TO HOW MANY FACTORS IT
79 // CAN HAVE. CURRENTLY THE LIMIT IS 16, BUT THE LIMIT CAN BE READILY
80 // RAISED IF NECESSARY.
81 // A SECOND RESTRICTION OF THE PROGRAM IS THAT LOCAL STORAGE OF THE
82 // THE ORDER P**2 IS REQUIRED BY THE GENERAL KERNEL RPFTK, SO SOME
83 // LIMIT MUST BE SET ON P. CURRENTLY THIS IS 19, BUT IT CAN BE INCRE
84 // INCREASED BY TRIVIAL CHANGES.
85 // OTHER COMMENTS.
86 //(1) THE ROUTINE IS ADAPTED TO CHECK WHETHER A GIVEN N WILL MEET THE
87 // ABOVE FACTORING REQUIREMENTS AN, IF NOT, TO RETURN THE NEXT HIGHER
88 // NUMBER, NX, SAY, WHICH WILL MEET THESE REQUIREMENTS.
89 // THIS CAN BE ACCHIEVED BY A STATEMENT OF THE FORM
90 // CALL FACTR(N,X,Y).
91 // IF A DIFFERENT N, SAY NX, IS RETURNED THEN THE TRANSFORMS COULD BE
92 // OBTAINED BY EXTENDING THE SIZE OF THE X-ARRAY AND Y-ARRAY TO NX,
93 // AND SETTING X(I) = Y(I) = 0., FOR I = N+1, NX.
94 //(2) IF THE SEQUENCE Z IS ONLY A REAL SEQUENCE, THEN THE IMAGINARY PART
95 // Y(I)=0., THIS WILL RETURN THE COSINE TRANSFORM OF THE REAL SEQUENCE
96 // IN X, AND THE SINE TRANSFORM IN Y.
97 
98 
99 #define WANT_STREAM
100 
101 #define WANT_MATH
102 
103 #include "newmatap.h"
104 
105 #ifdef use_namespace
106 namespace NEWMAT {
107 #endif
108 
109 #ifdef DO_REPORT
110 #define REPORT { static ExeCounter ExeCount(__LINE__,20); ++ExeCount; }
111 #else
112 #define REPORT {}
113 #endif
114 
115 inline Real square(Real x) { return x*x; }
116 inline int square(int x) { return x*x; }
117 
118 static void GR_1D_FS (int PTS, int N_SYM, int N_UN_SYM,
119  const SimpleIntArray& SYM, int P_SYM, const SimpleIntArray& UN_SYM,
120  Real* X, Real* Y);
121 static void GR_1D_FT (int N, int N_FACTOR, const SimpleIntArray& FACTOR,
122  Real* X, Real* Y);
123 static void R_P_FTK (int N, int M, int P, Real* X, Real* Y);
124 static void R_2_FTK (int N, int M, Real* X0, Real* Y0, Real* X1, Real* Y1);
125 static void R_3_FTK (int N, int M, Real* X0, Real* Y0, Real* X1, Real* Y1,
126  Real* X2, Real* Y2);
127 static void R_4_FTK (int N, int M,
128  Real* X0, Real* Y0, Real* X1, Real* Y1,
129  Real* X2, Real* Y2, Real* X3, Real* Y3);
130 static void R_5_FTK (int N, int M,
131  Real* X0, Real* Y0, Real* X1, Real* Y1, Real* X2, Real* Y2,
132  Real* X3, Real* Y3, Real* X4, Real* Y4);
133 static void R_8_FTK (int N, int M,
134  Real* X0, Real* Y0, Real* X1, Real* Y1,
135  Real* X2, Real* Y2, Real* X3, Real* Y3,
136  Real* X4, Real* Y4, Real* X5, Real* Y5,
137  Real* X6, Real* Y6, Real* X7, Real* Y7);
138 static void R_16_FTK (int N, int M,
139  Real* X0, Real* Y0, Real* X1, Real* Y1,
140  Real* X2, Real* Y2, Real* X3, Real* Y3,
141  Real* X4, Real* Y4, Real* X5, Real* Y5,
142  Real* X6, Real* Y6, Real* X7, Real* Y7,
143  Real* X8, Real* Y8, Real* X9, Real* Y9,
144  Real* X10, Real* Y10, Real* X11, Real* Y11,
145  Real* X12, Real* Y12, Real* X13, Real* Y13,
146  Real* X14, Real* Y14, Real* X15, Real* Y15);
147 static int BitReverse(int x, int prod, int n, const SimpleIntArray& f);
148 
149 
150 bool FFT_Controller::ar_1d_ft (int PTS, Real* X, Real *Y)
151 {
152 // ARBITRARY RADIX ONE DIMENSIONAL FOURIER TRANSFORM
153 
154  REPORT
155 
156  int F,J,N,NF,P,PMAX,P_SYM,P_TWO,Q,R,TWO_GRP;
157 
158  // NP is maximum number of squared factors allows PTS up to 2**32 at least
159  // NQ is number of not-squared factors - increase if we increase PMAX
160  const int NP = 16, NQ = 10;
161  SimpleIntArray PP(NP), QQ(NQ);
162 
163  TWO_GRP=16; PMAX=19;
164 
165  // PMAX is the maximum factor size
166  // TWO_GRP is the maximum power of 2 handled as a single factor
167  // Doesn't take advantage of combining powers of 2 when calculating
168  // number of factors
169 
170  if (PTS<=1) return true;
171  N=PTS; P_SYM=1; F=2; P=0; Q=0;
172 
173  // P counts the number of squared factors
174  // Q counts the number of the rest
175  // R = 0 for no non-squared factors; 1 otherwise
176 
177  // FACTOR holds all the factors - non-squared ones in the middle
178  // - length is 2*P+Q
179  // SYM also holds all the factors but with the non-squared ones
180  // multiplied together - length is 2*P+R
181  // PP holds the values of the squared factors - length is P
182  // QQ holds the values of the rest - length is Q
183 
184  // P_SYM holds the product of the squared factors
185 
186  // find the factors - load into PP and QQ
187  while (N > 1)
188  {
189  bool fail = true;
190  for (J=F; J<=PMAX; J++)
191  if (N % J == 0) { fail = false; F=J; break; }
192  if (fail || P >= NP || Q >= NQ) return false; // can't factor
193  N /= F;
194  if (N % F != 0) QQ[Q++] = F;
195  else { N /= F; PP[P++] = F; P_SYM *= F; }
196  }
197 
198  R = (Q == 0) ? 0 : 1; // R = 0 if no not-squared factors, 1 otherwise
199 
200  NF = 2*P + Q;
201  SimpleIntArray FACTOR(NF + 1), SYM(2*P + R);
202  FACTOR[NF] = 0; // we need this in the "combine powers of 2"
203 
204  // load into SYM and FACTOR
205  for (J=0; J<P; J++)
206  { SYM[J]=FACTOR[J]=PP[P-1-J]; FACTOR[P+Q+J]=SYM[P+R+J]=PP[J]; }
207 
208  if (Q>0)
209  {
210  REPORT
211  for (J=0; J<Q; J++) FACTOR[P+J]=QQ[J];
212  SYM[P]=PTS/square(P_SYM);
213  }
214 
215  // combine powers of 2
216  P_TWO = 1;
217  for (J=0; J < NF; J++)
218  {
219  if (FACTOR[J]!=2) continue;
220  P_TWO=P_TWO*2; FACTOR[J]=1;
221  if (P_TWO<TWO_GRP && FACTOR[J+1]==2) continue;
222  FACTOR[J]=P_TWO; P_TWO=1;
223  }
224 
225  if (P==0) R=0;
226  if (Q<=1) Q=0;
227 
228  // do the analysis
229  GR_1D_FT(PTS,NF,FACTOR,X,Y); // the transform
230  GR_1D_FS(PTS,2*P+R,Q,SYM,P_SYM,QQ,X,Y); // the reshuffling
231 
232  return true;
233 
234 }
235 
236 static void GR_1D_FS (int PTS, int N_SYM, int N_UN_SYM,
237  const SimpleIntArray& SYM, int P_SYM, const SimpleIntArray& UN_SYM,
238  Real* X, Real* Y)
239 {
240 // GENERAL RADIX ONE DIMENSIONAL FOURIER SORT
241 
242 // PTS = number of points
243 // N_SYM = length of SYM
244 // N_UN_SYM = length of UN_SYM
245 // SYM: squared factors + product of non-squared factors + squared factors
246 // P_SYM = product of squared factors (each included only once)
247 // UN_SYM: not-squared factors
248 
249  REPORT
250 
251  Real T;
252  int JJ,KK,P_UN_SYM;
253 
254  // I have replaced the multiple for-loop used by Sande-Gentleman code
255  // by the following code which does not limit the number of factors
256 
257  if (N_SYM > 0)
258  {
259  REPORT
260  SimpleIntArray U(N_SYM);
261  for(MultiRadixCounter MRC(N_SYM, SYM, U); !MRC.Finish(); ++MRC)
262  {
263  if (MRC.Swap())
264  {
265  int P = MRC.Reverse(); int JJ = MRC.Counter(); Real T;
266  T=X[JJ]; X[JJ]=X[P]; X[P]=T; T=Y[JJ]; Y[JJ]=Y[P]; Y[P]=T;
267  }
268  }
269  }
270 
271  int J,JL,K,L,M,MS;
272 
273  // UN_SYM contains the non-squared factors
274  // I have replaced the Sande-Gentleman code as it runs into
275  // integer overflow problems
276  // My code (and theirs) would be improved by using a bit array
277  // as suggested by Van Loan
278 
279  if (N_UN_SYM==0) { REPORT return; }
280  P_UN_SYM=PTS/square(P_SYM); JL=(P_UN_SYM-3)*P_SYM; MS=P_UN_SYM*P_SYM;
281 
282  for (J = P_SYM; J<=JL; J+=P_SYM)
283  {
284  K=J;
285  do K = P_SYM * BitReverse(K / P_SYM, P_UN_SYM, N_UN_SYM, UN_SYM);
286  while (K<J);
287 
288  if (K!=J)
289  {
290  REPORT
291  for (L=0; L<P_SYM; L++) for (M=L; M<PTS; M+=MS)
292  {
293  JJ=M+J; KK=M+K;
294  T=X[JJ]; X[JJ]=X[KK]; X[KK]=T; T=Y[JJ]; Y[JJ]=Y[KK]; Y[KK]=T;
295  }
296  }
297  }
298 
299  return;
300 }
301 
302 static void GR_1D_FT (int N, int N_FACTOR, const SimpleIntArray& FACTOR,
303  Real* X, Real* Y)
304 {
305 // GENERAL RADIX ONE DIMENSIONAL FOURIER TRANSFORM;
306 
307  REPORT
308 
309  int M = N;
310 
311  for (int i = 0; i < N_FACTOR; i++)
312  {
313  int P = FACTOR[i]; M /= P;
314 
315  switch(P)
316  {
317  case 1: REPORT break;
318  case 2: REPORT R_2_FTK (N,M,X,Y,X+M,Y+M); break;
319  case 3: REPORT R_3_FTK (N,M,X,Y,X+M,Y+M,X+2*M,Y+2*M); break;
320  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;
321  case 5:
322  REPORT
323  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);
324  break;
325  case 8:
326  REPORT
327  R_8_FTK (N,M,X,Y,X+M,Y+M,X+2*M,Y+2*M,
328  X+3*M,Y+3*M,X+4*M,Y+4*M,X+5*M,Y+5*M,
329  X+6*M,Y+6*M,X+7*M,Y+7*M);
330  break;
331  case 16:
332  REPORT
333  R_16_FTK (N,M,X,Y,X+M,Y+M,X+2*M,Y+2*M,
334  X+3*M,Y+3*M,X+4*M,Y+4*M,X+5*M,Y+5*M,
335  X+6*M,Y+6*M,X+7*M,Y+7*M,X+8*M,Y+8*M,
336  X+9*M,Y+9*M,X+10*M,Y+10*M,X+11*M,Y+11*M,
337  X+12*M,Y+12*M,X+13*M,Y+13*M,X+14*M,Y+14*M,
338  X+15*M,Y+15*M);
339  break;
340  default: REPORT R_P_FTK (N,M,P,X,Y); break;
341  }
342  }
343 
344 }
345 
346 static void R_P_FTK (int N, int M, int P, Real* X, Real* Y)
347 // RADIX PRIME FOURIER TRANSFORM KERNEL;
348 // X and Y are treated as M * P matrices with Fortran storage
349 {
350  REPORT
351  bool NO_FOLD,ZERO;
352  Real ANGLE,IS,IU,RS,RU,T,TWOPI,XT,YT;
353  int J,JJ,K0,K,M_OVER_2,MP,PM,PP,U,V;
354 
355  Real AA [9][9], BB [9][9];
356  Real A [18], B [18], C [18], S [18];
357  Real IA [9], IB [9], RA [9], RB [9];
358 
359  TWOPI=8.0*atan(1.0);
360  M_OVER_2=M/2+1; MP=M*P; PP=P/2; PM=P-1;
361 
362  for (U=0; U<PP; U++)
363  {
364  ANGLE=TWOPI*Real(U+1)/Real(P);
365  JJ=P-U-2;
366  A[U]=cos(ANGLE); B[U]=sin(ANGLE);
367  A[JJ]=A[U]; B[JJ]= -B[U];
368  }
369 
370  for (U=1; U<=PP; U++)
371  {
372  for (V=1; V<=PP; V++)
373  { JJ=U*V-U*V/P*P; AA[V-1][U-1]=A[JJ-1]; BB[V-1][U-1]=B[JJ-1]; }
374  }
375 
376  for (J=0; J<M_OVER_2; J++)
377  {
378  NO_FOLD = (J==0 || 2*J==M);
379  K0=J;
380  ANGLE=TWOPI*Real(J)/Real(MP); ZERO=ANGLE==0.0;
381  C[0]=cos(ANGLE); S[0]=sin(ANGLE);
382  for (U=1; U<PM; U++)
383  {
384  C[U]=C[U-1]*C[0]-S[U-1]*S[0];
385  S[U]=S[U-1]*C[0]+C[U-1]*S[0];
386  }
387  goto L700;
388  L500:
389  REPORT
390  if (NO_FOLD) { REPORT goto L1500; }
391  REPORT
392  NO_FOLD=true; K0=M-J;
393  for (U=0; U<PM; U++)
394  { T=C[U]*A[U]+S[U]*B[U]; S[U]= -S[U]*A[U]+C[U]*B[U]; C[U]=T; }
395  L700:
396  REPORT
397  for (K=K0; K<N; K+=MP)
398  {
399  XT=X[K]; YT=Y[K];
400  for (U=1; U<=PP; U++)
401  {
402  RA[U-1]=XT; IA[U-1]=YT;
403  RB[U-1]=0.0; IB[U-1]=0.0;
404  }
405  for (U=1; U<=PP; U++)
406  {
407  JJ=P-U;
408  RS=X[K+M*U]+X[K+M*JJ]; IS=Y[K+M*U]+Y[K+M*JJ];
409  RU=X[K+M*U]-X[K+M*JJ]; IU=Y[K+M*U]-Y[K+M*JJ];
410  XT=XT+RS; YT=YT+IS;
411  for (V=0; V<PP; V++)
412  {
413  RA[V]=RA[V]+RS*AA[V][U-1]; IA[V]=IA[V]+IS*AA[V][U-1];
414  RB[V]=RB[V]+RU*BB[V][U-1]; IB[V]=IB[V]+IU*BB[V][U-1];
415  }
416  }
417  X[K]=XT; Y[K]=YT;
418  for (U=1; U<=PP; U++)
419  {
420  if (!ZERO)
421  {
422  REPORT
423  XT=RA[U-1]+IB[U-1]; YT=IA[U-1]-RB[U-1];
424  X[K+M*U]=XT*C[U-1]+YT*S[U-1]; Y[K+M*U]=YT*C[U-1]-XT*S[U-1];
425  JJ=P-U;
426  XT=RA[U-1]-IB[U-1]; YT=IA[U-1]+RB[U-1];
427  X[K+M*JJ]=XT*C[JJ-1]+YT*S[JJ-1];
428  Y[K+M*JJ]=YT*C[JJ-1]-XT*S[JJ-1];
429  }
430  else
431  {
432  REPORT
433  X[K+M*U]=RA[U-1]+IB[U-1]; Y[K+M*U]=IA[U-1]-RB[U-1];
434  JJ=P-U;
435  X[K+M*JJ]=RA[U-1]-IB[U-1]; Y[K+M*JJ]=IA[U-1]+RB[U-1];
436  }
437  }
438  }
439  goto L500;
440 L1500: ;
441  }
442  return;
443 }
444 
445 static void R_2_FTK (int N, int M, Real* X0, Real* Y0, Real* X1, Real* Y1)
446 // RADIX TWO FOURIER TRANSFORM KERNEL;
447 {
448  REPORT
449  bool NO_FOLD,ZERO;
450  int J,K,K0,M2,M_OVER_2;
451  Real ANGLE,C,IS,IU,RS,RU,S,TWOPI;
452 
453  M2=M*2; M_OVER_2=M/2+1;
454  TWOPI=8.0*atan(1.0);
455 
456  for (J=0; J<M_OVER_2; J++)
457  {
458  NO_FOLD = (J==0 || 2*J==M);
459  K0=J;
460  ANGLE=TWOPI*Real(J)/Real(M2); ZERO=ANGLE==0.0;
461  C=cos(ANGLE); S=sin(ANGLE);
462  goto L200;
463  L100:
464  REPORT
465  if (NO_FOLD) { REPORT goto L600; }
466  REPORT
467  NO_FOLD=true; K0=M-J; C= -C;
468  L200:
469  REPORT
470  for (K=K0; K<N; K+=M2)
471  {
472  RS=X0[K]+X1[K]; IS=Y0[K]+Y1[K];
473  RU=X0[K]-X1[K]; IU=Y0[K]-Y1[K];
474  X0[K]=RS; Y0[K]=IS;
475  if (!ZERO) { X1[K]=RU*C+IU*S; Y1[K]=IU*C-RU*S; }
476  else { X1[K]=RU; Y1[K]=IU; }
477  }
478  goto L100;
479  L600: ;
480  }
481 
482  return;
483 }
484 
485 static void R_3_FTK (int N, int M, Real* X0, Real* Y0, Real* X1, Real* Y1,
486  Real* X2, Real* Y2)
487 // RADIX THREE FOURIER TRANSFORM KERNEL
488 {
489  REPORT
490  bool NO_FOLD,ZERO;
491  int J,K,K0,M3,M_OVER_2;
492  Real ANGLE,A,B,C1,C2,S1,S2,T,TWOPI;
493  Real I0,I1,I2,IA,IB,IS,R0,R1,R2,RA,RB,RS;
494 
495  M3=M*3; M_OVER_2=M/2+1; TWOPI=8.0*atan(1.0);
496  A=cos(TWOPI/3.0); B=sin(TWOPI/3.0);
497 
498  for (J=0; J<M_OVER_2; J++)
499  {
500  NO_FOLD = (J==0 || 2*J==M);
501  K0=J;
502  ANGLE=TWOPI*Real(J)/Real(M3); ZERO=ANGLE==0.0;
503  C1=cos(ANGLE); S1=sin(ANGLE);
504  C2=C1*C1-S1*S1; S2=S1*C1+C1*S1;
505  goto L200;
506  L100:
507  REPORT
508  if (NO_FOLD) { REPORT goto L600; }
509  REPORT
510  NO_FOLD=true; K0=M-J;
511  T=C1*A+S1*B; S1=C1*B-S1*A; C1=T;
512  T=C2*A-S2*B; S2= -C2*B-S2*A; C2=T;
513  L200:
514  REPORT
515  for (K=K0; K<N; K+=M3)
516  {
517  R0 = X0[K]; I0 = Y0[K];
518  RS=X1[K]+X2[K]; IS=Y1[K]+Y2[K];
519  X0[K]=R0+RS; Y0[K]=I0+IS;
520  RA=R0+RS*A; IA=I0+IS*A;
521  RB=(X1[K]-X2[K])*B; IB=(Y1[K]-Y2[K])*B;
522  if (!ZERO)
523  {
524  REPORT
525  R1=RA+IB; I1=IA-RB; R2=RA-IB; I2=IA+RB;
526  X1[K]=R1*C1+I1*S1; Y1[K]=I1*C1-R1*S1;
527  X2[K]=R2*C2+I2*S2; Y2[K]=I2*C2-R2*S2;
528  }
529  else { REPORT X1[K]=RA+IB; Y1[K]=IA-RB; X2[K]=RA-IB; Y2[K]=IA+RB; }
530  }
531  goto L100;
532  L600: ;
533  }
534 
535  return;
536 }
537 
538 static void R_4_FTK (int N, int M,
539  Real* X0, Real* Y0, Real* X1, Real* Y1,
540  Real* X2, Real* Y2, Real* X3, Real* Y3)
541 // RADIX FOUR FOURIER TRANSFORM KERNEL
542 {
543  REPORT
544  bool NO_FOLD,ZERO;
545  int J,K,K0,M4,M_OVER_2;
546  Real ANGLE,C1,C2,C3,S1,S2,S3,T,TWOPI;
547  Real I1,I2,I3,IS0,IS1,IU0,IU1,R1,R2,R3,RS0,RS1,RU0,RU1;
548 
549  M4=M*4; M_OVER_2=M/2+1;
550  TWOPI=8.0*atan(1.0);
551 
552  for (J=0; J<M_OVER_2; J++)
553  {
554  NO_FOLD = (J==0 || 2*J==M);
555  K0=J;
556  ANGLE=TWOPI*Real(J)/Real(M4); ZERO=ANGLE==0.0;
557  C1=cos(ANGLE); S1=sin(ANGLE);
558  C2=C1*C1-S1*S1; S2=S1*C1+C1*S1;
559  C3=C2*C1-S2*S1; S3=S2*C1+C2*S1;
560  goto L200;
561  L100:
562  REPORT
563  if (NO_FOLD) { REPORT goto L600; }
564  REPORT
565  NO_FOLD=true; K0=M-J;
566  T=C1; C1=S1; S1=T;
567  C2= -C2;
568  T=C3; C3= -S3; S3= -T;
569  L200:
570  REPORT
571  for (K=K0; K<N; K+=M4)
572  {
573  RS0=X0[K]+X2[K]; IS0=Y0[K]+Y2[K];
574  RU0=X0[K]-X2[K]; IU0=Y0[K]-Y2[K];
575  RS1=X1[K]+X3[K]; IS1=Y1[K]+Y3[K];
576  RU1=X1[K]-X3[K]; IU1=Y1[K]-Y3[K];
577  X0[K]=RS0+RS1; Y0[K]=IS0+IS1;
578  if (!ZERO)
579  {
580  REPORT
581  R1=RU0+IU1; I1=IU0-RU1;
582  R2=RS0-RS1; I2=IS0-IS1;
583  R3=RU0-IU1; I3=IU0+RU1;
584  X2[K]=R1*C1+I1*S1; Y2[K]=I1*C1-R1*S1;
585  X1[K]=R2*C2+I2*S2; Y1[K]=I2*C2-R2*S2;
586  X3[K]=R3*C3+I3*S3; Y3[K]=I3*C3-R3*S3;
587  }
588  else
589  {
590  REPORT
591  X2[K]=RU0+IU1; Y2[K]=IU0-RU1;
592  X1[K]=RS0-RS1; Y1[K]=IS0-IS1;
593  X3[K]=RU0-IU1; Y3[K]=IU0+RU1;
594  }
595  }
596  goto L100;
597  L600: ;
598  }
599 
600  return;
601 }
602 
603 static void R_5_FTK (int N, int M,
604  Real* X0, Real* Y0, Real* X1, Real* Y1, Real* X2, Real* Y2,
605  Real* X3, Real* Y3, Real* X4, Real* Y4)
606 // RADIX FIVE FOURIER TRANSFORM KERNEL
607 
608 {
609  REPORT
610  bool NO_FOLD,ZERO;
611  int J,K,K0,M5,M_OVER_2;
612  Real ANGLE,A1,A2,B1,B2,C1,C2,C3,C4,S1,S2,S3,S4,T,TWOPI;
613  Real R0,R1,R2,R3,R4,RA1,RA2,RB1,RB2,RS1,RS2,RU1,RU2;
614  Real I0,I1,I2,I3,I4,IA1,IA2,IB1,IB2,IS1,IS2,IU1,IU2;
615 
616  M5=M*5; M_OVER_2=M/2+1;
617  TWOPI=8.0*atan(1.0);
618  A1=cos(TWOPI/5.0); B1=sin(TWOPI/5.0);
619  A2=cos(2.0*TWOPI/5.0); B2=sin(2.0*TWOPI/5.0);
620 
621  for (J=0; J<M_OVER_2; J++)
622  {
623  NO_FOLD = (J==0 || 2*J==M);
624  K0=J;
625  ANGLE=TWOPI*Real(J)/Real(M5); ZERO=ANGLE==0.0;
626  C1=cos(ANGLE); S1=sin(ANGLE);
627  C2=C1*C1-S1*S1; S2=S1*C1+C1*S1;
628  C3=C2*C1-S2*S1; S3=S2*C1+C2*S1;
629  C4=C2*C2-S2*S2; S4=S2*C2+C2*S2;
630  goto L200;
631  L100:
632  REPORT
633  if (NO_FOLD) { REPORT goto L600; }
634  REPORT
635  NO_FOLD=true; K0=M-J;
636  T=C1*A1+S1*B1; S1=C1*B1-S1*A1; C1=T;
637  T=C2*A2+S2*B2; S2=C2*B2-S2*A2; C2=T;
638  T=C3*A2-S3*B2; S3= -C3*B2-S3*A2; C3=T;
639  T=C4*A1-S4*B1; S4= -C4*B1-S4*A1; C4=T;
640  L200:
641  REPORT
642  for (K=K0; K<N; K+=M5)
643  {
644  R0=X0[K]; I0=Y0[K];
645  RS1=X1[K]+X4[K]; IS1=Y1[K]+Y4[K];
646  RU1=X1[K]-X4[K]; IU1=Y1[K]-Y4[K];
647  RS2=X2[K]+X3[K]; IS2=Y2[K]+Y3[K];
648  RU2=X2[K]-X3[K]; IU2=Y2[K]-Y3[K];
649  X0[K]=R0+RS1+RS2; Y0[K]=I0+IS1+IS2;
650  RA1=R0+RS1*A1+RS2*A2; IA1=I0+IS1*A1+IS2*A2;
651  RA2=R0+RS1*A2+RS2*A1; IA2=I0+IS1*A2+IS2*A1;
652  RB1=RU1*B1+RU2*B2; IB1=IU1*B1+IU2*B2;
653  RB2=RU1*B2-RU2*B1; IB2=IU1*B2-IU2*B1;
654  if (!ZERO)
655  {
656  REPORT
657  R1=RA1+IB1; I1=IA1-RB1;
658  R2=RA2+IB2; I2=IA2-RB2;
659  R3=RA2-IB2; I3=IA2+RB2;
660  R4=RA1-IB1; I4=IA1+RB1;
661  X1[K]=R1*C1+I1*S1; Y1[K]=I1*C1-R1*S1;
662  X2[K]=R2*C2+I2*S2; Y2[K]=I2*C2-R2*S2;
663  X3[K]=R3*C3+I3*S3; Y3[K]=I3*C3-R3*S3;
664  X4[K]=R4*C4+I4*S4; Y4[K]=I4*C4-R4*S4;
665  }
666  else
667  {
668  REPORT
669  X1[K]=RA1+IB1; Y1[K]=IA1-RB1;
670  X2[K]=RA2+IB2; Y2[K]=IA2-RB2;
671  X3[K]=RA2-IB2; Y3[K]=IA2+RB2;
672  X4[K]=RA1-IB1; Y4[K]=IA1+RB1;
673  }
674  }
675  goto L100;
676  L600: ;
677  }
678 
679  return;
680 }
681 
682 static void R_8_FTK (int N, int M,
683  Real* X0, Real* Y0, Real* X1, Real* Y1,
684  Real* X2, Real* Y2, Real* X3, Real* Y3,
685  Real* X4, Real* Y4, Real* X5, Real* Y5,
686  Real* X6, Real* Y6, Real* X7, Real* Y7)
687 // RADIX EIGHT FOURIER TRANSFORM KERNEL
688 {
689  REPORT
690  bool NO_FOLD,ZERO;
691  int J,K,K0,M8,M_OVER_2;
692  Real ANGLE,C1,C2,C3,C4,C5,C6,C7,E,S1,S2,S3,S4,S5,S6,S7,T,TWOPI;
693  Real R1,R2,R3,R4,R5,R6,R7,RS0,RS1,RS2,RS3,RU0,RU1,RU2,RU3;
694  Real I1,I2,I3,I4,I5,I6,I7,IS0,IS1,IS2,IS3,IU0,IU1,IU2,IU3;
695  Real RSS0,RSS1,RSU0,RSU1,RUS0,RUS1,RUU0,RUU1;
696  Real ISS0,ISS1,ISU0,ISU1,IUS0,IUS1,IUU0,IUU1;
697 
698  M8=M*8; M_OVER_2=M/2+1;
699  TWOPI=8.0*atan(1.0); E=cos(TWOPI/8.0);
700 
701  for (J=0;J<M_OVER_2;J++)
702  {
703  NO_FOLD= (J==0 || 2*J==M);
704  K0=J;
705  ANGLE=TWOPI*Real(J)/Real(M8); ZERO=ANGLE==0.0;
706  C1=cos(ANGLE); S1=sin(ANGLE);
707  C2=C1*C1-S1*S1; S2=C1*S1+S1*C1;
708  C3=C2*C1-S2*S1; S3=S2*C1+C2*S1;
709  C4=C2*C2-S2*S2; S4=S2*C2+C2*S2;
710  C5=C4*C1-S4*S1; S5=S4*C1+C4*S1;
711  C6=C4*C2-S4*S2; S6=S4*C2+C4*S2;
712  C7=C4*C3-S4*S3; S7=S4*C3+C4*S3;
713  goto L200;
714  L100:
715  REPORT
716  if (NO_FOLD) { REPORT goto L600; }
717  REPORT
718  NO_FOLD=true; K0=M-J;
719  T=(C1+S1)*E; S1=(C1-S1)*E; C1=T;
720  T=S2; S2=C2; C2=T;
721  T=(-C3+S3)*E; S3=(C3+S3)*E; C3=T;
722  C4= -C4;
723  T= -(C5+S5)*E; S5=(-C5+S5)*E; C5=T;
724  T= -S6; S6= -C6; C6=T;
725  T=(C7-S7)*E; S7= -(C7+S7)*E; C7=T;
726  L200:
727  REPORT
728  for (K=K0; K<N; K+=M8)
729  {
730  RS0=X0[K]+X4[K]; IS0=Y0[K]+Y4[K];
731  RU0=X0[K]-X4[K]; IU0=Y0[K]-Y4[K];
732  RS1=X1[K]+X5[K]; IS1=Y1[K]+Y5[K];
733  RU1=X1[K]-X5[K]; IU1=Y1[K]-Y5[K];
734  RS2=X2[K]+X6[K]; IS2=Y2[K]+Y6[K];
735  RU2=X2[K]-X6[K]; IU2=Y2[K]-Y6[K];
736  RS3=X3[K]+X7[K]; IS3=Y3[K]+Y7[K];
737  RU3=X3[K]-X7[K]; IU3=Y3[K]-Y7[K];
738  RSS0=RS0+RS2; ISS0=IS0+IS2;
739  RSU0=RS0-RS2; ISU0=IS0-IS2;
740  RSS1=RS1+RS3; ISS1=IS1+IS3;
741  RSU1=RS1-RS3; ISU1=IS1-IS3;
742  RUS0=RU0-IU2; IUS0=IU0+RU2;
743  RUU0=RU0+IU2; IUU0=IU0-RU2;
744  RUS1=RU1-IU3; IUS1=IU1+RU3;
745  RUU1=RU1+IU3; IUU1=IU1-RU3;
746  T=(RUS1+IUS1)*E; IUS1=(IUS1-RUS1)*E; RUS1=T;
747  T=(RUU1+IUU1)*E; IUU1=(IUU1-RUU1)*E; RUU1=T;
748  X0[K]=RSS0+RSS1; Y0[K]=ISS0+ISS1;
749  if (!ZERO)
750  {
751  REPORT
752  R1=RUU0+RUU1; I1=IUU0+IUU1;
753  R2=RSU0+ISU1; I2=ISU0-RSU1;
754  R3=RUS0+IUS1; I3=IUS0-RUS1;
755  R4=RSS0-RSS1; I4=ISS0-ISS1;
756  R5=RUU0-RUU1; I5=IUU0-IUU1;
757  R6=RSU0-ISU1; I6=ISU0+RSU1;
758  R7=RUS0-IUS1; I7=IUS0+RUS1;
759  X4[K]=R1*C1+I1*S1; Y4[K]=I1*C1-R1*S1;
760  X2[K]=R2*C2+I2*S2; Y2[K]=I2*C2-R2*S2;
761  X6[K]=R3*C3+I3*S3; Y6[K]=I3*C3-R3*S3;
762  X1[K]=R4*C4+I4*S4; Y1[K]=I4*C4-R4*S4;
763  X5[K]=R5*C5+I5*S5; Y5[K]=I5*C5-R5*S5;
764  X3[K]=R6*C6+I6*S6; Y3[K]=I6*C6-R6*S6;
765  X7[K]=R7*C7+I7*S7; Y7[K]=I7*C7-R7*S7;
766  }
767  else
768  {
769  REPORT
770  X4[K]=RUU0+RUU1; Y4[K]=IUU0+IUU1;
771  X2[K]=RSU0+ISU1; Y2[K]=ISU0-RSU1;
772  X6[K]=RUS0+IUS1; Y6[K]=IUS0-RUS1;
773  X1[K]=RSS0-RSS1; Y1[K]=ISS0-ISS1;
774  X5[K]=RUU0-RUU1; Y5[K]=IUU0-IUU1;
775  X3[K]=RSU0-ISU1; Y3[K]=ISU0+RSU1;
776  X7[K]=RUS0-IUS1; Y7[K]=IUS0+RUS1;
777  }
778  }
779  goto L100;
780  L600: ;
781  }
782 
783  return;
784 }
785 
786 static void R_16_FTK (int N, int M,
787  Real* X0, Real* Y0, Real* X1, Real* Y1,
788  Real* X2, Real* Y2, Real* X3, Real* Y3,
789  Real* X4, Real* Y4, Real* X5, Real* Y5,
790  Real* X6, Real* Y6, Real* X7, Real* Y7,
791  Real* X8, Real* Y8, Real* X9, Real* Y9,
792  Real* X10, Real* Y10, Real* X11, Real* Y11,
793  Real* X12, Real* Y12, Real* X13, Real* Y13,
794  Real* X14, Real* Y14, Real* X15, Real* Y15)
795 // RADIX SIXTEEN FOURIER TRANSFORM KERNEL
796 {
797  REPORT
798  bool NO_FOLD,ZERO;
799  int J,K,K0,M16,M_OVER_2;
800  Real ANGLE,EI1,ER1,E2,EI3,ER3,EI5,ER5,T,TWOPI;
801  Real RS0,RS1,RS2,RS3,RS4,RS5,RS6,RS7;
802  Real IS0,IS1,IS2,IS3,IS4,IS5,IS6,IS7;
803  Real RU0,RU1,RU2,RU3,RU4,RU5,RU6,RU7;
804  Real IU0,IU1,IU2,IU3,IU4,IU5,IU6,IU7;
805  Real RUS0,RUS1,RUS2,RUS3,RUU0,RUU1,RUU2,RUU3;
806  Real ISS0,ISS1,ISS2,ISS3,ISU0,ISU1,ISU2,ISU3;
807  Real RSS0,RSS1,RSS2,RSS3,RSU0,RSU1,RSU2,RSU3;
808  Real IUS0,IUS1,IUS2,IUS3,IUU0,IUU1,IUU2,IUU3;
809  Real RSSS0,RSSS1,RSSU0,RSSU1,RSUS0,RSUS1,RSUU0,RSUU1;
810  Real ISSS0,ISSS1,ISSU0,ISSU1,ISUS0,ISUS1,ISUU0,ISUU1;
811  Real RUSS0,RUSS1,RUSU0,RUSU1,RUUS0,RUUS1,RUUU0,RUUU1;
812  Real IUSS0,IUSS1,IUSU0,IUSU1,IUUS0,IUUS1,IUUU0,IUUU1;
813  Real R1,R2,R3,R4,R5,R6,R7,R8,R9,R10,R11,R12,R13,R14,R15;
814  Real I1,I2,I3,I4,I5,I6,I7,I8,I9,I10,I11,I12,I13,I14,I15;
815  Real C1,C2,C3,C4,C5,C6,C7,C8,C9,C10,C11,C12,C13,C14,C15;
816  Real S1,S2,S3,S4,S5,S6,S7,S8,S9,S10,S11,S12,S13,S14,S15;
817 
818  M16=M*16; M_OVER_2=M/2+1;
819  TWOPI=8.0*atan(1.0);
820  ER1=cos(TWOPI/16.0); EI1=sin(TWOPI/16.0);
821  E2=cos(TWOPI/8.0);
822  ER3=cos(3.0*TWOPI/16.0); EI3=sin(3.0*TWOPI/16.0);
823  ER5=cos(5.0*TWOPI/16.0); EI5=sin(5.0*TWOPI/16.0);
824 
825  for (J=0; J<M_OVER_2; J++)
826  {
827  NO_FOLD = (J==0 || 2*J==M);
828  K0=J;
829  ANGLE=TWOPI*Real(J)/Real(M16);
830  ZERO=ANGLE==0.0;
831  C1=cos(ANGLE); S1=sin(ANGLE);
832  C2=C1*C1-S1*S1; S2=C1*S1+S1*C1;
833  C3=C2*C1-S2*S1; S3=S2*C1+C2*S1;
834  C4=C2*C2-S2*S2; S4=S2*C2+C2*S2;
835  C5=C4*C1-S4*S1; S5=S4*C1+C4*S1;
836  C6=C4*C2-S4*S2; S6=S4*C2+C4*S2;
837  C7=C4*C3-S4*S3; S7=S4*C3+C4*S3;
838  C8=C4*C4-S4*S4; S8=C4*S4+S4*C4;
839  C9=C8*C1-S8*S1; S9=S8*C1+C8*S1;
840  C10=C8*C2-S8*S2; S10=S8*C2+C8*S2;
841  C11=C8*C3-S8*S3; S11=S8*C3+C8*S3;
842  C12=C8*C4-S8*S4; S12=S8*C4+C8*S4;
843  C13=C8*C5-S8*S5; S13=S8*C5+C8*S5;
844  C14=C8*C6-S8*S6; S14=S8*C6+C8*S6;
845  C15=C8*C7-S8*S7; S15=S8*C7+C8*S7;
846  goto L200;
847  L100:
848  REPORT
849  if (NO_FOLD) { REPORT goto L600; }
850  REPORT
851  NO_FOLD=true; K0=M-J;
852  T=C1*ER1+S1*EI1; S1= -S1*ER1+C1*EI1; C1=T;
853  T=(C2+S2)*E2; S2=(C2-S2)*E2; C2=T;
854  T=C3*ER3+S3*EI3; S3= -S3*ER3+C3*EI3; C3=T;
855  T=S4; S4=C4; C4=T;
856  T=S5*ER1-C5*EI1; S5=C5*ER1+S5*EI1; C5=T;
857  T=(-C6+S6)*E2; S6=(C6+S6)*E2; C6=T;
858  T=S7*ER3-C7*EI3; S7=C7*ER3+S7*EI3; C7=T;
859  C8= -C8;
860  T= -(C9*ER1+S9*EI1); S9=S9*ER1-C9*EI1; C9=T;
861  T= -(C10+S10)*E2; S10=(-C10+S10)*E2; C10=T;
862  T= -(C11*ER3+S11*EI3); S11=S11*ER3-C11*EI3; C11=T;
863  T= -S12; S12= -C12; C12=T;
864  T= -S13*ER1+C13*EI1; S13= -(C13*ER1+S13*EI1); C13=T;
865  T=(C14-S14)*E2; S14= -(C14+S14)*E2; C14=T;
866  T= -S15*ER3+C15*EI3; S15= -(C15*ER3+S15*EI3); C15=T;
867  L200:
868  REPORT
869  for (K=K0; K<N; K+=M16)
870  {
871  RS0=X0[K]+X8[K]; IS0=Y0[K]+Y8[K];
872  RU0=X0[K]-X8[K]; IU0=Y0[K]-Y8[K];
873  RS1=X1[K]+X9[K]; IS1=Y1[K]+Y9[K];
874  RU1=X1[K]-X9[K]; IU1=Y1[K]-Y9[K];
875  RS2=X2[K]+X10[K]; IS2=Y2[K]+Y10[K];
876  RU2=X2[K]-X10[K]; IU2=Y2[K]-Y10[K];
877  RS3=X3[K]+X11[K]; IS3=Y3[K]+Y11[K];
878  RU3=X3[K]-X11[K]; IU3=Y3[K]-Y11[K];
879  RS4=X4[K]+X12[K]; IS4=Y4[K]+Y12[K];
880  RU4=X4[K]-X12[K]; IU4=Y4[K]-Y12[K];
881  RS5=X5[K]+X13[K]; IS5=Y5[K]+Y13[K];
882  RU5=X5[K]-X13[K]; IU5=Y5[K]-Y13[K];
883  RS6=X6[K]+X14[K]; IS6=Y6[K]+Y14[K];
884  RU6=X6[K]-X14[K]; IU6=Y6[K]-Y14[K];
885  RS7=X7[K]+X15[K]; IS7=Y7[K]+Y15[K];
886  RU7=X7[K]-X15[K]; IU7=Y7[K]-Y15[K];
887  RSS0=RS0+RS4; ISS0=IS0+IS4;
888  RSS1=RS1+RS5; ISS1=IS1+IS5;
889  RSS2=RS2+RS6; ISS2=IS2+IS6;
890  RSS3=RS3+RS7; ISS3=IS3+IS7;
891  RSU0=RS0-RS4; ISU0=IS0-IS4;
892  RSU1=RS1-RS5; ISU1=IS1-IS5;
893  RSU2=RS2-RS6; ISU2=IS2-IS6;
894  RSU3=RS3-RS7; ISU3=IS3-IS7;
895  RUS0=RU0-IU4; IUS0=IU0+RU4;
896  RUS1=RU1-IU5; IUS1=IU1+RU5;
897  RUS2=RU2-IU6; IUS2=IU2+RU6;
898  RUS3=RU3-IU7; IUS3=IU3+RU7;
899  RUU0=RU0+IU4; IUU0=IU0-RU4;
900  RUU1=RU1+IU5; IUU1=IU1-RU5;
901  RUU2=RU2+IU6; IUU2=IU2-RU6;
902  RUU3=RU3+IU7; IUU3=IU3-RU7;
903  T=(RSU1+ISU1)*E2; ISU1=(ISU1-RSU1)*E2; RSU1=T;
904  T=(RSU3+ISU3)*E2; ISU3=(ISU3-RSU3)*E2; RSU3=T;
905  T=RUS1*ER3+IUS1*EI3; IUS1=IUS1*ER3-RUS1*EI3; RUS1=T;
906  T=(RUS2+IUS2)*E2; IUS2=(IUS2-RUS2)*E2; RUS2=T;
907  T=RUS3*ER5+IUS3*EI5; IUS3=IUS3*ER5-RUS3*EI5; RUS3=T;
908  T=RUU1*ER1+IUU1*EI1; IUU1=IUU1*ER1-RUU1*EI1; RUU1=T;
909  T=(RUU2+IUU2)*E2; IUU2=(IUU2-RUU2)*E2; RUU2=T;
910  T=RUU3*ER3+IUU3*EI3; IUU3=IUU3*ER3-RUU3*EI3; RUU3=T;
911  RSSS0=RSS0+RSS2; ISSS0=ISS0+ISS2;
912  RSSS1=RSS1+RSS3; ISSS1=ISS1+ISS3;
913  RSSU0=RSS0-RSS2; ISSU0=ISS0-ISS2;
914  RSSU1=RSS1-RSS3; ISSU1=ISS1-ISS3;
915  RSUS0=RSU0-ISU2; ISUS0=ISU0+RSU2;
916  RSUS1=RSU1-ISU3; ISUS1=ISU1+RSU3;
917  RSUU0=RSU0+ISU2; ISUU0=ISU0-RSU2;
918  RSUU1=RSU1+ISU3; ISUU1=ISU1-RSU3;
919  RUSS0=RUS0-IUS2; IUSS0=IUS0+RUS2;
920  RUSS1=RUS1-IUS3; IUSS1=IUS1+RUS3;
921  RUSU0=RUS0+IUS2; IUSU0=IUS0-RUS2;
922  RUSU1=RUS1+IUS3; IUSU1=IUS1-RUS3;
923  RUUS0=RUU0+RUU2; IUUS0=IUU0+IUU2;
924  RUUS1=RUU1+RUU3; IUUS1=IUU1+IUU3;
925  RUUU0=RUU0-RUU2; IUUU0=IUU0-IUU2;
926  RUUU1=RUU1-RUU3; IUUU1=IUU1-IUU3;
927  X0[K]=RSSS0+RSSS1; Y0[K]=ISSS0+ISSS1;
928  if (!ZERO)
929  {
930  REPORT
931  R1=RUUS0+RUUS1; I1=IUUS0+IUUS1;
932  R2=RSUU0+RSUU1; I2=ISUU0+ISUU1;
933  R3=RUSU0+RUSU1; I3=IUSU0+IUSU1;
934  R4=RSSU0+ISSU1; I4=ISSU0-RSSU1;
935  R5=RUUU0+IUUU1; I5=IUUU0-RUUU1;
936  R6=RSUS0+ISUS1; I6=ISUS0-RSUS1;
937  R7=RUSS0+IUSS1; I7=IUSS0-RUSS1;
938  R8=RSSS0-RSSS1; I8=ISSS0-ISSS1;
939  R9=RUUS0-RUUS1; I9=IUUS0-IUUS1;
940  R10=RSUU0-RSUU1; I10=ISUU0-ISUU1;
941  R11=RUSU0-RUSU1; I11=IUSU0-IUSU1;
942  R12=RSSU0-ISSU1; I12=ISSU0+RSSU1;
943  R13=RUUU0-IUUU1; I13=IUUU0+RUUU1;
944  R14=RSUS0-ISUS1; I14=ISUS0+RSUS1;
945  R15=RUSS0-IUSS1; I15=IUSS0+RUSS1;
946  X8[K]=R1*C1+I1*S1; Y8[K]=I1*C1-R1*S1;
947  X4[K]=R2*C2+I2*S2; Y4[K]=I2*C2-R2*S2;
948  X12[K]=R3*C3+I3*S3; Y12[K]=I3*C3-R3*S3;
949  X2[K]=R4*C4+I4*S4; Y2[K]=I4*C4-R4*S4;
950  X10[K]=R5*C5+I5*S5; Y10[K]=I5*C5-R5*S5;
951  X6[K]=R6*C6+I6*S6; Y6[K]=I6*C6-R6*S6;
952  X14[K]=R7*C7+I7*S7; Y14[K]=I7*C7-R7*S7;
953  X1[K]=R8*C8+I8*S8; Y1[K]=I8*C8-R8*S8;
954  X9[K]=R9*C9+I9*S9; Y9[K]=I9*C9-R9*S9;
955  X5[K]=R10*C10+I10*S10; Y5[K]=I10*C10-R10*S10;
956  X13[K]=R11*C11+I11*S11; Y13[K]=I11*C11-R11*S11;
957  X3[K]=R12*C12+I12*S12; Y3[K]=I12*C12-R12*S12;
958  X11[K]=R13*C13+I13*S13; Y11[K]=I13*C13-R13*S13;
959  X7[K]=R14*C14+I14*S14; Y7[K]=I14*C14-R14*S14;
960  X15[K]=R15*C15+I15*S15; Y15[K]=I15*C15-R15*S15;
961  }
962  else
963  {
964  REPORT
965  X8[K]=RUUS0+RUUS1; Y8[K]=IUUS0+IUUS1;
966  X4[K]=RSUU0+RSUU1; Y4[K]=ISUU0+ISUU1;
967  X12[K]=RUSU0+RUSU1; Y12[K]=IUSU0+IUSU1;
968  X2[K]=RSSU0+ISSU1; Y2[K]=ISSU0-RSSU1;
969  X10[K]=RUUU0+IUUU1; Y10[K]=IUUU0-RUUU1;
970  X6[K]=RSUS0+ISUS1; Y6[K]=ISUS0-RSUS1;
971  X14[K]=RUSS0+IUSS1; Y14[K]=IUSS0-RUSS1;
972  X1[K]=RSSS0-RSSS1; Y1[K]=ISSS0-ISSS1;
973  X9[K]=RUUS0-RUUS1; Y9[K]=IUUS0-IUUS1;
974  X5[K]=RSUU0-RSUU1; Y5[K]=ISUU0-ISUU1;
975  X13[K]=RUSU0-RUSU1; Y13[K]=IUSU0-IUSU1;
976  X3[K]=RSSU0-ISSU1; Y3[K]=ISSU0+RSSU1;
977  X11[K]=RUUU0-IUUU1; Y11[K]=IUUU0+RUUU1;
978  X7[K]=RSUS0-ISUS1; Y7[K]=ISUS0+RSUS1;
979  X15[K]=RUSS0-IUSS1; Y15[K]=IUSS0+RUSS1;
980  }
981  }
982  goto L100;
983  L600: ;
984  }
985 
986  return;
987 }
988 
989 // can the number of points be factorised sufficiently
990 // for the fft to run
991 
993 {
994  REPORT
995  const int NP = 16, NQ = 10, PMAX=19;
996 
997  if (PTS<=1) { REPORT return true; }
998 
999  int N = PTS, F = 2, P = 0, Q = 0;
1000 
1001  while (N > 1)
1002  {
1003  bool fail = true;
1004  for (int J = F; J <= PMAX; J++)
1005  if (N % J == 0) { fail = false; F=J; break; }
1006  if (fail || P >= NP || Q >= NQ) { REPORT return false; }
1007  N /= F;
1008  if (N % F != 0) Q++; else { N /= F; P++; }
1009  }
1010 
1011  return true; // can factorise
1012 
1013 }
1014 
1015 bool FFT_Controller::OnlyOldFFT; // static variable
1016 
1017 // **************************** multi radix counter **********************
1018 
1020  SimpleIntArray& vx)
1021  : Radix(rx), Value(vx), n(nx), reverse(0),
1022  product(1), counter(0), finish(false)
1023 {
1024  REPORT for (int k = 0; k < n; k++) { Value[k] = 0; product *= Radix[k]; }
1025 }
1026 
1028 {
1029  REPORT
1030  counter++; int p = product;
1031  for (int k = 0; k < n; k++)
1032  {
1033  Value[k]++; int p1 = p / Radix[k]; reverse += p1;
1034  if (Value[k] == Radix[k]) { REPORT Value[k] = 0; reverse -= p; p = p1; }
1035  else { REPORT return; }
1036  }
1037  finish = true;
1038 }
1039 
1040 
1041 static int BitReverse(int x, int prod, int n, const SimpleIntArray& f)
1042 {
1043  // x = c[0]+f[0]*(c[1]+f[1]*(c[2]+...
1044  // return c[n-1]+f[n-1]*(c[n-2]+f[n-2]*(c[n-3]+...
1045  // prod is the product of the f[i]
1046  // n is the number of f[i] (don't assume f has the correct length)
1047 
1048  REPORT
1049  const int* d = f.Data() + n; int sum = 0; int q = 1;
1050  while (n--)
1051  {
1052  prod /= *(--d);
1053  int c = x / prod; x-= c * prod;
1054  sum += q * c; q *= *d;
1055  }
1056  return sum;
1057 }
1058 
1059 
1060 #ifdef use_namespace
1061 }
1062 #endif
1063 
1065 
const SimpleIntArray & Radix
Definition: newmatap.h:197
static void R_3_FTK(int N, int M, Real *X0, Real *Y0, Real *X1, Real *Y1, Real *X2, Real *Y2)
Definition: newfft.cpp:485
Real square(Real x)
Definition: newfft.cpp:115
int * Data()
pointer to the data
Definition: newmat.h:1873
static void R_P_FTK(int N, int M, int P, Real *X, Real *Y)
Definition: newfft.cpp:346
Matrix K
Definition: demo.cpp:228
const int n
Definition: newmatap.h:201
bool Finish() const
Definition: newmatap.h:211
double Real
Definition: include.h:307
static void GR_1D_FS(int PTS, int N_SYM, int N_UN_SYM, const SimpleIntArray &SYM, int P_SYM, const SimpleIntArray &UN_SYM, Real *X, Real *Y)
Definition: newfft.cpp:236
MultiRadixCounter(int nx, const SimpleIntArray &rx, SimpleIntArray &vx)
Definition: newfft.cpp:1019
SimpleIntArray & Value
Definition: newmatap.h:200
static void R_2_FTK(int N, int M, Real *X0, Real *Y0, Real *X1, Real *Y1)
Definition: newfft.cpp:445
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)
Definition: newfft.cpp:786
FloatVector * d
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)
Definition: newfft.cpp:603
#define REPORT
Definition: newfft.cpp:112
static bool ar_1d_ft(int PTS, Real *X, Real *Y)
Definition: newfft.cpp:150
static bool OnlyOldFFT
Definition: newmatap.h:145
static int BitReverse(int x, int prod, int n, const SimpleIntArray &f)
Definition: newfft.cpp:1041
Real sum(const BaseMatrix &B)
Definition: newmat.h:2105
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)
Definition: newfft.cpp:538
static bool CanFactor(int PTS)
Definition: newfft.cpp:992
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)
Definition: newfft.cpp:682
static void GR_1D_FT(int N, int N_FACTOR, const SimpleIntArray &FACTOR, Real *X, Real *Y)
Definition: newfft.cpp:302


kni
Author(s): Martin Günther
autogenerated on Fri Jan 3 2020 04:01:16