1 #include "common.h"
2 #include "encoder.h"
3 #include "fft.h"
4
5 /*****************************************************************************
6 * FFT computes fast fourier transform of BLKSIZE samples of data *
7 * uses decimation-in-frequency algorithm described in "Digital *
8 * Signal Processing" by Oppenheim and Schafer, refer to pages 304 *
9 * (flow graph) and 330-332 (Fortran program in problem 5) *
10 * to get the inverse fft, change line 20 from *
11 * w_imag[L] = -sin(PI/le1); *
12 * to *
13 * w_imag[L] = sin(PI/le1); *
14 * *
15 * required constants: *
16 * #define PI 3.14159265358979 *
17 * #define BLKSIZE 1024 *
18 * #define LOGBLKSIZE 10 *
19 * *
20 *****************************************************************************/
21
fft(float * x_real,float * x_imag,float * energy,float * phi)22 void fft (float *x_real, float *x_imag, float *energy, float *phi)
23 {
24 static int M, MM1;
25 static int init = 0, N, NV2, NM1;
26 static double w_real[LOGBLKSIZE], w_imag[LOGBLKSIZE];
27 int i, j, k, ll;
28 int ip, le, le1;
29 double t_real, t_imag, u_real, u_imag;
30
31 if (init == 0) {
32 memset ((char *) w_real, 0, sizeof (w_real)); /* preset statics to 0 */
33 memset ((char *) w_imag, 0, sizeof (w_imag)); /* preset statics to 0 */
34 M = LOGBLKSIZE;
35 MM1 = LOGBLKSIZE - 1;
36 N = BLKSIZE;
37 NV2 = BLKSIZE >> 1;
38 NM1 = BLKSIZE - 1;
39 for (ll = 0; ll < M; ll++) {
40 le = 1 << (M - ll);
41 le1 = le >> 1;
42 w_real[ll] = cos (PI / le1);
43 w_imag[ll] = -sin (PI / le1);
44 }
45 init++;
46 }
47 for (ll = 0; ll < MM1; ll++) {
48 le = 1 << (M - ll);
49 le1 = le >> 1;
50 u_real = 1;
51 u_imag = 0;
52 for (j = 0; j < le1; j++) {
53 for (i = j; i < N; i += le) {
54 ip = i + le1;
55 t_real = x_real[i] + x_real[ip];
56 t_imag = x_imag[i] + x_imag[ip];
57 x_real[ip] = x_real[i] - x_real[ip];
58 x_imag[ip] = x_imag[i] - x_imag[ip];
59 x_real[i] = t_real;
60 x_imag[i] = t_imag;
61 t_real = x_real[ip];
62 x_real[ip] = x_real[ip] * u_real - x_imag[ip] * u_imag;
63 x_imag[ip] = x_imag[ip] * u_real + t_real * u_imag;
64 }
65 t_real = u_real;
66 u_real = u_real * w_real[ll] - u_imag * w_imag[ll];
67 u_imag = u_imag * w_real[ll] + t_real * w_imag[ll];
68 }
69 }
70 /* special case: ll = M-1; all Wn = 1 */
71 for (i = 0; i < N; i += 2) {
72 ip = i + 1;
73 t_real = x_real[i] + x_real[ip];
74 t_imag = x_imag[i] + x_imag[ip];
75 x_real[ip] = x_real[i] - x_real[ip];
76 x_imag[ip] = x_imag[i] - x_imag[ip];
77 x_real[i] = t_real;
78 x_imag[i] = t_imag;
79 energy[i] = x_real[i] * x_real[i] + x_imag[i] * x_imag[i];
80 if (energy[i] <= 0.0005) {
81 phi[i] = 0;
82 energy[i] = 0.0005;
83 } else
84 phi[i] = atan2 ((double) x_imag[i], (double) x_real[i]);
85 energy[ip] = x_real[ip] * x_real[ip] + x_imag[ip] * x_imag[ip];
86 if (energy[ip] == 0)
87 phi[ip] = 0;
88 else
89 phi[ip] = atan2 ((double) x_imag[ip], (double) x_real[ip]);
90 }
91 /* this section reorders the data to the correct ordering */
92 j = 0;
93 for (i = 0; i < NM1; i++) {
94 if (i < j) {
95 /* use this section only if you need the FFT in complex number form *
96 * (and in the correct ordering) */
97 t_real = x_real[j];
98 t_imag = x_imag[j];
99 x_real[j] = x_real[i];
100 x_imag[j] = x_imag[i];
101 x_real[i] = t_real;
102 x_imag[i] = t_imag;
103 /* reorder the energy and phase, phi */
104 t_real = energy[j];
105 energy[j] = energy[i];
106 energy[i] = t_real;
107 t_real = phi[j];
108 phi[j] = phi[i];
109 phi[i] = t_real;
110 }
111 k = NV2;
112 while (k <= j) {
113 j = j - k;
114 k = k >> 1;
115 }
116 j = j + k;
117 }
118 }
119