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