43 #include "sigpr/EST_fft.h"
45 #include "EST_error.h"
47 #define PI8 0.392699081698724
48 #define RT2 1.4142135623731
49 #define IRT2 0.707106781186548
51 static void FR2TR(
int,
float*,
float*);
52 static void FR4TR(
int,
int,
float*,
float*,
float*,
float*);
53 static void FORD1(
int,
float*);
54 static void FORD2(
int,
float*);
83 float tmp_real,tmp_imag;
88 M = fastlog2(real.
n());
89 N = (int)pow(
float(2.0),(float)M);
93 EST_warning(
"Illegal FFT order %d", real.
n());
99 int le = (int)pow(
float(2.0),(float)(M+1-l));
106 w_imag=f * sin(PI/le1);
110 for (i=j;i<=N-le1;i+=le)
119 real.
a_no_check(ip-1) = tmp_real*u_real - tmp_imag*u_imag;
120 imag.
a_no_check(ip-1) = tmp_real*u_imag + tmp_imag*u_real;
126 tmp_real = u_real*w_real - u_imag*w_imag;
127 tmp_imag = u_real*w_imag + u_imag*w_real;
142 for (i=1; i<=NM1;i++)
149 real[j-1] = real(i-1);
150 imag[j-1] = imag(i-1);
175 return slowFFTsub(real,imag,-1.0);
185 if (slowFFTsub(real,imag,1.0) != 0)
188 for(
int i=1;i<=N;i++){
189 real[i-1] /= (float)N;
190 imag[i-1] /= (float)N;
199 if (slowFFT(real, imag) != 0)
203 for(i=0;i<real.
n();i++)
204 real[i] = imag[i] = (real(i)*real(i) + imag(i)*imag(i));
212 if (slowFFT(real,imag) != 0)
216 for(i=0;i<real.
n();i++)
217 real[i] = imag[i] = sqrt((real(i)*real(i) + imag(i)*imag(i)));
225 if (fastFFT(real) == 0)
230 for(i=0,j=0, k=1;i<n;i+=2,j++,k+=2)
254 #define signum(i) (i < 0 ? -1 : i == 0 ? 0 : 1)
260 int i, in, nn, n2pow, n4pow, nthpo;
272 float *b =
new float[n];
279 if (n2pow <= 0)
return 0;
289 FR2TR(in, b, b + in );
294 for(i = 1; i <= n4pow; i++) {
297 FR4TR(in, nn, b, b + in, b + 2 * in, b + 3 * in);
305 for(i = 3; i < n; i += 2) b[i] = -b[i];
317 void FR2TR(
int in,
float *b0,
float *b1)
321 for (k = 0; k < in; k++)
324 b1[k] = b0[k] - b1[k];
330 void FR4TR(
int in,
int nn,
float *b0,
float *b1,
float *b2,
float* b3) {
331 float arg, piovn, th2;
332 float *b4 = b0, *b5 = b1, *b6 = b2, *b7 = b3;
333 float t0, t1, t2, t3, t4, t5, t6, t7;
334 float r1, r5, pr, pi;
335 float c1, c2, c3, s1, s2, s3;
337 int j, k, jj, kk, jthet, jlast, ji, jl, jr, int4;
338 int L[16], L1, L2, L3, L4, L5, L6, L7, L8, L9, L10, L11, L12, L13, L14, L15;
339 int j0, j1, j2, j3, j4, j5, j6, j7, j8, j9, j10, j11, j12, j13, j14;
343 for(k = 2; k < 16; k++) {
344 switch (signum(L[k-1] - 2)) {
355 L15=L[1]; L14=L[2]; L13=L[3]; L12=L[4]; L11=L[5]; L10=L[6]; L9=L[7];
356 L8=L[8]; L7=L[9]; L6=L[10]; L5=L[11]; L4=L[12]; L3=L[13]; L2=L[14];
364 for(j1=2;j1<=L1;j1+=2)
365 for(j2=j1;j2<=L2;j2+=L1)
366 for(j3=j2;j3<=L3;j3+=L2)
367 for(j4=j3;j4<=L4;j4+=L3)
368 for(j5=j4;j5<=L5;j5+=L4)
369 for(j6=j5;j6<=L6;j6+=L5)
370 for(j7=j6;j7<=L7;j7+=L6)
371 for(j8=j7;j8<=L8;j8+=L7)
372 for(j9=j8;j9<=L9;j9+=L8)
373 for(j10=j9;j10<=L10;j10+=L9)
374 for(j11=j10;j11<=L11;j11+=L10)
375 for(j12=j11;j12<=L12;j12+=L11)
376 for(j13=j12;j13<=L13;j13+=L12)
377 for(j14=j13;j14<=L14;j14+=L13)
378 for(jthet=j14;jthet<=L15;jthet+=L14)
387 b2[k] = b0[k] - b2[k];
388 b3[k] = b1[k] - b3[k];
399 pr = IRT2 * (b1[kk]-b3[kk]);
400 pi = IRT2 * (b1[kk]+b3[kk]);
401 b3[kk] = b2[kk] + pi;
402 b1[kk] = pi - b2[kk];
403 b2[kk] = b0[kk] - pr;
404 b0[kk] = b0[kk] + pr;
422 for(j=j0;j<=jlast;j++)
426 r1 = b1[jj]*c1 - b5[kk]*s1;
427 r5 = b1[jj]*s1 + b5[kk]*c1;
428 t2 = b2[jj]*c2 - b6[kk]*s2;
429 t6 = b2[jj]*s2 + b6[kk]*c2;
430 t3 = b3[jj]*c3 - b7[kk]*s3;
431 t7 = b3[jj]*s3 + b7[kk]*c3;
460 void FORD1(
int m,
float *b) {
461 int j, k = 4, kl = 2, n = 0x1 << m;
464 for(j = 4; j <= n; j += 2) {
479 void FORD2(
int m,
float *b)
483 int n = 0x1<<m, k, ij, ji, ij1, ji1;
485 int l[16], l1, l2, l3, l4, l5, l6, l7, l8, l9, l10, l11, l12, l13, l14, l15;
486 int j1, j2, j3, j4, j5, j6, j7, j8, j9, j10, j11, j12, j13, j14;
490 for(k=2;k<=m;k++) l[k]=l[k-1]/2;
491 for(k=m;k<=14;k++) l[k+1]=2;
493 l15=l[1];l14=l[2];l13=l[3];l12=l[4];l11=l[5];l10=l[6];l9=l[7];
494 l8=l[8];l7=l[9];l6=l[10];l5=l[11];l4=l[12];l3=l[13];l2=l[14];l1=l[15];
498 for(j1=2;j1<=l1;j1+=2)
499 for(j2=j1;j2<=l2;j2+=l1)
500 for(j3=j2;j3<=l3;j3+=l2)
501 for(j4=j3;j4<=l4;j4+=l3)
502 for(j5=j4;j5<=l5;j5+=l4)
503 for(j6=j5;j6<=l6;j6+=l5)
504 for(j7=j6;j7<=l7;j7+=l6)
505 for(j8=j7;j8<=l8;j8+=l7)
506 for(j9=j8;j9<=l9;j9+=l8)
507 for(j10=j9;j10<=l10;j10+=l9)
508 for(j11=j10;j11<=l11;j11+=l10)
509 for(j12=j11;j12<=l12;j12+=l11)
510 for(j13=j12;j13<=l13;j13+=l12)
511 for(j14=j13;j14<=l14;j14+=l13)
512 for(ji=j14;ji<=l15;ji+=l14) {
513 ij1 = ij-1; ji1 = ji - 1;
527 int fastlog2(
int n) {
528 int num_bits, power = 0;
530 if ((n < 2) || (n % 2 != 0))
return(0);
531 num_bits =
sizeof(int) * 8;
533 while(power <= num_bits) {
537 if (n > 1)
return(0);