removed test and main from all
[swftools.git] / lib / bladeenc / subs.c
1 /*
2                         (c) Copyright 1998-2001 - Tord Jansson
3                         ======================================
4
5                 This file is part of the BladeEnc MP3 Encoder, based on
6                 ISO's reference code for MPEG Layer 3 compression, and might
7                 contain smaller or larger sections that are directly taken
8                 from ISO's reference code.
9
10                 All changes to the ISO reference code herein are either
11                 copyrighted by Tord Jansson (tord.jansson@swipnet.se)
12                 or sublicensed to Tord Jansson by a third party.
13
14         BladeEnc is free software; you can redistribute this file
15         and/or modify it under the terms of the GNU Lesser General Public
16         License as published by the Free Software Foundation; either
17         version 2.1 of the License, or (at your option) any later version.
18
19
20
21         ------------    Changes    ------------
22
23         2000-11-04  Andre Piotrowski
24
25         -       speed up: don«t calculate the unneeded phi values!
26
27         2000-11-22      ap
28
29         -       speed up: slightly improved way to handle the 'special case'
30
31         2000-12-05  ap
32
33         -       speed up: implemented prepacking
34
35         2000-12-11  ap
36
37         -       speed up: avoid data reordering
38
39         2000-12-12  ap
40
41         -       moved fft configuration switches to "l3psy.h"
42
43         2001-01-12  ap
44
45         -       use some explicit type casting to avoid compiler warnings
46 */
47
48 /*****************************************************************************
49  * FFT computes fast fourier transform of BLKSIZE samples of data            *
50  *   based on the decimation-in-frequency algorithm described in "Digital    *
51  *   Signal Processing" by Oppenheim and Schafer, refer to pages 304         *
52  *   (flow graph) and 330-332 (Fortran program in problem 5).                *
53  *                                                                           *
54  *   required constants:                                                     *
55  *         PI          3.14159265358979                                      *
56  *         BLKSIZE     should be 2^^M for a positive integer M               *
57  *                                                                           *
58  *****************************************************************************/
59
60 #include        "common.h"
61 #include        "encoder.h"
62 #include        "l3psy.h"
63
64
65
66
67
68 /* The switches have been moved to "l3psy.h" */
69
70
71
72
73
74 int                                             fInit_fft;
75
76
77
78
79
80 void fft (FLOAT x_real[], FLOAT x_imag[], FLOAT energy[], FLOAT phi[], int N)
81 {
82         int                                             i, j, k, off;
83 #if USED_VALUES_ONLY
84         int                                             min_i, max_i;
85 #endif
86         int                                             ip, le, le1;
87         static double                   w_real[BLKSIZE/2], w_imag[BLKSIZE/2];
88         double                                  t_real, t_imag, u_real, u_imag;
89         int                                             N_ORG;
90
91 #if !REORDER_DATA
92         static  int                             swap_l[BLKSIZE/2+1];
93         static  int                             swap_s[BLKSIZE_s/2+1];
94         int                                             *pSwap, a, b;
95 #endif
96
97         double  t1, t2, t3, t4, t5, t6;
98
99         if (fInit_fft == 0)
100         {
101                 for (i = 0;  i < BLKSIZE/2;  i++)
102                 {
103                         w_real[i] =  cos(PI*i/(BLKSIZE/2));
104 #if NORMAL_FFT
105                         w_imag[i] = -sin(PI*i/(BLKSIZE/2));
106 #else
107                         w_imag[i] =  sin(PI*i/(BLKSIZE/2));
108 #endif
109                 }
110
111 #if !REORDER_DATA
112                 j = 0;
113                 for (i = 0;  i < BLKSIZE/2-1;  i++)
114                 {
115                         swap_l[i] = j;  k = BLKSIZE/4;  while (k <= j) {j -= k;  k >>= 1;}  j += k;
116                 }
117                 swap_l[i] = i;  swap_l[i+1] = i+1;
118
119                 j = 0;
120                 for (i = 0;  i < BLKSIZE_s/2-1;  i++)
121                 {
122                         swap_s[i] = j;  k = BLKSIZE_s/4;  while (k <= j) {j -= k;  k >>= 1;}  j += k;
123                 }
124                 swap_s[i] = i;  swap_s[i+1] = i+1;
125 #endif
126
127                 fInit_fft++;
128         }
129
130
131 #if REAL_SEQUENCE
132         N_ORG = N;
133         N >>= 1;
134 #if !PREPACKED
135         /* packing the sequence to the half length */
136         for (i = 0;  i < N;  i++)
137         {
138                 x_real[i] = x_real[2*i];
139                 x_imag[i] = x_real[2*i+1];
140         }
141 #endif
142 #endif
143
144
145         off = BLKSIZE/N;
146         for (le = N;  le > 1;  le = le1)
147         {
148                 le1 = le >> 1;
149
150                         /* special case: k=0; u_real=1.0; u_imag=0.0; j=0 */
151                         for (i = 0;  i < N;  i += le)
152                         {
153                                 ip = i + le1;
154                                 t_real = x_real[i] - x_real[ip];
155                                 t_imag = x_imag[i] - x_imag[ip];
156                                 x_real[i] += x_real[ip];
157                                 x_imag[i] += x_imag[ip];
158                                 x_real[ip] = (FLOAT) t_real;
159                                 x_imag[ip] = (FLOAT) t_imag;
160                         }
161
162                 k = off;
163                 for (j = 1;  j < le1;  j++)
164                 {
165                         u_real = w_real[k];
166                         u_imag = w_imag[k];
167                         for (i = j;  i < N;  i += le)
168                         {
169                                 ip = i + le1;
170                                 t_real = x_real[i] - x_real[ip];
171                                 t_imag = x_imag[i] - x_imag[ip];
172                                 x_real[i] += x_real[ip];
173                                 x_imag[i] += x_imag[ip];
174                                 x_real[ip] = (FLOAT) (t_real*u_real - t_imag*u_imag);
175                                 x_imag[ip] = (FLOAT) (t_imag*u_real + t_real*u_imag);
176                         }
177                         k += off;
178                 }
179
180                 off += off;
181         }
182
183
184 #if REORDER_DATA
185         /* this section reorders the data to the correct ordering */
186         j = 0;
187         for (i = 0;  i < N-1;  i++)
188         {
189                 if (i < j)
190                 {
191                         t_real = x_real[j];
192                         t_imag = x_imag[j];
193                         x_real[j] = x_real[i];
194                         x_imag[j] = x_imag[i];
195                         x_real[i] = (FLOAT) t_real;
196                         x_imag[i] = (FLOAT) t_imag;
197                 }
198                 k=N/2;
199                 while (k <= j)
200                 {
201                         j -= k;
202                         k >>= 1;
203                 }
204                 j += k;
205         }
206 #else
207         /*
208                 We don«t reorder the data to the correct ordering,
209                 but access the data by the bit reverse order index array.
210         */
211         pSwap = (N_ORG == BLKSIZE) ? swap_l : swap_s;
212 #endif
213
214
215 #if REAL_SEQUENCE
216         /* unpacking the sequence */
217         t_real = x_real[0];
218         t_imag = x_imag[0];
219         x_real[0] = (FLOAT) (t_real+t_imag);
220         x_imag[0] = 0.0;
221         x_real[N] = (FLOAT) (t_real-t_imag);
222         x_imag[N] = 0.0;
223
224         k = off = BLKSIZE/N_ORG;
225         for (i = 1;  i < N/2;  i++)
226         {
227 #if REORDER_DATA
228 #define         a       i
229 #define         b       (N-i)
230 #else
231                 a = pSwap[i];
232                 b = pSwap[N-i];
233 #endif
234                 t1 = x_real[a] + x_real[b];
235                 t2 = x_real[a] - x_real[b];
236                 t3 = x_imag[a] + x_imag[b];
237                 t4 = x_imag[a] - x_imag[b];
238                 t5 = t2*w_imag[k] + t3*w_real[k];
239                 t6 = t3*w_imag[k] - t2*w_real[k];
240
241                 x_real[a] = (FLOAT) (t1+t5) / 2.0;
242                 x_imag[a] = (FLOAT) (t6+t4) / 2.0;
243                 x_real[b] = (FLOAT) (t1-t5) / 2.0;
244                 x_imag[b] = (FLOAT) (t6-t4) / 2.0;
245
246                 k += off;
247         }
248         /* x_real[N/2] doesn«t change */
249         /* x_real[N/2] changes the sign in case of a normal fft */
250 #if (NORMAL_FFT)
251 #if REORDER_DATA
252         x_imag[i] = -x_imag[i];
253 #else
254         x_imag[pSwap[i]] *= -1.0;
255 #endif   /* REORDER_DATA */
256 #endif   /* NORMAL_FFT */
257         N = N_ORG;
258 #endif   /* REAL_SEQUENCE */
259
260
261         /* calculating the energy and phase, phi */
262 #if USED_VALUES_ONLY
263         if (N == BLKSIZE)
264         {
265                 min_i = LONG_FFT_MIN_IDX;
266                 max_i = LONG_FFT_MAX_IDX;
267         }
268         else
269         {
270                 min_i = SHORT_FFT_MIN_IDX;
271                 max_i = SHORT_FFT_MAX_IDX;
272         }
273 #endif
274 #if REAL_SEQUENCE
275         for (i = 0;  i <= N/2;  i++)
276 #else
277         for (i = 0;  i < N;  i++)
278 #endif
279         {
280 #if REORDER_DATA
281 #define         a       i
282 #else
283                 a = pSwap[i];
284 #endif
285                 energy[i] = x_real[a]*x_real[a] + x_imag[a]*x_imag[a];
286                 if(energy[i] <= 0.0005)
287                 {
288                         energy[i] = (FLOAT) 0.0005;        /* keep the identity */
289                         x_real[a] = (FLOAT) sqrt(0.0005);  /* energy[i] * cos(phi[i]) == x_real[i] */
290                         x_imag[a] = 0.0;                   /* energy[i] * sin(phi[i]) == x_imag[i] */
291                 }
292 #if USED_VALUES_ONLY
293                 if (i >= min_i  &&  i <= max_i)
294 #endif
295                 phi[i] = (FLOAT) atan2((double) x_imag[a], (double) x_real[a]);
296         }
297 }
298
299
300