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