43#include "sigpr/EST_fft.h"
47#define PI8 0.392699081698724
48#define RT2 1.4142135623731
49#define IRT2 0.707106781186548
51static void FR2TR(
int,
float*,
float*);
52static void FR4TR(
int,
int,
float*,
float*,
float*,
float*);
53static void FORD1(
int,
float*);
54static void FORD2(
int,
float*);
88 M = fastlog2(real.
n());
93 EST_warning(
"Illegal FFT order %d", real.
n());
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++){
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)
271 float *b =
new float[n];
275 float *b=
invec.memory();
278 if (
n2pow <= 0)
return 0;
286 FR2TR(
in, b, b +
in );
291 for(i = 1; i <=
n4pow; i++) {
294 FR4TR(
in,
nn, b, b +
in, b + 2 *
in, b + 3 *
in);
302 for(i = 3; i < n; i += 2) b[i] = -b[i];
314void FR2TR(
int in,
float *
b0,
float *
b1)
318 for (k = 0; k <
in; k++)
327void FR4TR(
int in,
int nn,
float *
b0,
float *
b1,
float *
b2,
float*
b3) {
332 float c1, c2, c3, s1, s2,
s3;
335 int L[16],
L1,
L2,
L3,
L4,
L5,
L6,
L7,
L8,
L9,
L10,
L11,
L12,
L13,
L14,
L15;
336 int j0,
j1,
j2,
j3,
j4,
j5,
j6,
j7,
j8,
j9,
j10,
j11,
j12,
j13,
j14;
340 for(k = 2; k < 16; k++) {
341 switch (signum(L[k-1] - 2)) {
353 L8=L[8];
L7=L[9];
L6=L[10];
L5=L[11];
L4=L[12];
L3=L[13];
L2=L[14];
457void FORD1(
int m,
float *b) {
458 int j, k = 4,
kl = 2, n = 0x1 <<
m;
461 for(
j = 4;
j <= n;
j += 2) {
476void FORD2(
int m,
float *b)
482 int l[16],
l1,
l2,
l3,
l4,
l5,
l6,
l7,
l8,
l9,
l10,
l11,
l12,
l13,
l14,
l15;
483 int j1,
j2,
j3,
j4,
j5,
j6,
j7,
j8,
j9,
j10,
j11,
j12,
j13,
j14;
487 for(k=2;k<=
m;k++) l[k]=l[k-1]/2;
488 for(k=
m;k<=14;k++) l[k+1]=2;
527 if ((n < 2) || (n % 2 != 0))
return(0);
534 if (n > 1)
return(0);
INLINE int n() const
number of items in vector.
INLINE const T & a_no_check(int n) const
read-only const access operator: without bounds checking