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 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
