1 2 /* ----------------------------------------------------------------------------------------------------------- 3 Software License for The Fraunhofer FDK AAC Codec Library for Android 4 5 Copyright 1995 - 2012 Fraunhofer-Gesellschaft zur Frderung der angewandten Forschung e.V. 6 All rights reserved. 7 8 1. INTRODUCTION 9 The Fraunhofer FDK AAC Codec Library for Android ("FDK AAC Codec") is software that implements 10 the MPEG Advanced Audio Coding ("AAC") encoding and decoding scheme for digital audio. 11 This FDK AAC Codec software is intended to be used on a wide variety of Android devices. 12 13 AAC's HE-AAC and HE-AAC v2 versions are regarded as today's most efficient general perceptual 14 audio codecs. AAC-ELD is considered the best-performing full-bandwidth communications codec by 15 independent studies and is widely deployed. AAC has been standardized by ISO and IEC as part 16 of the MPEG specifications. 17 18 Patent licenses for necessary patent claims for the FDK AAC Codec (including those of Fraunhofer) 19 may be obtained through Via Licensing (www.vialicensing.com) or through the respective patent owners 20 individually for the purpose of encoding or decoding bit streams in products that are compliant with 21 the ISO/IEC MPEG audio standards. Please note that most manufacturers of Android devices already license 22 these patent claims through Via Licensing or directly from the patent owners, and therefore FDK AAC Codec 23 software may already be covered under those patent licenses when it is used for those licensed purposes only. 24 25 Commercially-licensed AAC software libraries, including floating-point versions with enhanced sound quality, 26 are also available from Fraunhofer. Users are encouraged to check the Fraunhofer website for additional 27 applications information and documentation. 28 29 2. COPYRIGHT LICENSE 30 31 Redistribution and use in source and binary forms, with or without modification, are permitted without 32 payment of copyright license fees provided that you satisfy the following conditions: 33 34 You must retain the complete text of this software license in redistributions of the FDK AAC Codec or 35 your modifications thereto in source code form. 36 37 You must retain the complete text of this software license in the documentation and/or other materials 38 provided with redistributions of the FDK AAC Codec or your modifications thereto in binary form. 39 You must make available free of charge copies of the complete source code of the FDK AAC Codec and your 40 modifications thereto to recipients of copies in binary form. 41 42 The name of Fraunhofer may not be used to endorse or promote products derived from this library without 43 prior written permission. 44 45 You may not charge copyright license fees for anyone to use, copy or distribute the FDK AAC Codec 46 software or your modifications thereto. 47 48 Your modified versions of the FDK AAC Codec must carry prominent notices stating that you changed the software 49 and the date of any change. For modified versions of the FDK AAC Codec, the term 50 "Fraunhofer FDK AAC Codec Library for Android" must be replaced by the term 51 "Third-Party Modified Version of the Fraunhofer FDK AAC Codec Library for Android." 52 53 3. NO PATENT LICENSE 54 55 NO EXPRESS OR IMPLIED LICENSES TO ANY PATENT CLAIMS, including without limitation the patents of Fraunhofer, 56 ARE GRANTED BY THIS SOFTWARE LICENSE. Fraunhofer provides no warranty of patent non-infringement with 57 respect to this software. 58 59 You may use this FDK AAC Codec software or modifications thereto only for purposes that are authorized 60 by appropriate patent licenses. 61 62 4. DISCLAIMER 63 64 This FDK AAC Codec software is provided by Fraunhofer on behalf of the copyright holders and contributors 65 "AS IS" and WITHOUT ANY EXPRESS OR IMPLIED WARRANTIES, including but not limited to the implied warranties 66 of merchantability and fitness for a particular purpose. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR 67 CONTRIBUTORS BE LIABLE for any direct, indirect, incidental, special, exemplary, or consequential damages, 68 including but not limited to procurement of substitute goods or services; loss of use, data, or profits, 69 or business interruption, however caused and on any theory of liability, whether in contract, strict 70 liability, or tort (including negligence), arising in any way out of the use of this software, even if 71 advised of the possibility of such damage. 72 73 5. CONTACT INFORMATION 74 75 Fraunhofer Institute for Integrated Circuits IIS 76 Attention: Audio and Multimedia Departments - FDK AAC LL 77 Am Wolfsmantel 33 78 91058 Erlangen, Germany 79 80 www.iis.fraunhofer.de/amm 81 amm-info (at) iis.fraunhofer.de 82 ----------------------------------------------------------------------------------------------------------- */ 83 84 /*************************** Fraunhofer IIS FDK Tools ********************** 85 86 Author(s): Josef Hoepfl, DSP Solutions 87 Description: Fix point FFT 88 89 ******************************************************************************/ 90 91 #include "fft.h" 92 93 #include "fft_rad2.h" 94 #include "FDK_tools_rom.h" 95 96 97 98 99 100 #define F3C(x) STC(x) 101 102 #define C31 (F3C(0x91261468)) /* FL2FXCONST_DBL(-0.86602540) */ 103 104 /* Performs the FFT of length 3 according to the algorithm after winograd. 105 No scaling of the input vector because the scaling is already done in the rotation vector. */ 106 static FORCEINLINE void fft3(FIXP_DBL *RESTRICT pDat) 107 { 108 FIXP_DBL r1,r2; 109 FIXP_DBL s1,s2; 110 /* real part */ 111 r1 = pDat[2] + pDat[4]; 112 r2 = fMult((pDat[2] - pDat[4]), C31); 113 pDat[0] = pDat[0] + r1; 114 r1 = pDat[0] - r1 - (r1>>1); 115 116 /* imaginary part */ 117 s1 = pDat[3] + pDat[5]; 118 s2 = fMult((pDat[3] - pDat[5]), C31); 119 pDat[1] = pDat[1] + s1; 120 s1 = pDat[1] - s1 - (s1>>1); 121 122 /* combination */ 123 pDat[2] = r1 - s2; 124 pDat[4] = r1 + s2; 125 pDat[3] = s1 + r2; 126 pDat[5] = s1 - r2; 127 } 128 129 130 #define F5C(x) STC(x) 131 132 #define C51 (F5C(0x79bc3854)) /* FL2FXCONST_DBL( 0.95105652) */ 133 #define C52 (F5C(0x9d839db0)) /* FL2FXCONST_DBL(-1.53884180/2) */ 134 #define C53 (F5C(0xd18053ce)) /* FL2FXCONST_DBL(-0.36327126) */ 135 #define C54 (F5C(0x478dde64)) /* FL2FXCONST_DBL( 0.55901699) */ 136 #define C55 (F5C(0xb0000001)) /* FL2FXCONST_DBL(-1.25/2) */ 137 138 /* performs the FFT of length 5 according to the algorithm after winograd */ 139 static FORCEINLINE void fft5(FIXP_DBL *RESTRICT pDat) 140 { 141 FIXP_DBL r1,r2,r3,r4; 142 FIXP_DBL s1,s2,s3,s4; 143 FIXP_DBL t; 144 145 /* real part */ 146 r1 = pDat[2] + pDat[8]; 147 r4 = pDat[2] - pDat[8]; 148 r3 = pDat[4] + pDat[6]; 149 r2 = pDat[4] - pDat[6]; 150 t = fMult((r1-r3), C54); 151 r1 = r1 + r3; 152 pDat[0] = pDat[0] + r1; 153 /* Bit shift left because of the constant C55 which was scaled with the factor 0.5 because of the representation of 154 the values as fracts */ 155 r1 = pDat[0] + (fMultDiv2(r1, C55) <<(2)); 156 r3 = r1 - t; 157 r1 = r1 + t; 158 t = fMult((r4 + r2), C51); 159 /* Bit shift left because of the constant C55 which was scaled with the factor 0.5 because of the representation of 160 the values as fracts */ 161 r4 = t + (fMultDiv2(r4, C52) <<(2)); 162 r2 = t + fMult(r2, C53); 163 164 /* imaginary part */ 165 s1 = pDat[3] + pDat[9]; 166 s4 = pDat[3] - pDat[9]; 167 s3 = pDat[5] + pDat[7]; 168 s2 = pDat[5] - pDat[7]; 169 t = fMult((s1 - s3), C54); 170 s1 = s1 + s3; 171 pDat[1] = pDat[1] + s1; 172 /* Bit shift left because of the constant C55 which was scaled with the factor 0.5 because of the representation of 173 the values as fracts */ 174 s1 = pDat[1] + (fMultDiv2(s1, C55) <<(2)); 175 s3 = s1 - t; 176 s1 = s1 + t; 177 t = fMult((s4 + s2), C51); 178 /* Bit shift left because of the constant C55 which was scaled with the factor 0.5 because of the representation of 179 the values as fracts */ 180 s4 = t + (fMultDiv2(s4, C52) <<(2)); 181 s2 = t + fMult(s2, C53); 182 183 /* combination */ 184 pDat[2] = r1 + s2; 185 pDat[8] = r1 - s2; 186 pDat[4] = r3 - s4; 187 pDat[6] = r3 + s4; 188 189 pDat[3] = s1 - r2; 190 pDat[9] = s1 + r2; 191 pDat[5] = s3 + r4; 192 pDat[7] = s3 - r4; 193 } 194 195 196 197 198 #define N3 3 199 #define N5 5 200 #define N6 6 201 #define N15 15 202 203 /* Performs the FFT of length 15. It is split into FFTs of length 3 and length 5. */ 204 static inline void fft15(FIXP_DBL *pInput) 205 { 206 FIXP_DBL aDst[2*N15]; 207 FIXP_DBL aDst1[2*N15]; 208 int i,k,l; 209 210 /* Sort input vector for fft's of length 3 211 input3(0:2) = [input(0) input(5) input(10)]; 212 input3(3:5) = [input(3) input(8) input(13)]; 213 input3(6:8) = [input(6) input(11) input(1)]; 214 input3(9:11) = [input(9) input(14) input(4)]; 215 input3(12:14) = [input(12) input(2) input(7)]; */ 216 { 217 const FIXP_DBL *pSrc = pInput; 218 FIXP_DBL *RESTRICT pDst = aDst; 219 /* Merge 3 loops into one, skip call of fft3 */ 220 for(i=0,l=0,k=0; i<N5; i++, k+=6) 221 { 222 pDst[k+0] = pSrc[l]; 223 pDst[k+1] = pSrc[l+1]; 224 l += 2*N5; 225 if (l >= (2*N15)) 226 l -= (2*N15); 227 228 pDst[k+2] = pSrc[l]; 229 pDst[k+3] = pSrc[l+1]; 230 l += 2*N5; 231 if (l >= (2*N15)) 232 l -= (2*N15); 233 pDst[k+4] = pSrc[l]; 234 pDst[k+5] = pSrc[l+1]; 235 l += (2*N5) + (2*N3); 236 if (l >= (2*N15)) 237 l -= (2*N15); 238 239 /* fft3 merged with shift right by 2 loop */ 240 FIXP_DBL r1,r2,r3; 241 FIXP_DBL s1,s2; 242 /* real part */ 243 r1 = pDst[k+2] + pDst[k+4]; 244 r2 = fMult((pDst[k+2] - pDst[k+4]), C31); 245 s1 = pDst[k+0]; 246 pDst[k+0] = (s1 + r1)>>2; 247 r1 = s1 - (r1>>1); 248 249 /* imaginary part */ 250 s1 = pDst[k+3] + pDst[k+5]; 251 s2 = fMult((pDst[k+3] - pDst[k+5]), C31); 252 r3 = pDst[k+1]; 253 pDst[k+1] = (r3 + s1)>>2; 254 s1 = r3 - (s1>>1); 255 256 /* combination */ 257 pDst[k+2] = (r1 - s2)>>2; 258 pDst[k+4] = (r1 + s2)>>2; 259 pDst[k+3] = (s1 + r2)>>2; 260 pDst[k+5] = (s1 - r2)>>2; 261 } 262 } 263 /* Sort input vector for fft's of length 5 264 input5(0:4) = [output3(0) output3(3) output3(6) output3(9) output3(12)]; 265 input5(5:9) = [output3(1) output3(4) output3(7) output3(10) output3(13)]; 266 input5(10:14) = [output3(2) output3(5) output3(8) output3(11) output3(14)]; */ 267 /* Merge 2 loops into one, brings about 10% */ 268 { 269 const FIXP_DBL *pSrc = aDst; 270 FIXP_DBL *RESTRICT pDst = aDst1; 271 for(i=0,l=0,k=0; i<N3; i++, k+=10) 272 { 273 l = 2*i; 274 pDst[k+0] = pSrc[l+0]; 275 pDst[k+1] = pSrc[l+1]; 276 pDst[k+2] = pSrc[l+0+(2*N3)]; 277 pDst[k+3] = pSrc[l+1+(2*N3)]; 278 pDst[k+4] = pSrc[l+0+(4*N3)]; 279 pDst[k+5] = pSrc[l+1+(4*N3)]; 280 pDst[k+6] = pSrc[l+0+(6*N3)]; 281 pDst[k+7] = pSrc[l+1+(6*N3)]; 282 pDst[k+8] = pSrc[l+0+(8*N3)]; 283 pDst[k+9] = pSrc[l+1+(8*N3)]; 284 fft5(&pDst[k]); 285 } 286 } 287 /* Sort output vector of length 15 288 output = [out5(0) out5(6) out5(12) out5(3) out5(9) 289 out5(10) out5(1) out5(7) out5(13) out5(4) 290 out5(5) out5(11) out5(2) out5(8) out5(14)]; */ 291 /* optimize clumsy loop, brings about 5% */ 292 { 293 const FIXP_DBL *pSrc = aDst1; 294 FIXP_DBL *RESTRICT pDst = pInput; 295 for(i=0,l=0,k=0; i<N3; i++, k+=10) 296 { 297 pDst[k+0] = pSrc[l]; 298 pDst[k+1] = pSrc[l+1]; 299 l += (2*N6); 300 if (l >= (2*N15)) 301 l -= (2*N15); 302 pDst[k+2] = pSrc[l]; 303 pDst[k+3] = pSrc[l+1]; 304 l += (2*N6); 305 if (l >= (2*N15)) 306 l -= (2*N15); 307 pDst[k+4] = pSrc[l]; 308 pDst[k+5] = pSrc[l+1]; 309 l += (2*N6); 310 if (l >= (2*N15)) 311 l -= (2*N15); 312 pDst[k+6] = pSrc[l]; 313 pDst[k+7] = pSrc[l+1]; 314 l += (2*N6); 315 if (l >= (2*N15)) 316 l -= (2*N15); 317 pDst[k+8] = pSrc[l]; 318 pDst[k+9] = pSrc[l+1]; 319 l += 2; /* no modulo check needed, it cannot occur */ 320 } 321 } 322 } 323 324 #define W_PiFOURTH STC(0x5a82799a) 325 #ifndef SUMDIFF_PIFOURTH 326 #define SUMDIFF_PIFOURTH(diff,sum,a,b) \ 327 { \ 328 FIXP_DBL wa, wb;\ 329 wa = fMultDiv2(a, W_PiFOURTH);\ 330 wb = fMultDiv2(b, W_PiFOURTH);\ 331 diff = wb - wa;\ 332 sum = wb + wa;\ 333 } 334 #endif 335 336 /* This version is more overflow save, but less cycle optimal. */ 337 #define SUMDIFF_EIGTH(x, y, ix, iy, vr, vi, ur, ui) \ 338 vr = (x[ 0 + ix]>>1) + (x[16 + ix]>>1); /* Re A + Re B */ \ 339 vi = (x[ 8 + ix]>>1) + (x[24 + ix]>>1); /* Re C + Re D */ \ 340 ur = (x[ 1 + ix]>>1) + (x[17 + ix]>>1); /* Im A + Im B */ \ 341 ui = (x[ 9 + ix]>>1) + (x[25 + ix]>>1); /* Im C + Im D */ \ 342 y[ 0 + iy] = vr + vi; /* Re A' = ReA + ReB +ReC + ReD */ \ 343 y[ 4 + iy] = vr - vi; /* Re C' = -(ReC+ReD) + (ReA+ReB) */ \ 344 y[ 1 + iy] = ur + ui; /* Im A' = sum of imag values */ \ 345 y[ 5 + iy] = ur - ui; /* Im C' = -Im C -Im D +Im A +Im B */ \ 346 vr -= x[16 + ix]; /* Re A - Re B */ \ 347 vi = vi - x[24 + ix]; /* Re C - Re D */ \ 348 ur -= x[17 + ix]; /* Im A - Im B */ \ 349 ui = ui - x[25 + ix]; /* Im C - Im D */ \ 350 y[ 2 + iy] = ui + vr; /* Re B' = Im C - Im D + Re A - Re B */ \ 351 y[ 6 + iy] = vr - ui; /* Re D' = -Im C + Im D + Re A - Re B */ \ 352 y[ 3 + iy] = ur - vi; /* Im B'= -Re C + Re D + Im A - Im B */ \ 353 y[ 7 + iy] = vi + ur; /* Im D'= Re C - Re D + Im A - Im B */ 354 355 static const FIXP_STP fft16_w16[2] = { STCP(0x7641af3d, 0x30fbc54d), STCP(0x30fbc54d, 0x7641af3d) }; 356 357 LNK_SECTION_CODE_L1 358 inline void fft_16(FIXP_DBL *RESTRICT x) 359 { 360 FIXP_DBL vr, vi, ur, ui; 361 FIXP_DBL y[32]; 362 363 SUMDIFF_EIGTH(x, y, 0, 0, vr, vi, ur, ui); 364 SUMDIFF_EIGTH(x, y, 4, 8, vr, vi, ur, ui); 365 SUMDIFF_EIGTH(x, y, 2, 16, vr, vi, ur, ui); 366 SUMDIFF_EIGTH(x, y, 6, 24, vr, vi, ur, ui); 367 368 // xt1 = 0 369 // xt2 = 8 370 vr = y[ 8]; 371 vi = y[ 9]; 372 ur = y[ 0]>>1; 373 ui = y[ 1]>>1; 374 x[ 0] = ur + (vr>>1); 375 x[ 1] = ui + (vi>>1); 376 x[ 8] = ur - (vr>>1); 377 x[ 9] = ui - (vi>>1); 378 379 // xt1 = 4 380 // xt2 = 12 381 vr = y[13]; 382 vi = y[12]; 383 ur = y[ 4]>>1; 384 ui = y[ 5]>>1; 385 x[ 4] = ur + (vr>>1); 386 x[ 5] = ui - (vi>>1); 387 x[12] = ur - (vr>>1); 388 x[13] = ui + (vi>>1); 389 390 // xt1 = 16 391 // xt2 = 24 392 vr = y[24]; 393 vi = y[25]; 394 ur = y[16]>>1; 395 ui = y[17]>>1; 396 x[16] = ur + (vr>>1); 397 x[17] = ui + (vi>>1); 398 x[24] = ur - (vr>>1); 399 x[25] = ui - (vi>>1); 400 401 // xt1 = 20 402 // xt2 = 28 403 vr = y[29]; 404 vi = y[28]; 405 ur = y[20]>>1; 406 ui = y[21]>>1; 407 x[20] = ur + (vr>>1); 408 x[21] = ui - (vi>>1); 409 x[28] = ur - (vr>>1); 410 x[29] = ui + (vi>>1); 411 412 // xt1 = 2 413 // xt2 = 10 414 SUMDIFF_PIFOURTH(vi, vr, y[10], y[11]) 415 //vr = fMultDiv2((y[11] + y[10]),W_PiFOURTH); 416 //vi = fMultDiv2((y[11] - y[10]),W_PiFOURTH); 417 ur = y[ 2]; 418 ui = y[ 3]; 419 x[ 2] = (ur>>1) + vr; 420 x[ 3] = (ui>>1) + vi; 421 x[10] = (ur>>1) - vr; 422 x[11] = (ui>>1) - vi; 423 424 // xt1 = 6 425 // xt2 = 14 426 SUMDIFF_PIFOURTH(vr, vi, y[14], y[15]) 427 ur = y[ 6]; 428 ui = y[ 7]; 429 x[ 6] = (ur>>1) + vr; 430 x[ 7] = (ui>>1) - vi; 431 x[14] = (ur>>1) - vr; 432 x[15] = (ui>>1) + vi; 433 434 // xt1 = 18 435 // xt2 = 26 436 SUMDIFF_PIFOURTH(vi, vr, y[26], y[27]) 437 ur = y[18]; 438 ui = y[19]; 439 x[18] = (ur>>1) + vr; 440 x[19] = (ui>>1) + vi; 441 x[26] = (ur>>1) - vr; 442 x[27] = (ui>>1) - vi; 443 444 // xt1 = 22 445 // xt2 = 30 446 SUMDIFF_PIFOURTH(vr, vi, y[30], y[31]) 447 ur = y[22]; 448 ui = y[23]; 449 x[22] = (ur>>1) + vr; 450 x[23] = (ui>>1) - vi; 451 x[30] = (ur>>1) - vr; 452 x[31] = (ui>>1) + vi; 453 454 // xt1 = 0 455 // xt2 = 16 456 vr = x[16]; 457 vi = x[17]; 458 ur = x[ 0]>>1; 459 ui = x[ 1]>>1; 460 x[ 0] = ur + (vr>>1); 461 x[ 1] = ui + (vi>>1); 462 x[16] = ur - (vr>>1); 463 x[17] = ui - (vi>>1); 464 465 // xt1 = 8 466 // xt2 = 24 467 vi = x[24]; 468 vr = x[25]; 469 ur = x[ 8]>>1; 470 ui = x[ 9]>>1; 471 x[ 8] = ur + (vr>>1); 472 x[ 9] = ui - (vi>>1); 473 x[24] = ur - (vr>>1); 474 x[25] = ui + (vi>>1); 475 476 // xt1 = 2 477 // xt2 = 18 478 cplxMultDiv2(&vi, &vr, x[19], x[18], fft16_w16[0]); 479 ur = x[ 2]; 480 ui = x[ 3]; 481 x[ 2] = (ur>>1) + vr; 482 x[ 3] = (ui>>1) + vi; 483 x[18] = (ur>>1) - vr; 484 x[19] = (ui>>1) - vi; 485 486 // xt1 = 10 487 // xt2 = 26 488 cplxMultDiv2(&vr, &vi, x[27], x[26], fft16_w16[0]); 489 ur = x[10]; 490 ui = x[11]; 491 x[10] = (ur>>1) + vr; 492 x[11] = (ui>>1) - vi; 493 x[26] = (ur>>1) - vr; 494 x[27] = (ui>>1) + vi; 495 496 // xt1 = 4 497 // xt2 = 20 498 SUMDIFF_PIFOURTH(vi, vr, x[20], x[21]) 499 ur = x[ 4]; 500 ui = x[ 5]; 501 x[ 4] = (ur>>1) + vr; 502 x[ 5] = (ui>>1) + vi; 503 x[20] = (ur>>1) - vr; 504 x[21] = (ui>>1) - vi; 505 506 // xt1 = 12 507 // xt2 = 28 508 SUMDIFF_PIFOURTH(vr, vi, x[28], x[29]) 509 ur = x[12]; 510 ui = x[13]; 511 x[12] = (ur>>1) + vr; 512 x[13] = (ui>>1) - vi; 513 x[28] = (ur>>1) - vr; 514 x[29] = (ui>>1) + vi; 515 516 // xt1 = 6 517 // xt2 = 22 518 cplxMultDiv2(&vi, &vr, x[23], x[22], fft16_w16[1]); 519 ur = x[ 6]; 520 ui = x[ 7]; 521 x[ 6] = (ur>>1) + vr; 522 x[ 7] = (ui>>1) + vi; 523 x[22] = (ur>>1) - vr; 524 x[23] = (ui>>1) - vi; 525 526 // xt1 = 14 527 // xt2 = 30 528 cplxMultDiv2(&vr, &vi, x[31], x[30], fft16_w16[1]); 529 ur = x[14]; 530 ui = x[15]; 531 x[14] = (ur>>1) + vr; 532 x[15] = (ui>>1) - vi; 533 x[30] = (ur>>1) - vr; 534 x[31] = (ui>>1) + vi; 535 } 536 537 #ifndef FUNCTION_fft_32 538 static const FIXP_STP fft32_w32[6] = 539 { 540 STCP (0x7641af3d, 0x30fbc54d), STCP(0x30fbc54d, 0x7641af3d), STCP(0x7d8a5f40, 0x18f8b83c), 541 STCP (0x6a6d98a4, 0x471cece7), STCP(0x471cece7, 0x6a6d98a4), STCP(0x18f8b83c, 0x7d8a5f40) 542 }; 543 544 LNK_SECTION_CODE_L1 545 inline void fft_32(FIXP_DBL *x) 546 { 547 548 #define W_PiFOURTH STC(0x5a82799a) 549 550 FIXP_DBL vr,vi,ur,ui; 551 FIXP_DBL y[64]; 552 553 /* 554 * 1+2 stage radix 4 555 */ 556 557 ///////////////////////////////////////////////////////////////////////////////////////// 558 559 // i = 0 560 vr = (x[ 0] + x[32])>>1; /* Re A + Re B */ 561 vi = (x[16] + x[48]); /* Re C + Re D */ 562 ur = (x[ 1] + x[33])>>1; /* Im A + Im B */ 563 ui = (x[17] + x[49]); /* Im C + Im D */ 564 565 y[ 0] = vr + (vi>>1); /* Re A' = ReA + ReB +ReC + ReD */ 566 y[ 4] = vr - (vi>>1); /* Re C' = -(ReC+ReD) + (ReA+ReB) */ 567 y[ 1] = ur + (ui>>1); /* Im A' = sum of imag values */ 568 y[ 5] = ur - (ui>>1); /* Im C' = -Im C -Im D +Im A +Im B */ 569 570 vr -= x[32]; /* Re A - Re B */ 571 vi = (vi>>1) - x[48]; /* Re C - Re D */ 572 ur -= x[33]; /* Im A - Im B */ 573 ui = (ui>>1) - x[49]; /* Im C - Im D */ 574 575 y[ 2] = ui + vr; /* Re B' = Im C - Im D + Re A - Re B */ 576 y[ 6] = vr - ui; /* Re D' = -Im C + Im D + Re A - Re B */ 577 y[ 3] = ur - vi; /* Im B'= -Re C + Re D + Im A - Im B */ 578 y[ 7] = vi + ur; /* Im D'= Re C - Re D + Im A - Im B */ 579 580 //i=8 581 vr = (x[ 8] + x[40])>>1; /* Re A + Re B */ 582 vi = (x[24] + x[56]); /* Re C + Re D */ 583 ur = (x[ 9] + x[41])>>1; /* Im A + Im B */ 584 ui = (x[25] + x[57]); /* Im C + Im D */ 585 586 y[ 8] = vr + (vi>>1); /* Re A' = ReA + ReB +ReC + ReD */ 587 y[12] = vr - (vi>>1); /* Re C' = -(ReC+ReD) + (ReA+ReB) */ 588 y[ 9] = ur + (ui>>1); /* Im A' = sum of imag values */ 589 y[13] = ur - (ui>>1); /* Im C' = -Im C -Im D +Im A +Im B */ 590 591 vr -= x[40]; /* Re A - Re B */ 592 vi = (vi>>1) - x[56]; /* Re C - Re D */ 593 ur -= x[41]; /* Im A - Im B */ 594 ui = (ui>>1) - x[57]; /* Im C - Im D */ 595 596 y[10] = ui + vr; /* Re B' = Im C - Im D + Re A - Re B */ 597 y[14] = vr - ui; /* Re D' = -Im C + Im D + Re A - Re B */ 598 y[11] = ur - vi; /* Im B'= -Re C + Re D + Im A - Im B */ 599 y[15] = vi + ur; /* Im D'= Re C - Re D + Im A - Im B */ 600 601 //i=16 602 vr = (x[ 4] + x[36])>>1; /* Re A + Re B */ 603 vi = (x[20] + x[52]); /* Re C + Re D */ 604 ur = (x[ 5] + x[37])>>1; /* Im A + Im B */ 605 ui = (x[21] + x[53]); /* Im C + Im D */ 606 607 y[16] = vr + (vi>>1); /* Re A' = ReA + ReB +ReC + ReD */ 608 y[20] = vr - (vi>>1); /* Re C' = -(ReC+ReD) + (ReA+ReB) */ 609 y[17] = ur + (ui>>1); /* Im A' = sum of imag values */ 610 y[21] = ur - (ui>>1); /* Im C' = -Im C -Im D +Im A +Im B */ 611 612 vr -= x[36]; /* Re A - Re B */ 613 vi = (vi>>1) - x[52]; /* Re C - Re D */ 614 ur -= x[37]; /* Im A - Im B */ 615 ui = (ui>>1) - x[53]; /* Im C - Im D */ 616 617 y[18] = ui + vr; /* Re B' = Im C - Im D + Re A - Re B */ 618 y[22] = vr - ui; /* Re D' = -Im C + Im D + Re A - Re B */ 619 y[19] = ur - vi; /* Im B'= -Re C + Re D + Im A - Im B */ 620 y[23] = vi + ur; /* Im D'= Re C - Re D + Im A - Im B */ 621 622 //i=24 623 vr = (x[12] + x[44])>>1; /* Re A + Re B */ 624 vi = (x[28] + x[60]); /* Re C + Re D */ 625 ur = (x[13] + x[45])>>1; /* Im A + Im B */ 626 ui = (x[29] + x[61]); /* Im C + Im D */ 627 628 y[24] = vr + (vi>>1); /* Re A' = ReA + ReB +ReC + ReD */ 629 y[28] = vr - (vi>>1); /* Re C' = -(ReC+ReD) + (ReA+ReB) */ 630 y[25] = ur + (ui>>1); /* Im A' = sum of imag values */ 631 y[29] = ur - (ui>>1); /* Im C' = -Im C -Im D +Im A +Im B */ 632 633 vr -= x[44]; /* Re A - Re B */ 634 vi = (vi>>1) - x[60]; /* Re C - Re D */ 635 ur -= x[45]; /* Im A - Im B */ 636 ui = (ui>>1) - x[61]; /* Im C - Im D */ 637 638 y[26] = ui + vr; /* Re B' = Im C - Im D + Re A - Re B */ 639 y[30] = vr - ui; /* Re D' = -Im C + Im D + Re A - Re B */ 640 y[27] = ur - vi; /* Im B'= -Re C + Re D + Im A - Im B */ 641 y[31] = vi + ur; /* Im D'= Re C - Re D + Im A - Im B */ 642 643 // i = 32 644 vr = (x[ 2] + x[34])>>1; /* Re A + Re B */ 645 vi = (x[18] + x[50]); /* Re C + Re D */ 646 ur = (x[ 3] + x[35])>>1; /* Im A + Im B */ 647 ui = (x[19] + x[51]); /* Im C + Im D */ 648 649 y[32] = vr + (vi>>1); /* Re A' = ReA + ReB +ReC + ReD */ 650 y[36] = vr - (vi>>1); /* Re C' = -(ReC+ReD) + (ReA+ReB) */ 651 y[33] = ur + (ui>>1); /* Im A' = sum of imag values */ 652 y[37] = ur - (ui>>1); /* Im C' = -Im C -Im D +Im A +Im B */ 653 654 vr -= x[34]; /* Re A - Re B */ 655 vi = (vi>>1) - x[50]; /* Re C - Re D */ 656 ur -= x[35]; /* Im A - Im B */ 657 ui = (ui>>1) - x[51]; /* Im C - Im D */ 658 659 y[34] = ui + vr; /* Re B' = Im C - Im D + Re A - Re B */ 660 y[38] = vr - ui; /* Re D' = -Im C + Im D + Re A - Re B */ 661 y[35] = ur - vi; /* Im B'= -Re C + Re D + Im A - Im B */ 662 y[39] = vi + ur; /* Im D'= Re C - Re D + Im A - Im B */ 663 664 //i=40 665 vr = (x[10] + x[42])>>1; /* Re A + Re B */ 666 vi = (x[26] + x[58]); /* Re C + Re D */ 667 ur = (x[11] + x[43])>>1; /* Im A + Im B */ 668 ui = (x[27] + x[59]); /* Im C + Im D */ 669 670 y[40] = vr + (vi>>1); /* Re A' = ReA + ReB +ReC + ReD */ 671 y[44] = vr - (vi>>1); /* Re C' = -(ReC+ReD) + (ReA+ReB) */ 672 y[41] = ur + (ui>>1); /* Im A' = sum of imag values */ 673 y[45] = ur - (ui>>1); /* Im C' = -Im C -Im D +Im A +Im B */ 674 675 vr -= x[42]; /* Re A - Re B */ 676 vi = (vi>>1) - x[58]; /* Re C - Re D */ 677 ur -= x[43]; /* Im A - Im B */ 678 ui = (ui>>1) - x[59]; /* Im C - Im D */ 679 680 y[42] = ui + vr; /* Re B' = Im C - Im D + Re A - Re B */ 681 y[46] = vr - ui; /* Re D' = -Im C + Im D + Re A - Re B */ 682 y[43] = ur - vi; /* Im B'= -Re C + Re D + Im A - Im B */ 683 y[47] = vi + ur; /* Im D'= Re C - Re D + Im A - Im B */ 684 685 //i=48 686 vr = (x[ 6] + x[38])>>1; /* Re A + Re B */ 687 vi = (x[22] + x[54]); /* Re C + Re D */ 688 ur = (x[ 7] + x[39])>>1; /* Im A + Im B */ 689 ui = (x[23] + x[55]); /* Im C + Im D */ 690 691 y[48] = vr + (vi>>1); /* Re A' = ReA + ReB +ReC + ReD */ 692 y[52] = vr - (vi>>1); /* Re C' = -(ReC+ReD) + (ReA+ReB) */ 693 y[49] = ur + (ui>>1); /* Im A' = sum of imag values */ 694 y[53] = ur - (ui>>1); /* Im C' = -Im C -Im D +Im A +Im B */ 695 696 vr -= x[38]; /* Re A - Re B */ 697 vi = (vi>>1) - x[54]; /* Re C - Re D */ 698 ur -= x[39]; /* Im A - Im B */ 699 ui = (ui>>1) - x[55]; /* Im C - Im D */ 700 701 y[50] = ui + vr; /* Re B' = Im C - Im D + Re A - Re B */ 702 y[54] = vr - ui; /* Re D' = -Im C + Im D + Re A - Re B */ 703 y[51] = ur - vi; /* Im B'= -Re C + Re D + Im A - Im B */ 704 y[55] = vi + ur; /* Im D'= Re C - Re D + Im A - Im B */ 705 706 //i=56 707 vr = (x[14] + x[46])>>1; /* Re A + Re B */ 708 vi = (x[30] + x[62]); /* Re C + Re D */ 709 ur = (x[15] + x[47])>>1; /* Im A + Im B */ 710 ui = (x[31] + x[63]); /* Im C + Im D */ 711 712 y[56] = vr + (vi>>1); /* Re A' = ReA + ReB +ReC + ReD */ 713 y[60] = vr - (vi>>1); /* Re C' = -(ReC+ReD) + (ReA+ReB) */ 714 y[57] = ur + (ui>>1); /* Im A' = sum of imag values */ 715 y[61] = ur - (ui>>1); /* Im C' = -Im C -Im D +Im A +Im B */ 716 717 vr -= x[46]; /* Re A - Re B */ 718 vi = (vi>>1) - x[62]; /* Re C - Re D */ 719 ur -= x[47]; /* Im A - Im B */ 720 ui = (ui>>1) - x[63]; /* Im C - Im D */ 721 722 y[58] = ui + vr; /* Re B' = Im C - Im D + Re A - Re B */ 723 y[62] = vr - ui; /* Re D' = -Im C + Im D + Re A - Re B */ 724 y[59] = ur - vi; /* Im B'= -Re C + Re D + Im A - Im B */ 725 y[63] = vi + ur; /* Im D'= Re C - Re D + Im A - Im B */ 726 727 728 FIXP_DBL *xt = &x[0]; 729 FIXP_DBL *yt = &y[0]; 730 731 int j = 4; 732 do 733 { 734 vr = yt[8]; 735 vi = yt[9]; 736 ur = yt[0]>>1; 737 ui = yt[1]>>1; 738 xt[ 0] = ur + (vr>>1); 739 xt[ 1] = ui + (vi>>1); 740 xt[ 8] = ur - (vr>>1); 741 xt[ 9] = ui - (vi>>1); 742 743 vr = yt[13]; 744 vi = yt[12]; 745 ur = yt[4]>>1; 746 ui = yt[5]>>1; 747 xt[ 4] = ur + (vr>>1); 748 xt[ 5] = ui - (vi>>1); 749 xt[12] = ur - (vr>>1); 750 xt[13] = ui + (vi>>1); 751 752 SUMDIFF_PIFOURTH(vi, vr, yt[10], yt[11]) 753 ur = yt[2]; 754 ui = yt[3]; 755 xt[ 2] = (ur>>1) + vr; 756 xt[ 3] = (ui>>1) + vi; 757 xt[10] = (ur>>1) - vr; 758 xt[11] = (ui>>1) - vi; 759 760 SUMDIFF_PIFOURTH(vr, vi, yt[14], yt[15]) 761 ur = yt[6]; 762 ui = yt[7]; 763 764 xt[ 6] = (ur>>1) + vr; 765 xt[ 7] = (ui>>1) - vi; 766 xt[14] = (ur>>1) - vr; 767 xt[15] = (ui>>1) + vi; 768 xt += 16; 769 yt += 16; 770 } while (--j != 0); 771 772 vr = x[16]; 773 vi = x[17]; 774 ur = x[ 0]>>1; 775 ui = x[ 1]>>1; 776 x[ 0] = ur + (vr>>1); 777 x[ 1] = ui + (vi>>1); 778 x[16] = ur - (vr>>1); 779 x[17] = ui - (vi>>1); 780 781 vi = x[24]; 782 vr = x[25]; 783 ur = x[ 8]>>1; 784 ui = x[ 9]>>1; 785 x[ 8] = ur + (vr>>1); 786 x[ 9] = ui - (vi>>1); 787 x[24] = ur - (vr>>1); 788 x[25] = ui + (vi>>1); 789 790 vr = x[48]; 791 vi = x[49]; 792 ur = x[32]>>1; 793 ui = x[33]>>1; 794 x[32] = ur + (vr>>1); 795 x[33] = ui + (vi>>1); 796 x[48] = ur - (vr>>1); 797 x[49] = ui - (vi>>1); 798 799 vi = x[56]; 800 vr = x[57]; 801 ur = x[40]>>1; 802 ui = x[41]>>1; 803 x[40] = ur + (vr>>1); 804 x[41] = ui - (vi>>1); 805 x[56] = ur - (vr>>1); 806 x[57] = ui + (vi>>1); 807 808 cplxMultDiv2(&vi, &vr, x[19], x[18], fft32_w32[0]); 809 ur = x[ 2]; 810 ui = x[ 3]; 811 x[ 2] = (ur>>1) + vr; 812 x[ 3] = (ui>>1) + vi; 813 x[18] = (ur>>1) - vr; 814 x[19] = (ui>>1) - vi; 815 816 cplxMultDiv2(&vr, &vi, x[27], x[26], fft32_w32[0]); 817 ur = x[10]; 818 ui = x[11]; 819 x[10] = (ur>>1) + vr; 820 x[11] = (ui>>1) - vi; 821 x[26] = (ur>>1) - vr; 822 x[27] = (ui>>1) + vi; 823 824 cplxMultDiv2(&vi, &vr, x[51], x[50], fft32_w32[0]); 825 ur = x[34]; 826 ui = x[35]; 827 x[34] = (ur>>1) + vr; 828 x[35] = (ui>>1) + vi; 829 x[50] = (ur>>1) - vr; 830 x[51] = (ui>>1) - vi; 831 832 cplxMultDiv2(&vr, &vi, x[59], x[58], fft32_w32[0]); 833 ur = x[42]; 834 ui = x[43]; 835 x[42] = (ur>>1) + vr; 836 x[43] = (ui>>1) - vi; 837 x[58] = (ur>>1) - vr; 838 x[59] = (ui>>1) + vi; 839 840 SUMDIFF_PIFOURTH(vi, vr, x[20], x[21]) 841 ur = x[ 4]; 842 ui = x[ 5]; 843 x[ 4] = (ur>>1) + vr; 844 x[ 5] = (ui>>1) + vi; 845 x[20] = (ur>>1) - vr; 846 x[21] = (ui>>1) - vi; 847 848 SUMDIFF_PIFOURTH(vr, vi, x[28], x[29]) 849 ur = x[12]; 850 ui = x[13]; 851 x[12] = (ur>>1) + vr; 852 x[13] = (ui>>1) - vi; 853 x[28] = (ur>>1) - vr; 854 x[29] = (ui>>1) + vi; 855 856 SUMDIFF_PIFOURTH(vi, vr, x[52], x[53]) 857 ur = x[36]; 858 ui = x[37]; 859 x[36] = (ur>>1) + vr; 860 x[37] = (ui>>1) + vi; 861 x[52] = (ur>>1) - vr; 862 x[53] = (ui>>1) - vi; 863 864 SUMDIFF_PIFOURTH(vr, vi, x[60], x[61]) 865 ur = x[44]; 866 ui = x[45]; 867 x[44] = (ur>>1) + vr; 868 x[45] = (ui>>1) - vi; 869 x[60] = (ur>>1) - vr; 870 x[61] = (ui>>1) + vi; 871 872 873 cplxMultDiv2(&vi, &vr, x[23], x[22], fft32_w32[1]); 874 ur = x[ 6]; 875 ui = x[ 7]; 876 x[ 6] = (ur>>1) + vr; 877 x[ 7] = (ui>>1) + vi; 878 x[22] = (ur>>1) - vr; 879 x[23] = (ui>>1) - vi; 880 881 cplxMultDiv2(&vr, &vi, x[31], x[30], fft32_w32[1]); 882 ur = x[14]; 883 ui = x[15]; 884 x[14] = (ur>>1) + vr; 885 x[15] = (ui>>1) - vi; 886 x[30] = (ur>>1) - vr; 887 x[31] = (ui>>1) + vi; 888 889 cplxMultDiv2(&vi, &vr, x[55], x[54], fft32_w32[1]); 890 ur = x[38]; 891 ui = x[39]; 892 x[38] = (ur>>1) + vr; 893 x[39] = (ui>>1) + vi; 894 x[54] = (ur>>1) - vr; 895 x[55] = (ui>>1) - vi; 896 897 cplxMultDiv2(&vr, &vi, x[63], x[62], fft32_w32[1]); 898 ur = x[46]; 899 ui = x[47]; 900 901 x[46] = (ur>>1) + vr; 902 x[47] = (ui>>1) - vi; 903 x[62] = (ur>>1) - vr; 904 x[63] = (ui>>1) + vi; 905 906 vr = x[32]; 907 vi = x[33]; 908 ur = x[ 0]>>1; 909 ui = x[ 1]>>1; 910 x[ 0] = ur + (vr>>1); 911 x[ 1] = ui + (vi>>1); 912 x[32] = ur - (vr>>1); 913 x[33] = ui - (vi>>1); 914 915 vi = x[48]; 916 vr = x[49]; 917 ur = x[16]>>1; 918 ui = x[17]>>1; 919 x[16] = ur + (vr>>1); 920 x[17] = ui - (vi>>1); 921 x[48] = ur - (vr>>1); 922 x[49] = ui + (vi>>1); 923 924 cplxMultDiv2(&vi, &vr, x[35], x[34], fft32_w32[2]); 925 ur = x[ 2]; 926 ui = x[ 3]; 927 x[ 2] = (ur>>1) + vr; 928 x[ 3] = (ui>>1) + vi; 929 x[34] = (ur>>1) - vr; 930 x[35] = (ui>>1) - vi; 931 932 cplxMultDiv2(&vr, &vi, x[51], x[50], fft32_w32[2]); 933 ur = x[18]; 934 ui = x[19]; 935 x[18] = (ur>>1) + vr; 936 x[19] = (ui>>1) - vi; 937 x[50] = (ur>>1) - vr; 938 x[51] = (ui>>1) + vi; 939 940 cplxMultDiv2(&vi, &vr, x[37], x[36], fft32_w32[0]); 941 ur = x[ 4]; 942 ui = x[ 5]; 943 x[ 4] = (ur>>1) + vr; 944 x[ 5] = (ui>>1) + vi; 945 x[36] = (ur>>1) - vr; 946 x[37] = (ui>>1) - vi; 947 948 cplxMultDiv2(&vr, &vi, x[53], x[52], fft32_w32[0]); 949 ur = x[20]; 950 ui = x[21]; 951 x[20] = (ur>>1) + vr; 952 x[21] = (ui>>1) - vi; 953 x[52] = (ur>>1) - vr; 954 x[53] = (ui>>1) + vi; 955 956 cplxMultDiv2(&vi, &vr, x[39], x[38], fft32_w32[3]); 957 ur = x[ 6]; 958 ui = x[ 7]; 959 x[ 6] = (ur>>1) + vr; 960 x[ 7] = (ui>>1) + vi; 961 x[38] = (ur>>1) - vr; 962 x[39] = (ui>>1) - vi; 963 964 cplxMultDiv2(&vr, &vi, x[55], x[54], fft32_w32[3]); 965 ur = x[22]; 966 ui = x[23]; 967 x[22] = (ur>>1) + vr; 968 x[23] = (ui>>1) - vi; 969 x[54] = (ur>>1) - vr; 970 x[55] = (ui>>1) + vi; 971 972 SUMDIFF_PIFOURTH(vi, vr, x[40], x[41]) 973 ur = x[ 8]; 974 ui = x[ 9]; 975 x[ 8] = (ur>>1) + vr; 976 x[ 9] = (ui>>1) + vi; 977 x[40] = (ur>>1) - vr; 978 x[41] = (ui>>1) - vi; 979 980 SUMDIFF_PIFOURTH(vr, vi, x[56], x[57]) 981 ur = x[24]; 982 ui = x[25]; 983 x[24] = (ur>>1) + vr; 984 x[25] = (ui>>1) - vi; 985 x[56] = (ur>>1) - vr; 986 x[57] = (ui>>1) + vi; 987 988 cplxMultDiv2(&vi, &vr, x[43], x[42], fft32_w32[4]); 989 ur = x[10]; 990 ui = x[11]; 991 992 x[10] = (ur>>1) + vr; 993 x[11] = (ui>>1) + vi; 994 x[42] = (ur>>1) - vr; 995 x[43] = (ui>>1) - vi; 996 997 cplxMultDiv2(&vr, &vi, x[59], x[58], fft32_w32[4]); 998 ur = x[26]; 999 ui = x[27]; 1000 x[26] = (ur>>1) + vr; 1001 x[27] = (ui>>1) - vi; 1002 x[58] = (ur>>1) - vr; 1003 x[59] = (ui>>1) + vi; 1004 1005 cplxMultDiv2(&vi, &vr, x[45], x[44], fft32_w32[1]); 1006 ur = x[12]; 1007 ui = x[13]; 1008 x[12] = (ur>>1) + vr; 1009 x[13] = (ui>>1) + vi; 1010 x[44] = (ur>>1) - vr; 1011 x[45] = (ui>>1) - vi; 1012 1013 cplxMultDiv2(&vr, &vi, x[61], x[60], fft32_w32[1]); 1014 ur = x[28]; 1015 ui = x[29]; 1016 x[28] = (ur>>1) + vr; 1017 x[29] = (ui>>1) - vi; 1018 x[60] = (ur>>1) - vr; 1019 x[61] = (ui>>1) + vi; 1020 1021 cplxMultDiv2(&vi, &vr, x[47], x[46], fft32_w32[5]); 1022 ur = x[14]; 1023 ui = x[15]; 1024 x[14] = (ur>>1) + vr; 1025 x[15] = (ui>>1) + vi; 1026 x[46] = (ur>>1) - vr; 1027 x[47] = (ui>>1) - vi; 1028 1029 cplxMultDiv2(&vr, &vi, x[63], x[62], fft32_w32[5]); 1030 ur = x[30]; 1031 ui = x[31]; 1032 x[30] = (ur>>1) + vr; 1033 x[31] = (ui>>1) - vi; 1034 x[62] = (ur>>1) - vr; 1035 x[63] = (ui>>1) + vi; 1036 } 1037 #endif /* #ifndef FUNCTION_fft_32 */ 1038 1039 1040 /** 1041 * \brief Apply rotation vectors to a data buffer. 1042 * \param cl length of each row of input data. 1043 * \param l total length of input data. 1044 * \param pVecRe real part of rotation ceofficient vector. 1045 * \param pVecIm imaginary part of rotation ceofficient vector. 1046 */ 1047 static inline void fft_apply_rot_vector(FIXP_DBL *RESTRICT pData, const int cl, const int l, const FIXP_STB *pVecRe, const FIXP_STB *pVecIm) 1048 { 1049 FIXP_DBL re, im; 1050 FIXP_STB vre, vim; 1051 1052 int i, c; 1053 1054 for(i=0; i<cl; i++) { 1055 re = pData[2*i]; 1056 im = pData[2*i+1]; 1057 1058 pData[2*i] = re>>2; /* * 0.25 */ 1059 pData[2*i+1] = im>>2; /* * 0.25 */ 1060 } 1061 for(; i<l; i+=cl) 1062 { 1063 re = pData[2*i]; 1064 im = pData[2*i+1]; 1065 1066 pData[2*i] = re>>2; /* * 0.25 */ 1067 pData[2*i+1] = im>>2; /* * 0.25 */ 1068 1069 for (c=i+1; c<i+cl; c++) 1070 { 1071 re = pData[2*c]>>1; 1072 im = pData[2*c+1]>>1; 1073 vre = *pVecRe++; 1074 vim = *pVecIm++; 1075 1076 cplxMultDiv2(&pData[2*c+1], &pData[2*c], im, re, vre, vim); 1077 } 1078 } 1079 } 1080 1081 #define FFT_TWO_STAGE_MACRO_ENABLE 1082 1083 1084 #ifdef FFT_TWO_STAGE_MACRO_ENABLE 1085 1086 #define fftN2(pInput, length, dim1, dim2, fft_func1, fft_func2, RotVectorReal, RotVectorImag) \ 1087 { \ 1088 int i, j; \ 1089 \ 1090 C_ALLOC_SCRATCH_START(aDst, FIXP_DBL, length*2); \ 1091 C_ALLOC_SCRATCH_START(aDst2, FIXP_DBL, dim2*2); \ 1092 \ 1093 FDK_ASSERT(length == dim1*dim2); \ 1094 \ 1095 /* Perform dim2 times the fft of length dim1. The input samples are at the address of pSrc and the \ 1096 output samples are at the address of pDst. The input vector for the fft of length dim1 is built \ 1097 of the interleaved samples in pSrc, the output samples are stored consecutively. \ 1098 */ \ 1099 { \ 1100 const FIXP_DBL* pSrc = pInput; \ 1101 FIXP_DBL *RESTRICT pDst = aDst; \ 1102 \ 1103 for(i=0; i<dim2; i++) \ 1104 { \ 1105 for(j=0; j<dim1; j++) \ 1106 { \ 1107 pDst[2*j] = pSrc[2*j*dim2]; \ 1108 pDst[2*j+1] = pSrc[2*j*dim2+1]; \ 1109 } \ 1110 \ 1111 fft_func1(pDst); \ 1112 pSrc += 2; \ 1113 pDst = pDst + 2*dim1; \ 1114 } \ 1115 } \ 1116 \ 1117 /* Perform the modulation of the output of the fft of length dim1 */ \ 1118 fft_apply_rot_vector(aDst, dim1, length, RotVectorReal, RotVectorImag); \ 1119 \ 1120 /* Perform dim1 times the fft of length dim2. The input samples are at the address of aDst and the \ 1121 output samples are at the address of pInput. The input vector for the fft of length dim2 is built \ 1122 of the interleaved samples in aDst, the output samples are stored consecutively at the address \ 1123 of pInput. \ 1124 */ \ 1125 { \ 1126 const FIXP_DBL* pSrc = aDst; \ 1127 FIXP_DBL *RESTRICT pDst = aDst2; \ 1128 FIXP_DBL *RESTRICT pDstOut = pInput; \ 1129 \ 1130 for(i=0; i<dim1; i++) \ 1131 { \ 1132 for(j=0; j<dim2; j++) \ 1133 { \ 1134 pDst[2*j] = pSrc[2*j*dim1]; \ 1135 pDst[2*j+1] = pSrc[2*j*dim1+1]; \ 1136 } \ 1137 \ 1138 fft_func2(pDst); \ 1139 \ 1140 for(j=0; j<dim2; j++) \ 1141 { \ 1142 pDstOut[2*j*dim1] = pDst[2*j]; \ 1143 pDstOut[2*j*dim1+1] = pDst[2*j+1]; \ 1144 } \ 1145 pSrc += 2; \ 1146 pDstOut += 2; \ 1147 } \ 1148 } \ 1149 \ 1150 C_ALLOC_SCRATCH_END(aDst2, FIXP_DBL, dim2*2); \ 1151 C_ALLOC_SCRATCH_END(aDst, FIXP_DBL, length*2); \ 1152 } \ 1153 1154 #else /* FFT_TWO_STAGE_MACRO_ENABLE */ 1155 1156 /* select either switch case of function pointer. */ 1157 //#define FFT_TWO_STAGE_SWITCH_CASE 1158 1159 static inline void fftN2( 1160 FIXP_DBL *pInput, 1161 const int length, 1162 const int dim1, 1163 const int dim2, 1164 void (* const fft1)(FIXP_DBL *), 1165 void (* const fft2)(FIXP_DBL *), 1166 const FIXP_STB *RotVectorReal, 1167 const FIXP_STB *RotVectorImag 1168 ) 1169 { 1170 /* The real part of the input samples are at the addresses with even indices and the imaginary 1171 part of the input samples are at the addresses with odd indices. The output samples are stored 1172 at the address of pInput 1173 */ 1174 FIXP_DBL *pSrc, *pDst, *pDstOut; 1175 int i, j; 1176 1177 C_ALLOC_SCRATCH_START(aDst, FIXP_DBL, length*2); 1178 C_ALLOC_SCRATCH_START(aDst2, FIXP_DBL, dim2*2); 1179 1180 FDK_ASSERT(length == dim1*dim2); 1181 1182 /* Perform dim2 times the fft of length dim1. The input samples are at the address of pSrc and the 1183 output samples are at the address of pDst. The input vector for the fft of length dim1 is built 1184 of the interleaved samples in pSrc, the output samples are stored consecutively. 1185 */ 1186 pSrc = pInput; 1187 pDst = aDst; 1188 for(i=0; i<length/dim1; i++) 1189 { 1190 for(j=0; j<length/dim2; j++) 1191 { 1192 pDst[2*j] = pSrc[2*j*dim2]; 1193 pDst[2*j+1] = pSrc[2*j*dim2+1]; 1194 } 1195 1196 /* fft of size dim1 */ 1197 #ifndef FFT_TWO_STAGE_SWITCH_CASE 1198 fft1(pDst); 1199 #else 1200 switch (dim1) { 1201 case 3: fft3(pDst); break; 1202 case 4: fft_4(pDst); break; 1203 case 5: fft5(pDst); break; 1204 case 8: fft_8(pDst); break; 1205 case 15: fft15(pDst); break; 1206 case 16: fft_16(pDst); break; 1207 case 32: fft_32(pDst); break; 1208 /*case 64: fft_64(pDst); break;*/ 1209 case 128: fft_128(pDst); break; 1210 } 1211 #endif 1212 pSrc += 2; 1213 pDst = pDst + 2*length/dim2; 1214 } 1215 1216 /* Perform the modulation of the output of the fft of length dim1 */ 1217 pSrc=aDst; 1218 fft_apply_rot_vector(pSrc, length/dim2, length, RotVectorReal, RotVectorImag); 1219 1220 /* Perform dim1 times the fft of length dim2. The input samples are at the address of aDst and the 1221 output samples are at the address of pInput. The input vector for the fft of length dim2 is built 1222 of the interleaved samples in aDst, the output samples are stored consecutively at the address 1223 of pInput. 1224 */ 1225 pSrc = aDst; 1226 pDst = aDst2; 1227 pDstOut = pInput; 1228 for(i=0; i<length/dim2; i++) 1229 { 1230 for(j=0; j<length/dim1; j++) 1231 { 1232 pDst[2*j] = pSrc[2*j*dim1]; 1233 pDst[2*j+1] = pSrc[2*j*dim1+1]; 1234 } 1235 1236 #ifndef FFT_TWO_STAGE_SWITCH_CASE 1237 fft2(pDst); 1238 #else 1239 switch (dim2) { 1240 case 3: fft3(pDst); break; 1241 case 4: fft_4(pDst); break; 1242 case 5: fft5(pDst); break; 1243 case 8: fft_8(pDst); break; 1244 case 15: fft15(pDst); break; 1245 case 16: fft_16(pDst); break; 1246 case 32: fft_32(pDst); break; 1247 /*case 64: fft_64(pDst); break;*/ 1248 case 128: fft_128(pDst); break; 1249 } 1250 #endif 1251 1252 for(j=0; j<length/dim1; j++) 1253 { 1254 pDstOut[2*j*dim1] = pDst[2*j]; 1255 pDstOut[2*j*dim1+1] = pDst[2*j+1]; 1256 } 1257 pSrc += 2; 1258 pDstOut += 2; 1259 } 1260 1261 C_ALLOC_SCRATCH_END(aDst2, FIXP_DBL, dim2*2); 1262 C_ALLOC_SCRATCH_END(aDst, FIXP_DBL, length*2); 1263 } 1264 1265 #endif /* FFT_TWO_STAGE_MACRO_ENABLE */ 1266 1267 1268 1269 1270 1271 1272 1273 1274 1275 1276 1277 1278 #define SCALEFACTOR60 5 1279 /** 1280 The function performs the fft of length 60. It is splittet into fft's of length 4 and fft's of 1281 length 15. Between the fft's a modolation is calculated. 1282 */ 1283 static inline void fft60(FIXP_DBL *pInput, INT *pScalefactor) 1284 { 1285 fftN2( 1286 pInput, 60, 4, 15, 1287 fft_4, fft15, 1288 RotVectorReal60, RotVectorImag60 1289 ); 1290 *pScalefactor += SCALEFACTOR60; 1291 } 1292 1293 1294 1295 /* Fallback implementation in case of no better implementation available. */ 1296 1297 #define SCALEFACTOR240 7 1298 1299 /** 1300 The function performs the fft of length 240. It is splittet into fft's of length 16 and fft's of 1301 length 15. Between the fft's a modulation is calculated. 1302 */ 1303 static inline void fft240(FIXP_DBL *pInput, INT *pScalefactor) 1304 { 1305 fftN2( 1306 pInput, 240, 16, 15, 1307 fft_16, fft15, 1308 RotVectorReal240, RotVectorImag240 1309 ); 1310 *pScalefactor += SCALEFACTOR240; 1311 } 1312 1313 1314 #define SCALEFACTOR480 8 1315 #define N32 32 1316 #define TABLE_SIZE_16 (32/2) 1317 1318 /** 1319 The function performs the fft of length 480. It is splittet into fft's of length 32 and fft's of 1320 length 15. Between the fft's a modulation is calculated. 1321 */ 1322 static inline void fft480(FIXP_DBL *pInput, INT *pScalefactor) 1323 { 1324 fftN2( 1325 pInput, 480, 32, 15, 1326 fft_32, fft15, 1327 RotVectorReal480, RotVectorImag480 1328 ); 1329 *pScalefactor += SCALEFACTOR480; 1330 } 1331 1332 void fft(int length, FIXP_DBL *pInput, INT *pScalefactor) 1333 { 1334 if (length == 32) 1335 { 1336 fft_32(pInput); 1337 *pScalefactor += SCALEFACTOR32; 1338 } 1339 else 1340 { 1341 1342 switch (length) { 1343 case 16: 1344 fft_16(pInput); 1345 *pScalefactor += SCALEFACTOR16; 1346 break; 1347 case 8: 1348 fft_8(pInput); 1349 *pScalefactor += SCALEFACTOR8; 1350 break; 1351 case 3: 1352 fft3(pInput); 1353 break; 1354 case 4: 1355 fft_4(pInput); 1356 *pScalefactor += SCALEFACTOR4; 1357 break; 1358 case 5: 1359 fft5(pInput); 1360 break; 1361 case 15: 1362 fft15(pInput); 1363 *pScalefactor += 2; 1364 break; 1365 case 60: 1366 fft60(pInput, pScalefactor); 1367 break; 1368 case 64: 1369 dit_fft(pInput, 6, SineTable512, 512); 1370 *pScalefactor += SCALEFACTOR64; 1371 break; 1372 case 240: 1373 fft240(pInput, pScalefactor); 1374 break; 1375 case 256: 1376 dit_fft(pInput, 8, SineTable512, 512); 1377 *pScalefactor += SCALEFACTOR256; 1378 break; 1379 case 480: 1380 fft480(pInput, pScalefactor); 1381 break; 1382 case 512: 1383 dit_fft(pInput, 9, SineTable512, 512); 1384 *pScalefactor += SCALEFACTOR512; 1385 break; 1386 default: 1387 FDK_ASSERT(0); /* FFT length not supported! */ 1388 break; 1389 } 1390 } 1391 } 1392 1393 1394 void ifft(int length, FIXP_DBL *pInput, INT *scalefactor) 1395 { 1396 switch (length) { 1397 default: 1398 FDK_ASSERT(0); /* IFFT length not supported! */ 1399 break; 1400 } 1401 } 1402 1403 1404
