1 /* ------------------------------------------------------------------
2 * Copyright (C) 1998-2009 PacketVideo
3 *
4 * Licensed under the Apache License, Version 2.0 (the "License");
5 * you may not use this file except in compliance with the License.
6 * You may obtain a copy of the License at
7 *
8 * http://www.apache.org/licenses/LICENSE-2.0
9 *
10 * Unless required by applicable law or agreed to in writing, software
11 * distributed under the License is distributed on an "AS IS" BASIS,
12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either
13 * express or implied.
14 * See the License for the specific language governing permissions
15 * and limitations under the License.
16 * -------------------------------------------------------------------
17 */
18 /*
19
20 Pathname: imdct_fxp.c
21 Funtions: imdct_fxp
22
23 ------------------------------------------------------------------------------
24 INPUT AND OUTPUT DEFINITIONS
25
26 Inputs:
27
28 data_quant = Input vector, with quantized spectral lines:
29 type Int32
30
31 freq_2_time_buffer = Scratch memory used for in-place FFT calculation,
32 min size required 1024,
33 type Int32
34
35 n = Length of input vector "data_quant". Currently 256 or 2048
36 type const Int
37
38 Q_format = Q_format of the input vector "data_quant"
39 type Int
40
41 max = Maximum value inside input vector "data_quant"
42 type Int32
43
44 Local Stores/Buffers/Pointers Needed:
45 None
46
47 Global Stores/Buffers/Pointers Needed:
48 None
49
50 Outputs:
51 shift = shift factor to reflect scaling introduced by IFFT and imdct_fxp,
52
53 Pointers and Buffers Modified:
54 Results are return in "Data_Int_precision"
55
56 Local Stores Modified:
57 None
58
59 Global Stores Modified:
60 None
61
62 ------------------------------------------------------------------------------
63 FUNCTION DESCRIPTION
64
65 The IMDCT is a linear orthogonal lapped transform, based on the idea of
66 time domain aliasing cancellation (TDAC).
67 IMDCT is critically sampled, which means that though it is 50% overlapped,
68 a sequence data after IMDCT has the same number of coefficients as samples
69 before the transform (after overlap-and-add). This means, that a single
70 block of IMDCT data does not correspond to the original block on which the
71 IMDCT was performed. When subsequent blocks of inverse transformed data
72 are added (still using 50% overlap), the errors introduced by the
73 transform cancels out.Thanks to the overlapping feature, the IMDCT is very
74 useful for quantization. It effectively removes the otherwise easily
75 detectable blocking artifact between transform blocks.
76
77 N = twice the length of input vector X
78 y = vector of length N, will hold fixed point IDCT
79 p = 0:1:N-1
80
81 2 N/2-1
82 y(p) = --- SUM X(m)*cos(pi/(2*N)*(2*p+1+N/2)*(2*m+1))
83 N m=0
84
85 The window that completes the TDAC is applied before calling this function.
86 The IMDCT can be calculated using an IFFT, for this, the IMDCT need be
87 rewritten as an odd-time odd-frequency discrete Fourier transform. Thus,
88 the IMDCT can be calculated using only one n/4 point FFT and some pre and
89 post-rotation of the sample points.
90
91
92 where X(k) is the input with N frequency lines
93
94 X(k) ----------------------------
95 |
96 |
97 Pre-rotation by exp(j(2pi/N)(k+1/8))
98 |
99 |
100 N/4- point IFFT
101 |
102 |
103 Post-rotation by exp(j(2pi/N)(n+1/8))
104 |
105 |
106 ------------- x(n) In the time domain
107
108
109 ------------------------------------------------------------------------------
110 REQUIREMENTS
111
112 This function should provide a fixed point IMDCT with an average
113 quantization error less than 1 % (variance and mean).
114
115 ------------------------------------------------------------------------------
116 REFERENCES
117
118 [1] Analysis/Synthesis Filter Bank design based on time domain
119 aliasing cancellation
120 Jhon Princen, et. al.
121 IEEE Transactions on ASSP, vol ASSP-34, No. 5 October 1986
122 Pg 1153 - 1161
123
124 [2] Regular FFT-related transform kernels for DCT/DST based
125 polyphase filterbanks
126 Rolf Gluth
127 Proc. ICASSP 1991, pg. 2205 - 2208
128
129 ------------------------------------------------------------------------------
130 PSEUDO-CODE
131
132 Cx, Cy are complex number
133
134
135 exp = log2(n)-1
136
137 FOR ( k=0; k< n/2; k +=2)
138
139 Cx = - data_quant[k] + j data_quant[n/2-1 - k]
140
141 freq_2_time_buffer = Cx * exp(j(2pi/n)(k+1/8))
142
143 ENDFOR
144
145 CALL IFFT( freq_2_time_buffer, n/4)
146
147 MODIFYING( freq_2_time_buffer )
148
149 RETURNING( shift )
150
151 FOR ( k=0; k< n/4; k +=2)
152
153 Cx = freq_2_time_buffer[ k] + j freq_2_time_buffer[ k+1]
154
155 Cy = Cx * exp(j(2pi/n)(k+1/8))
156
157 data_quant[3n/4-1 - k ] = Real(Cy)
158 data_quant[ n/4-1 - k ] = - Imag(Cy)
159 data_quant[3n/4 + k ] = Real(Cy)
160 data_quant[ n/4 + k ] = Imag(Cy)
161
162 ENDFOR
163
164 FOR ( k=n/4; k< n/2; k +=2)
165
166 Cx = freq_2_time_buffer[ k] + j freq_2_time_buffer[ k+1]
167
168 Cy = Cx * exp(j(2pi/n)(k+1/8))
169
170 data_quant[3n/4-1 - k ] = Real(Cy)
171 data_quant[ n/4 + k ] = - Real(Cy)
172 data_quant[5n/4 - k ] = Imag(Cy)
173 data_quant[ n/4 + k ] = Imag(Cy)
174
175 ENDFOR
176
177 MODIFIED data_quant[]
178
179 RETURN (exp - shift)
180
181 ------------------------------------------------------------------------------
182 RESOURCES USED
183 When the code is written for a specific target processor the
184 the resources used should be documented below.
185
186 STACK USAGE: [stack count for this module] + [variable to represent
187 stack usage for each subroutine called]
188
189 where: [stack usage variable] = stack usage for [subroutine
190 name] (see [filename].ext)
191
192 DATA MEMORY USED: x words
193
194 PROGRAM MEMORY USED: x words
195
196 CLOCK CYCLES: [cycle count equation for this module] + [variable
197 used to represent cycle count for each subroutine
198 called]
199
200 where: [cycle count variable] = cycle count for [subroutine
201 name] (see [filename].ext)
202
203 ------------------------------------------------------------------------------
204 */
205
206
207 /*----------------------------------------------------------------------------
208 ; INCLUDES
209 ----------------------------------------------------------------------------*/
210
211 #include "pv_audio_type_defs.h"
212 #include "imdct_fxp.h"
213
214
215 #include "mix_radix_fft.h"
216 #include "digit_reversal_tables.h"
217 #include "fft_rx4.h"
218 #include "inv_short_complex_rot.h"
219 #include "inv_long_complex_rot.h"
220 #include "pv_normalize.h"
221 #include "fxp_mul32.h"
222 #include "aac_mem_funcs.h"
223
224 #include "window_block_fxp.h"
225
226 /*----------------------------------------------------------------------------
227 ; MACROS
228 ; Define module specific macros here
229 ----------------------------------------------------------------------------*/
230
231 /*----------------------------------------------------------------------------
232 ; DEFINES
233 ; Include all pre-processor statements here. Include conditional
234 ; compile variables also.
235 ----------------------------------------------------------------------------*/
236 #define ERROR_IN_FRAME_SIZE 10
237
238 /*----------------------------------------------------------------------------
239 ; LOCAL FUNCTION DEFINITIONS
240 ; Function Prototype declaration
241 ----------------------------------------------------------------------------*/
242
243 /*----------------------------------------------------------------------------
244 ; LOCAL VARIABLE DEFINITIONS
245 ; Variable declaration - defined here and used outside this module
246 ----------------------------------------------------------------------------*/
247
248 /*----------------------------------------------------------------------------
249 ; EXTERNAL FUNCTION REFERENCES
250 ; Declare functions defined elsewhere and referenced in this module
251 ----------------------------------------------------------------------------*/
252
253 /*----------------------------------------------------------------------------
254 ; EXTERNAL VARIABLES REFERENCES
255 ; Declare variables used in this module but defined elsewhere
256 ----------------------------------------------------------------------------*/
257
258
259 /*----------------------------------------------------------------------------
260 ; EXTERNAL GLOBAL STORE/BUFFER/POINTER REFERENCES
261 ; Declare variables used in this module but defined elsewhere
262 ----------------------------------------------------------------------------*/
263
264 /*----------------------------------------------------------------------------
265 ; FUNCTION CODE
266 ----------------------------------------------------------------------------*/
267
268
imdct_fxp(Int32 data_quant[],Int32 freq_2_time_buffer[],const Int n,Int Q_format,Int32 max)269 Int imdct_fxp(Int32 data_quant[],
270 Int32 freq_2_time_buffer[],
271 const Int n,
272 Int Q_format,
273 Int32 max)
274 {
275
276 Int32 exp_jw;
277 Int shift = 0;
278
279 const Int32 *p_rotate;
280 const Int32 *p_rotate_2;
281
282 Int32 *p_data_1;
283 Int32 *p_data_2;
284
285 Int32 temp_re32;
286 Int32 temp_im32;
287
288 Int shift1 = 0;
289 Int32 temp1;
290 Int32 temp2;
291
292 Int k;
293 Int n_2 = n >> 1;
294 Int n_4 = n >> 2;
295
296
297
298 if (max != 0)
299 {
300
301 switch (n)
302 {
303 case SHORT_WINDOW_TYPE:
304 p_rotate = exp_rotation_N_256;
305 shift = 21; /* log2(n)-1 + 14 acomodates 2/N factor */
306 break;
307
308 case LONG_WINDOW_TYPE:
309 p_rotate = exp_rotation_N_2048;
310 shift = 24; /* log2(n)-1 +14 acomodates 2/N factor */
311 break;
312
313 default:
314 /*
315 * There is no defined behavior for a non supported frame
316 * size. By returning a fixed scaling factor, the input will
317 * scaled down and the will be heard as a low level noise
318 */
319 return(ERROR_IN_FRAME_SIZE);
320
321 }
322
323 /*
324 * p_data_1 p_data_2
325 * | |
326 * RIRIRIRIRIRIRIRIRIRIRIRIRIRIRI....RIRIRIRIRIRI
327 * | |
328 *
329 */
330
331 p_data_1 = data_quant; /* uses first half of buffer */
332 p_data_2 = &data_quant[n_2 - 1]; /* uses second half of buffer */
333
334 p_rotate_2 = &p_rotate[n_4-1];
335
336 shift1 = pv_normalize(max) - 1; /* -1 to leave room for addition */
337 Q_format -= (16 - shift1);
338 max = 0;
339
340
341 if (shift1 >= 0)
342 {
343 temp_re32 = *(p_data_1++) << shift1;
344 temp_im32 = *(p_data_2--) << shift1;
345
346 for (k = n_4 >> 1; k != 0; k--)
347 {
348 /*
349 * Real and Imag parts have been swaped to use FFT as IFFT
350 */
351 /*
352 * cos_n + j*sin_n == exp(j(2pi/N)(k+1/8))
353 */
354 exp_jw = *p_rotate++;
355
356 temp1 = cmplx_mul32_by_16(temp_im32, -temp_re32, exp_jw);
357 temp2 = -cmplx_mul32_by_16(temp_re32, temp_im32, exp_jw);
358
359 temp_im32 = *(p_data_1--) << shift1;
360 temp_re32 = *(p_data_2--) << shift1;
361 *(p_data_1++) = temp1;
362 *(p_data_1++) = temp2;
363 max |= (temp1 >> 31) ^ temp1;
364 max |= (temp2 >> 31) ^ temp2;
365
366
367 /*
368 * Real and Imag parts have been swaped to use FFT as IFFT
369 */
370
371 /*
372 * cos_n + j*sin_n == exp(j(2pi/N)(k+1/8))
373 */
374
375 exp_jw = *p_rotate_2--;
376
377 temp1 = cmplx_mul32_by_16(temp_im32, -temp_re32, exp_jw);
378 temp2 = -cmplx_mul32_by_16(temp_re32, temp_im32, exp_jw);
379
380
381 temp_re32 = *(p_data_1++) << shift1;
382 temp_im32 = *(p_data_2--) << shift1;
383
384 *(p_data_2 + 2) = temp1;
385 *(p_data_2 + 3) = temp2;
386 max |= (temp1 >> 31) ^ temp1;
387 max |= (temp2 >> 31) ^ temp2;
388
389 }
390 }
391 else
392 {
393 temp_re32 = *(p_data_1++) >> 1;
394 temp_im32 = *(p_data_2--) >> 1;
395
396 for (k = n_4 >> 1; k != 0; k--)
397 {
398 /*
399 * Real and Imag parts have been swaped to use FFT as IFFT
400 */
401 /*
402 * cos_n + j*sin_n == exp(j(2pi/N)(k+1/8))
403 */
404 exp_jw = *p_rotate++;
405
406 temp1 = cmplx_mul32_by_16(temp_im32, -temp_re32, exp_jw);
407 temp2 = -cmplx_mul32_by_16(temp_re32, temp_im32, exp_jw);
408
409 temp_im32 = *(p_data_1--) >> 1;
410 temp_re32 = *(p_data_2--) >> 1;
411 *(p_data_1++) = temp1;
412 *(p_data_1++) = temp2;
413
414 max |= (temp1 >> 31) ^ temp1;
415 max |= (temp2 >> 31) ^ temp2;
416
417
418 /*
419 * Real and Imag parts have been swaped to use FFT as IFFT
420 */
421
422 /*
423 * cos_n + j*sin_n == exp(j(2pi/N)(k+1/8))
424 */
425 exp_jw = *p_rotate_2--;
426
427 temp1 = cmplx_mul32_by_16(temp_im32, -temp_re32, exp_jw);
428 temp2 = -cmplx_mul32_by_16(temp_re32, temp_im32, exp_jw);
429
430 temp_re32 = *(p_data_1++) >> 1;
431 temp_im32 = *(p_data_2--) >> 1;
432
433 *(p_data_2 + 3) = temp2;
434 *(p_data_2 + 2) = temp1;
435
436 max |= (temp1 >> 31) ^ temp1;
437 max |= (temp2 >> 31) ^ temp2;
438
439 }
440 }
441
442
443 if (n != SHORT_WINDOW_TYPE)
444 {
445
446 shift -= mix_radix_fft(data_quant,
447 &max);
448
449 shift -= inv_long_complex_rot(data_quant,
450 max);
451
452 }
453 else /* n_4 is 64 */
454 {
455
456 shift -= fft_rx4_short(data_quant, &max);
457
458
459 shift -= inv_short_complex_rot(data_quant,
460 freq_2_time_buffer,
461 max);
462
463 pv_memcpy(data_quant,
464 freq_2_time_buffer,
465 SHORT_WINDOW*sizeof(*data_quant));
466 }
467
468 }
469 else
470 {
471 Q_format = ALL_ZEROS_BUFFER;
472 }
473
474 return(shift + Q_format);
475
476 } /* imdct_fxp */
477