1 /************************* MPEG-2 NBC Audio Decoder **************************
2  *                                                                           *
3 "This software module was originally developed by
4 AT&T, Dolby Laboratories, Fraunhofer Gesellschaft IIS in the course of
5 development of the MPEG-2 NBC/MPEG-4 Audio standard ISO/IEC 13818-7,
6 14496-1,2 and 3. This software module is an implementation of a part of one or more
7 MPEG-2 NBC/MPEG-4 Audio tools as specified by the MPEG-2 NBC/MPEG-4
8 Audio standard. ISO/IEC  gives users of the MPEG-2 NBC/MPEG-4 Audio
9 standards free license to this software module or modifications thereof for use in
10 hardware or software products claiming conformance to the MPEG-2 NBC/MPEG-4
11 Audio  standards. Those intending to use this software module in hardware or
12 software products are advised that this use may infringe existing patents.
13 The original developer of this software module and his/her company, the subsequent
14 editors and their companies, and ISO/IEC have no liability for use of this software
15 module or modifications thereof in an implementation. Copyright is not released for
16 non MPEG-2 NBC/MPEG-4 Audio conforming products.The original developer
17 retains full right to use the code for his/her  own purpose, assign or donate the
18 code to a third party and to inhibit third party from using the code for non
19 MPEG-2 NBC/MPEG-4 Audio conforming products. This copyright notice must
20 be included in all copies or derivative works."
21 Copyright(c)1996.
22  *                                                                           *
23  ****************************************************************************/
24 /*
25  * $Id: filtbank.c,v 1.14 2012/03/01 18:34:17 knik Exp $
26  */
27 
28 /*
29  * CHANGES:
30  *  2001/01/17: menno: Added frequency cut off filter.
31  *
32  */
33 
34 #include <math.h>
35 #include <stdio.h>
36 #include <stdlib.h>
37 
38 #include "coder.h"
39 #include "filtbank.h"
40 #include "frame.h"
41 #include "fft.h"
42 #include "util.h"
43 
44 #define  TWOPI       2*M_PI
45 
46 
47 static void		CalculateKBDWindow	( double* win, double alpha, int length );
48 static double	Izero				( double x);
49 static void		MDCT				( FFT_Tables *fft_tables, double *data, int N );
50 static void		IMDCT				( FFT_Tables *fft_tables, double *data, int N );
51 
52 
53 
FilterBankInit(faacEncStruct * hEncoder)54 void FilterBankInit(faacEncStruct* hEncoder)
55 {
56     unsigned int i, channel;
57 
58     for (channel = 0; channel < hEncoder->numChannels; channel++) {
59         hEncoder->freqBuff[channel] = (double*)AllocMemory(2*FRAME_LEN*sizeof(double));
60         hEncoder->overlapBuff[channel] = (double*)AllocMemory(FRAME_LEN*sizeof(double));
61         SetMemory(hEncoder->overlapBuff[channel], 0, FRAME_LEN*sizeof(double));
62     }
63 
64     hEncoder->sin_window_long = (double*)AllocMemory(BLOCK_LEN_LONG*sizeof(double));
65     hEncoder->sin_window_short = (double*)AllocMemory(BLOCK_LEN_SHORT*sizeof(double));
66     hEncoder->kbd_window_long = (double*)AllocMemory(BLOCK_LEN_LONG*sizeof(double));
67     hEncoder->kbd_window_short = (double*)AllocMemory(BLOCK_LEN_SHORT*sizeof(double));
68 
69     for( i=0; i<BLOCK_LEN_LONG; i++ )
70         hEncoder->sin_window_long[i] = sin((M_PI/(2*BLOCK_LEN_LONG)) * (i + 0.5));
71     for( i=0; i<BLOCK_LEN_SHORT; i++ )
72         hEncoder->sin_window_short[i] = sin((M_PI/(2*BLOCK_LEN_SHORT)) * (i + 0.5));
73 
74     CalculateKBDWindow(hEncoder->kbd_window_long, 4, BLOCK_LEN_LONG*2);
75     CalculateKBDWindow(hEncoder->kbd_window_short, 6, BLOCK_LEN_SHORT*2);
76 }
77 
FilterBankEnd(faacEncStruct * hEncoder)78 void FilterBankEnd(faacEncStruct* hEncoder)
79 {
80     unsigned int channel;
81 
82     for (channel = 0; channel < hEncoder->numChannels; channel++) {
83         if (hEncoder->freqBuff[channel]) FreeMemory(hEncoder->freqBuff[channel]);
84         if (hEncoder->overlapBuff[channel]) FreeMemory(hEncoder->overlapBuff[channel]);
85     }
86 
87     if (hEncoder->sin_window_long) FreeMemory(hEncoder->sin_window_long);
88     if (hEncoder->sin_window_short) FreeMemory(hEncoder->sin_window_short);
89     if (hEncoder->kbd_window_long) FreeMemory(hEncoder->kbd_window_long);
90     if (hEncoder->kbd_window_short) FreeMemory(hEncoder->kbd_window_short);
91 }
92 
FilterBank(faacEncStruct * hEncoder,CoderInfo * coderInfo,double * p_in_data,double * p_out_mdct,double * p_overlap,int overlap_select)93 void FilterBank(faacEncStruct* hEncoder,
94                 CoderInfo *coderInfo,
95                 double *p_in_data,
96                 double *p_out_mdct,
97                 double *p_overlap,
98                 int overlap_select)
99 {
100     double *p_o_buf, *first_window, *second_window;
101     double *transf_buf;
102     int k, i;
103     int block_type = coderInfo->block_type;
104 
105     transf_buf = (double*)AllocMemory(2*BLOCK_LEN_LONG*sizeof(double));
106 
107     /* create / shift old values */
108     /* We use p_overlap here as buffer holding the last frame time signal*/
109     if(overlap_select != MNON_OVERLAPPED) {
110         memcpy(transf_buf, p_overlap, FRAME_LEN*sizeof(double));
111         memcpy(transf_buf+BLOCK_LEN_LONG, p_in_data, FRAME_LEN*sizeof(double));
112         memcpy(p_overlap, p_in_data, FRAME_LEN*sizeof(double));
113     } else {
114         memcpy(transf_buf, p_in_data, 2*FRAME_LEN*sizeof(double));
115     }
116 
117     /*  Window shape processing */
118     if(overlap_select != MNON_OVERLAPPED) {
119         switch (coderInfo->prev_window_shape) {
120         case SINE_WINDOW:
121             if ( (block_type == ONLY_LONG_WINDOW) || (block_type == LONG_SHORT_WINDOW))
122                 first_window = hEncoder->sin_window_long;
123             else
124                 first_window = hEncoder->sin_window_short;
125             break;
126         default:
127         case KBD_WINDOW:
128             if ( (block_type == ONLY_LONG_WINDOW) || (block_type == LONG_SHORT_WINDOW))
129                 first_window = hEncoder->kbd_window_long;
130             else
131                 first_window = hEncoder->kbd_window_short;
132             break;
133         }
134 
135         switch (coderInfo->window_shape){
136         case SINE_WINDOW:
137         default:
138             if ( (block_type == ONLY_LONG_WINDOW) || (block_type == SHORT_LONG_WINDOW))
139                 second_window = hEncoder->sin_window_long;
140             else
141                 second_window = hEncoder->sin_window_short;
142             break;
143         case KBD_WINDOW:
144             if ( (block_type == ONLY_LONG_WINDOW) || (block_type == SHORT_LONG_WINDOW))
145                 second_window = hEncoder->kbd_window_long;
146             else
147                 second_window = hEncoder->kbd_window_short;
148             break;
149         }
150     } else {
151         /* Always long block and sine window for LTP */
152         first_window = hEncoder->sin_window_long;
153         second_window = hEncoder->sin_window_long;
154     }
155 
156     /* Set ptr to transf-Buffer */
157     p_o_buf = transf_buf;
158 
159     /* Separate action for each Block Type */
160     switch (block_type) {
161     case ONLY_LONG_WINDOW :
162         for ( i = 0 ; i < BLOCK_LEN_LONG ; i++){
163             p_out_mdct[i] = p_o_buf[i] * first_window[i];
164             p_out_mdct[i+BLOCK_LEN_LONG] = p_o_buf[i+BLOCK_LEN_LONG] * second_window[BLOCK_LEN_LONG-i-1];
165         }
166         MDCT( &hEncoder->fft_tables, p_out_mdct, 2*BLOCK_LEN_LONG );
167         break;
168 
169     case LONG_SHORT_WINDOW :
170         for ( i = 0 ; i < BLOCK_LEN_LONG ; i++)
171             p_out_mdct[i] = p_o_buf[i] * first_window[i];
172         memcpy(p_out_mdct+BLOCK_LEN_LONG,p_o_buf+BLOCK_LEN_LONG,NFLAT_LS*sizeof(double));
173         for ( i = 0 ; i < BLOCK_LEN_SHORT ; i++)
174             p_out_mdct[i+BLOCK_LEN_LONG+NFLAT_LS] = p_o_buf[i+BLOCK_LEN_LONG+NFLAT_LS] * second_window[BLOCK_LEN_SHORT-i-1];
175         SetMemory(p_out_mdct+BLOCK_LEN_LONG+NFLAT_LS+BLOCK_LEN_SHORT,0,NFLAT_LS*sizeof(double));
176         MDCT( &hEncoder->fft_tables, p_out_mdct, 2*BLOCK_LEN_LONG );
177         break;
178 
179     case SHORT_LONG_WINDOW :
180         SetMemory(p_out_mdct,0,NFLAT_LS*sizeof(double));
181         for ( i = 0 ; i < BLOCK_LEN_SHORT ; i++)
182             p_out_mdct[i+NFLAT_LS] = p_o_buf[i+NFLAT_LS] * first_window[i];
183         memcpy(p_out_mdct+NFLAT_LS+BLOCK_LEN_SHORT,p_o_buf+NFLAT_LS+BLOCK_LEN_SHORT,NFLAT_LS*sizeof(double));
184         for ( i = 0 ; i < BLOCK_LEN_LONG ; i++)
185             p_out_mdct[i+BLOCK_LEN_LONG] = p_o_buf[i+BLOCK_LEN_LONG] * second_window[BLOCK_LEN_LONG-i-1];
186         MDCT( &hEncoder->fft_tables, p_out_mdct, 2*BLOCK_LEN_LONG );
187         break;
188 
189     case ONLY_SHORT_WINDOW :
190         p_o_buf += NFLAT_LS;
191         for ( k=0; k < MAX_SHORT_WINDOWS; k++ ) {
192             for ( i = 0 ; i < BLOCK_LEN_SHORT ; i++ ){
193                 p_out_mdct[i] = p_o_buf[i] * first_window[i];
194                 p_out_mdct[i+BLOCK_LEN_SHORT] = p_o_buf[i+BLOCK_LEN_SHORT] * second_window[BLOCK_LEN_SHORT-i-1];
195             }
196             MDCT( &hEncoder->fft_tables, p_out_mdct, 2*BLOCK_LEN_SHORT );
197             p_out_mdct += BLOCK_LEN_SHORT;
198             p_o_buf += BLOCK_LEN_SHORT;
199             first_window = second_window;
200         }
201         break;
202     }
203 
204     if (transf_buf) FreeMemory(transf_buf);
205 }
206 
IFilterBank(faacEncStruct * hEncoder,CoderInfo * coderInfo,double * p_in_data,double * p_out_data,double * p_overlap,int overlap_select)207 void IFilterBank(faacEncStruct* hEncoder,
208                  CoderInfo *coderInfo,
209                  double *p_in_data,
210                  double *p_out_data,
211                  double *p_overlap,
212                  int overlap_select)
213 {
214     double *o_buf, *transf_buf, *overlap_buf;
215     double *first_window, *second_window;
216 
217     double  *fp;
218     int k, i;
219     int block_type = coderInfo->block_type;
220 
221     transf_buf = (double*)AllocMemory(2*BLOCK_LEN_LONG*sizeof(double));
222     overlap_buf = (double*)AllocMemory(2*BLOCK_LEN_LONG*sizeof(double));
223 
224     /*  Window shape processing */
225     if (overlap_select != MNON_OVERLAPPED) {
226 //      switch (coderInfo->prev_window_shape){
227 //      case SINE_WINDOW:
228             if ( (block_type == ONLY_LONG_WINDOW) || (block_type == LONG_SHORT_WINDOW))
229                 first_window = hEncoder->sin_window_long;
230             else
231                 first_window = hEncoder->sin_window_short;
232 //          break;
233 //      case KBD_WINDOW:
234 //          if ( (block_type == ONLY_LONG_WINDOW) || (block_type == LONG_SHORT_WINDOW))
235 //              first_window = hEncoder->kbd_window_long;
236 //          else
237 //              first_window = hEncoder->kbd_window_short;
238 //          break;
239 //      }
240 
241 //      switch (coderInfo->window_shape){
242 //      case SINE_WINDOW:
243             if ( (block_type == ONLY_LONG_WINDOW) || (block_type == SHORT_LONG_WINDOW))
244                 second_window = hEncoder->sin_window_long;
245             else
246                 second_window = hEncoder->sin_window_short;
247 //          break;
248 //      case KBD_WINDOW:
249 //          if ( (block_type == ONLY_LONG_WINDOW) || (block_type == SHORT_LONG_WINDOW))
250 //              second_window = hEncoder->kbd_window_long;
251 //          else
252 //              second_window = hEncoder->kbd_window_short;
253 //          break;
254 //      }
255     } else {
256         /* Always long block and sine window for LTP */
257         first_window  = hEncoder->sin_window_long;
258         second_window = hEncoder->sin_window_long;
259     }
260 
261     /* Assemble overlap buffer */
262     memcpy(overlap_buf,p_overlap,BLOCK_LEN_LONG*sizeof(double));
263     o_buf = overlap_buf;
264 
265     /* Separate action for each Block Type */
266     switch( block_type ) {
267     case ONLY_LONG_WINDOW :
268         memcpy(transf_buf, p_in_data,BLOCK_LEN_LONG*sizeof(double));
269         IMDCT( &hEncoder->fft_tables, transf_buf, 2*BLOCK_LEN_LONG );
270         for ( i = 0 ; i < BLOCK_LEN_LONG ; i++)
271             transf_buf[i] *= first_window[i];
272         if (overlap_select != MNON_OVERLAPPED) {
273             for ( i = 0 ; i < BLOCK_LEN_LONG; i++ ){
274                 o_buf[i] += transf_buf[i];
275                 o_buf[i+BLOCK_LEN_LONG] = transf_buf[i+BLOCK_LEN_LONG] * second_window[BLOCK_LEN_LONG-i-1];
276             }
277         } else { /* overlap_select == NON_OVERLAPPED */
278             for ( i = 0 ; i < BLOCK_LEN_LONG; i++ )
279                 transf_buf[i+BLOCK_LEN_LONG] *= second_window[BLOCK_LEN_LONG-i-1];
280         }
281         break;
282 
283     case LONG_SHORT_WINDOW :
284         memcpy(transf_buf, p_in_data,BLOCK_LEN_LONG*sizeof(double));
285         IMDCT( &hEncoder->fft_tables, transf_buf, 2*BLOCK_LEN_LONG );
286         for ( i = 0 ; i < BLOCK_LEN_LONG ; i++)
287             transf_buf[i] *= first_window[i];
288         if (overlap_select != MNON_OVERLAPPED) {
289             for ( i = 0 ; i < BLOCK_LEN_LONG; i++ )
290                 o_buf[i] += transf_buf[i];
291             memcpy(o_buf+BLOCK_LEN_LONG,transf_buf+BLOCK_LEN_LONG,NFLAT_LS*sizeof(double));
292             for ( i = 0 ; i < BLOCK_LEN_SHORT ; i++)
293                 o_buf[i+BLOCK_LEN_LONG+NFLAT_LS] = transf_buf[i+BLOCK_LEN_LONG+NFLAT_LS] * second_window[BLOCK_LEN_SHORT-i-1];
294             SetMemory(o_buf+BLOCK_LEN_LONG+NFLAT_LS+BLOCK_LEN_SHORT,0,NFLAT_LS*sizeof(double));
295         } else { /* overlap_select == NON_OVERLAPPED */
296             for ( i = 0 ; i < BLOCK_LEN_SHORT ; i++)
297                 transf_buf[i+BLOCK_LEN_LONG+NFLAT_LS] *= second_window[BLOCK_LEN_SHORT-i-1];
298             SetMemory(transf_buf+BLOCK_LEN_LONG+NFLAT_LS+BLOCK_LEN_SHORT,0,NFLAT_LS*sizeof(double));
299         }
300         break;
301 
302     case SHORT_LONG_WINDOW :
303         memcpy(transf_buf, p_in_data,BLOCK_LEN_LONG*sizeof(double));
304         IMDCT( &hEncoder->fft_tables, transf_buf, 2*BLOCK_LEN_LONG );
305         for ( i = 0 ; i < BLOCK_LEN_SHORT ; i++)
306             transf_buf[i+NFLAT_LS] *= first_window[i];
307         if (overlap_select != MNON_OVERLAPPED) {
308             for ( i = 0 ; i < BLOCK_LEN_SHORT; i++ )
309                 o_buf[i+NFLAT_LS] += transf_buf[i+NFLAT_LS];
310             memcpy(o_buf+BLOCK_LEN_SHORT+NFLAT_LS,transf_buf+BLOCK_LEN_SHORT+NFLAT_LS,NFLAT_LS*sizeof(double));
311             for ( i = 0 ; i < BLOCK_LEN_LONG ; i++)
312                 o_buf[i+BLOCK_LEN_LONG] = transf_buf[i+BLOCK_LEN_LONG] * second_window[BLOCK_LEN_LONG-i-1];
313         } else { /* overlap_select == NON_OVERLAPPED */
314             SetMemory(transf_buf,0,NFLAT_LS*sizeof(double));
315             for ( i = 0 ; i < BLOCK_LEN_LONG ; i++)
316                 transf_buf[i+BLOCK_LEN_LONG] *= second_window[BLOCK_LEN_LONG-i-1];
317         }
318         break;
319 
320     case ONLY_SHORT_WINDOW :
321         if (overlap_select != MNON_OVERLAPPED) {
322             fp = o_buf + NFLAT_LS;
323         } else { /* overlap_select == NON_OVERLAPPED */
324             fp = transf_buf;
325         }
326         for ( k=0; k < MAX_SHORT_WINDOWS; k++ ) {
327             memcpy(transf_buf,p_in_data,BLOCK_LEN_SHORT*sizeof(double));
328             IMDCT( &hEncoder->fft_tables, transf_buf, 2*BLOCK_LEN_SHORT );
329             p_in_data += BLOCK_LEN_SHORT;
330             if (overlap_select != MNON_OVERLAPPED) {
331                 for ( i = 0 ; i < BLOCK_LEN_SHORT ; i++){
332                     transf_buf[i] *= first_window[i];
333                     fp[i] += transf_buf[i];
334                     fp[i+BLOCK_LEN_SHORT] = transf_buf[i+BLOCK_LEN_SHORT] * second_window[BLOCK_LEN_SHORT-i-1];
335                 }
336                 fp += BLOCK_LEN_SHORT;
337             } else { /* overlap_select == NON_OVERLAPPED */
338                 for ( i = 0 ; i < BLOCK_LEN_SHORT ; i++){
339                     fp[i] *= first_window[i];
340                     fp[i+BLOCK_LEN_SHORT] *= second_window[BLOCK_LEN_SHORT-i-1];
341                 }
342                 fp += 2*BLOCK_LEN_SHORT;
343             }
344             first_window = second_window;
345         }
346         SetMemory(o_buf+BLOCK_LEN_LONG+NFLAT_LS+BLOCK_LEN_SHORT,0,NFLAT_LS*sizeof(double));
347         break;
348     }
349 
350     if (overlap_select != MNON_OVERLAPPED)
351         memcpy(p_out_data,o_buf,BLOCK_LEN_LONG*sizeof(double));
352     else  /* overlap_select == NON_OVERLAPPED */
353         memcpy(p_out_data,transf_buf,2*BLOCK_LEN_LONG*sizeof(double));
354 
355     /* save unused output data */
356     memcpy(p_overlap,o_buf+BLOCK_LEN_LONG,BLOCK_LEN_LONG*sizeof(double));
357 
358     if (overlap_buf) FreeMemory(overlap_buf);
359     if (transf_buf) FreeMemory(transf_buf);
360 }
361 
specFilter(double * freqBuff,int sampleRate,int lowpassFreq,int specLen)362 void specFilter(double *freqBuff,
363                 int sampleRate,
364                 int lowpassFreq,
365                 int specLen
366                 )
367 {
368     int lowpass,xlowpass;
369 
370     /* calculate the last line which is not zero */
371     lowpass = (lowpassFreq * specLen) / (sampleRate>>1) + 1;
372     xlowpass = (lowpass < specLen) ? lowpass : specLen ;
373 
374     SetMemory(freqBuff+xlowpass,0,(specLen-xlowpass)*sizeof(double));
375 }
376 
Izero(double x)377 static double Izero(double x)
378 {
379     const double IzeroEPSILON = 1E-41;  /* Max error acceptable in Izero */
380     double sum, u, halfx, temp;
381     int n;
382 
383     sum = u = n = 1;
384     halfx = x/2.0;
385     do {
386         temp = halfx/(double)n;
387         n += 1;
388         temp *= temp;
389         u *= temp;
390         sum += u;
391     } while (u >= IzeroEPSILON*sum);
392 
393     return(sum);
394 }
395 
CalculateKBDWindow(double * win,double alpha,int length)396 static void CalculateKBDWindow(double* win, double alpha, int length)
397 {
398     int i;
399     double IBeta;
400     double tmp;
401     double sum = 0.0;
402 
403     alpha *= M_PI;
404     IBeta = 1.0/Izero(alpha);
405 
406     /* calculate lower half of Kaiser Bessel window */
407     for(i=0; i<(length>>1); i++) {
408         tmp = 4.0*(double)i/(double)length - 1.0;
409         win[i] = Izero(alpha*sqrt(1.0-tmp*tmp))*IBeta;
410         sum += win[i];
411     }
412 
413     sum = 1.0/sum;
414     tmp = 0.0;
415 
416     /* calculate lower half of window */
417     for(i=0; i<(length>>1); i++) {
418         tmp += win[i];
419         win[i] = sqrt(tmp*sum);
420     }
421 }
422 
MDCT(FFT_Tables * fft_tables,double * data,int N)423 static void MDCT( FFT_Tables *fft_tables, double *data, int N )
424 {
425     double *xi, *xr;
426     double tempr, tempi, c, s, cold, cfreq, sfreq; /* temps for pre and post twiddle */
427     double freq = TWOPI / N;
428     double cosfreq8, sinfreq8;
429     int i, n;
430 
431     xi = (double*)AllocMemory((N >> 2)*sizeof(double));
432     xr = (double*)AllocMemory((N >> 2)*sizeof(double));
433 
434     /* prepare for recurrence relation in pre-twiddle */
435     cfreq = cos (freq);
436     sfreq = sin (freq);
437     cosfreq8 = cos (freq * 0.125);
438     sinfreq8 = sin (freq * 0.125);
439     c = cosfreq8;
440     s = sinfreq8;
441 
442     for (i = 0; i < (N >> 2); i++) {
443         /* calculate real and imaginary parts of g(n) or G(p) */
444         n = (N >> 1) - 1 - 2 * i;
445 
446         if (i < (N >> 3))
447             tempr = data [(N >> 2) + n] + data [N + (N >> 2) - 1 - n]; /* use second form of e(n) for n = N / 2 - 1 - 2i */
448         else
449             tempr = data [(N >> 2) + n] - data [(N >> 2) - 1 - n]; /* use first form of e(n) for n = N / 2 - 1 - 2i */
450 
451         n = 2 * i;
452         if (i < (N >> 3))
453             tempi = data [(N >> 2) + n] - data [(N >> 2) - 1 - n]; /* use first form of e(n) for n=2i */
454         else
455             tempi = data [(N >> 2) + n] + data [N + (N >> 2) - 1 - n]; /* use second form of e(n) for n=2i*/
456 
457         /* calculate pre-twiddled FFT input */
458         xr[i] = tempr * c + tempi * s;
459         xi[i] = tempi * c - tempr * s;
460 
461         /* use recurrence to prepare cosine and sine for next value of i */
462         cold = c;
463         c = c * cfreq - s * sfreq;
464         s = s * cfreq + cold * sfreq;
465     }
466 
467     /* Perform in-place complex FFT of length N/4 */
468     switch (N) {
469     case BLOCK_LEN_SHORT * 2:
470         fft( fft_tables, xr, xi, 6);
471         break;
472     case BLOCK_LEN_LONG * 2:
473         fft( fft_tables, xr, xi, 9);
474     }
475 
476     /* prepare for recurrence relations in post-twiddle */
477     c = cosfreq8;
478     s = sinfreq8;
479 
480     /* post-twiddle FFT output and then get output data */
481     for (i = 0; i < (N >> 2); i++) {
482         /* get post-twiddled FFT output  */
483         tempr = 2. * (xr[i] * c + xi[i] * s);
484         tempi = 2. * (xi[i] * c - xr[i] * s);
485 
486         /* fill in output values */
487         data [2 * i] = -tempr;   /* first half even */
488         data [(N >> 1) - 1 - 2 * i] = tempi;  /* first half odd */
489         data [(N >> 1) + 2 * i] = -tempi;  /* second half even */
490         data [N - 1 - 2 * i] = tempr;  /* second half odd */
491 
492         /* use recurrence to prepare cosine and sine for next value of i */
493         cold = c;
494         c = c * cfreq - s * sfreq;
495         s = s * cfreq + cold * sfreq;
496     }
497 
498     if (xr) FreeMemory(xr);
499     if (xi) FreeMemory(xi);
500 }
501 
IMDCT(FFT_Tables * fft_tables,double * data,int N)502 static void IMDCT( FFT_Tables *fft_tables, double *data, int N)
503 {
504     double *xi, *xr;
505     double tempr, tempi, c, s, cold, cfreq, sfreq; /* temps for pre and post twiddle */
506     double freq = 2.0 * M_PI / N;
507     double fac, cosfreq8, sinfreq8;
508     int i;
509 
510     xi = (double*)AllocMemory((N >> 2)*sizeof(double));
511     xr = (double*)AllocMemory((N >> 2)*sizeof(double));
512 
513     /* Choosing to allocate 2/N factor to Inverse Xform! */
514     fac = 2. / N; /* remaining 2/N from 4/N IFFT factor */
515 
516     /* prepare for recurrence relation in pre-twiddle */
517     cfreq = cos (freq);
518     sfreq = sin (freq);
519     cosfreq8 = cos (freq * 0.125);
520     sinfreq8 = sin (freq * 0.125);
521     c = cosfreq8;
522     s = sinfreq8;
523 
524     for (i = 0; i < (N >> 2); i++) {
525         /* calculate real and imaginary parts of g(n) or G(p) */
526         tempr = -data[2 * i];
527         tempi = data[(N >> 1) - 1 - 2 * i];
528 
529         /* calculate pre-twiddled FFT input */
530         xr[i] = tempr * c - tempi * s;
531         xi[i] = tempi * c + tempr * s;
532 
533         /* use recurrence to prepare cosine and sine for next value of i */
534         cold = c;
535         c = c * cfreq - s * sfreq;
536         s = s * cfreq + cold * sfreq;
537     }
538 
539     /* Perform in-place complex IFFT of length N/4 */
540     switch (N) {
541     case BLOCK_LEN_SHORT * 2:
542         ffti( fft_tables, xr, xi, 6);
543         break;
544     case BLOCK_LEN_LONG * 2:
545         ffti( fft_tables, xr, xi, 9);
546     }
547 
548     /* prepare for recurrence relations in post-twiddle */
549     c = cosfreq8;
550     s = sinfreq8;
551 
552     /* post-twiddle FFT output and then get output data */
553     for (i = 0; i < (N >> 2); i++) {
554 
555         /* get post-twiddled FFT output  */
556         tempr = fac * (xr[i] * c - xi[i] * s);
557         tempi = fac * (xi[i] * c + xr[i] * s);
558 
559         /* fill in output values */
560         data [(N >> 1) + (N >> 2) - 1 - 2 * i] = tempr;
561         if (i < (N >> 3))
562             data [(N >> 1) + (N >> 2) + 2 * i] = tempr;
563         else
564             data [2 * i - (N >> 2)] = -tempr;
565 
566         data [(N >> 2) + 2 * i] = tempi;
567         if (i < (N >> 3))
568             data [(N >> 2) - 1 - 2 * i] = -tempi;
569         else
570             data [(N >> 2) + N - 1 - 2*i] = tempi;
571 
572         /* use recurrence to prepare cosine and sine for next value of i */
573         cold = c;
574         c = c * cfreq - s * sfreq;
575         s = s * cfreq + cold * sfreq;
576     }
577 
578     if (xr) FreeMemory(xr);
579     if (xi) FreeMemory(xi);
580 }
581