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