1 /*Copyright (c) 2003-2004, Mark Borgerding 2 Lots of modifications by Jean-Marc Valin 3 Copyright (c) 2005-2007, Xiph.Org Foundation 4 Copyright (c) 2008, Xiph.Org Foundation, CSIRO 5 6 All rights reserved. 7 8 Redistribution and use in source and binary forms, with or without 9 modification, are permitted provided that the following conditions are met: 10 11 * Redistributions of source code must retain the above copyright notice, 12 this list of conditions and the following disclaimer. 13 * Redistributions in binary form must reproduce the above copyright notice, 14 this list of conditions and the following disclaimer in the 15 documentation and/or other materials provided with the distribution. 16 17 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 18 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 19 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 20 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE 21 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 22 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 23 SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 24 INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 25 CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 26 ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 27 POSSIBILITY OF SUCH DAMAGE.*/ 28 29 /* This code is originally from Mark Borgerding's KISS-FFT but has been 30 heavily modified to better suit Opus */ 31 32 #ifndef SKIP_CONFIG_H 33 # ifdef HAVE_CONFIG_H 34 # include "config.h" 35 # endif 36 #endif 37 38 #include "_kiss_fft_guts.h" 39 #include "arch.h" 40 #include "os_support.h" 41 #include "mathops.h" 42 #include "stack_alloc.h" 43 44 /* The guts header contains all the multiplication and addition macros that are defined for 45 complex numbers. It also delares the kf_ internal functions. 46 */ 47 48 static void kf_bfly2( 49 kiss_fft_cpx * Fout, 50 int m, 51 int N 52 ) 53 { 54 kiss_fft_cpx * Fout2; 55 int i; 56 (void)m; 57 #ifdef CUSTOM_MODES 58 if (m==1) 59 { 60 celt_assert(m==1); 61 for (i=0;i<N;i++) 62 { 63 kiss_fft_cpx t; 64 Fout2 = Fout + 1; 65 t = *Fout2; 66 C_SUB( *Fout2 , *Fout , t ); 67 C_ADDTO( *Fout , t ); 68 Fout += 2; 69 } 70 } else 71 #endif 72 { 73 opus_val16 tw; 74 tw = QCONST16(0.7071067812f, 15); 75 /* We know that m==4 here because the radix-2 is just after a radix-4 */ 76 celt_assert(m==4); 77 for (i=0;i<N;i++) 78 { 79 kiss_fft_cpx t; 80 Fout2 = Fout + 4; 81 t = Fout2[0]; 82 C_SUB( Fout2[0] , Fout[0] , t ); 83 C_ADDTO( Fout[0] , t ); 84 85 t.r = S_MUL(ADD32_ovflw(Fout2[1].r, Fout2[1].i), tw); 86 t.i = S_MUL(SUB32_ovflw(Fout2[1].i, Fout2[1].r), tw); 87 C_SUB( Fout2[1] , Fout[1] , t ); 88 C_ADDTO( Fout[1] , t ); 89 90 t.r = Fout2[2].i; 91 t.i = -Fout2[2].r; 92 C_SUB( Fout2[2] , Fout[2] , t ); 93 C_ADDTO( Fout[2] , t ); 94 95 t.r = S_MUL(SUB32_ovflw(Fout2[3].i, Fout2[3].r), tw); 96 t.i = S_MUL(NEG32_ovflw(ADD32_ovflw(Fout2[3].i, Fout2[3].r)), tw); 97 C_SUB( Fout2[3] , Fout[3] , t ); 98 C_ADDTO( Fout[3] , t ); 99 Fout += 8; 100 } 101 } 102 } 103 104 static void kf_bfly4( 105 kiss_fft_cpx * Fout, 106 const size_t fstride, 107 const kiss_fft_state *st, 108 int m, 109 int N, 110 int mm 111 ) 112 { 113 int i; 114 115 if (m==1) 116 { 117 /* Degenerate case where all the twiddles are 1. */ 118 for (i=0;i<N;i++) 119 { 120 kiss_fft_cpx scratch0, scratch1; 121 122 C_SUB( scratch0 , *Fout, Fout[2] ); 123 C_ADDTO(*Fout, Fout[2]); 124 C_ADD( scratch1 , Fout[1] , Fout[3] ); 125 C_SUB( Fout[2], *Fout, scratch1 ); 126 C_ADDTO( *Fout , scratch1 ); 127 C_SUB( scratch1 , Fout[1] , Fout[3] ); 128 129 Fout[1].r = ADD32_ovflw(scratch0.r, scratch1.i); 130 Fout[1].i = SUB32_ovflw(scratch0.i, scratch1.r); 131 Fout[3].r = SUB32_ovflw(scratch0.r, scratch1.i); 132 Fout[3].i = ADD32_ovflw(scratch0.i, scratch1.r); 133 Fout+=4; 134 } 135 } else { 136 int j; 137 kiss_fft_cpx scratch[6]; 138 const kiss_twiddle_cpx *tw1,*tw2,*tw3; 139 const int m2=2*m; 140 const int m3=3*m; 141 kiss_fft_cpx * Fout_beg = Fout; 142 for (i=0;i<N;i++) 143 { 144 Fout = Fout_beg + i*mm; 145 tw3 = tw2 = tw1 = st->twiddles; 146 /* m is guaranteed to be a multiple of 4. */ 147 for (j=0;j<m;j++) 148 { 149 C_MUL(scratch[0],Fout[m] , *tw1 ); 150 C_MUL(scratch[1],Fout[m2] , *tw2 ); 151 C_MUL(scratch[2],Fout[m3] , *tw3 ); 152 153 C_SUB( scratch[5] , *Fout, scratch[1] ); 154 C_ADDTO(*Fout, scratch[1]); 155 C_ADD( scratch[3] , scratch[0] , scratch[2] ); 156 C_SUB( scratch[4] , scratch[0] , scratch[2] ); 157 C_SUB( Fout[m2], *Fout, scratch[3] ); 158 tw1 += fstride; 159 tw2 += fstride*2; 160 tw3 += fstride*3; 161 C_ADDTO( *Fout , scratch[3] ); 162 163 Fout[m].r = ADD32_ovflw(scratch[5].r, scratch[4].i); 164 Fout[m].i = SUB32_ovflw(scratch[5].i, scratch[4].r); 165 Fout[m3].r = SUB32_ovflw(scratch[5].r, scratch[4].i); 166 Fout[m3].i = ADD32_ovflw(scratch[5].i, scratch[4].r); 167 ++Fout; 168 } 169 } 170 } 171 } 172 173 174 #ifndef RADIX_TWO_ONLY 175 176 static void kf_bfly3( 177 kiss_fft_cpx * Fout, 178 const size_t fstride, 179 const kiss_fft_state *st, 180 int m, 181 int N, 182 int mm 183 ) 184 { 185 int i; 186 size_t k; 187 const size_t m2 = 2*m; 188 const kiss_twiddle_cpx *tw1,*tw2; 189 kiss_fft_cpx scratch[5]; 190 kiss_twiddle_cpx epi3; 191 192 kiss_fft_cpx * Fout_beg = Fout; 193 #ifdef FIXED_POINT 194 /*epi3.r = -16384;*/ /* Unused */ 195 epi3.i = -28378; 196 #else 197 epi3 = st->twiddles[fstride*m]; 198 #endif 199 for (i=0;i<N;i++) 200 { 201 Fout = Fout_beg + i*mm; 202 tw1=tw2=st->twiddles; 203 /* For non-custom modes, m is guaranteed to be a multiple of 4. */ 204 k=m; 205 do { 206 207 C_MUL(scratch[1],Fout[m] , *tw1); 208 C_MUL(scratch[2],Fout[m2] , *tw2); 209 210 C_ADD(scratch[3],scratch[1],scratch[2]); 211 C_SUB(scratch[0],scratch[1],scratch[2]); 212 tw1 += fstride; 213 tw2 += fstride*2; 214 215 Fout[m].r = SUB32_ovflw(Fout->r, HALF_OF(scratch[3].r)); 216 Fout[m].i = SUB32_ovflw(Fout->i, HALF_OF(scratch[3].i)); 217 218 C_MULBYSCALAR( scratch[0] , epi3.i ); 219 220 C_ADDTO(*Fout,scratch[3]); 221 222 Fout[m2].r = ADD32_ovflw(Fout[m].r, scratch[0].i); 223 Fout[m2].i = SUB32_ovflw(Fout[m].i, scratch[0].r); 224 225 Fout[m].r = SUB32_ovflw(Fout[m].r, scratch[0].i); 226 Fout[m].i = ADD32_ovflw(Fout[m].i, scratch[0].r); 227 228 ++Fout; 229 } while(--k); 230 } 231 } 232 233 234 #ifndef OVERRIDE_kf_bfly5 235 static void kf_bfly5( 236 kiss_fft_cpx * Fout, 237 const size_t fstride, 238 const kiss_fft_state *st, 239 int m, 240 int N, 241 int mm 242 ) 243 { 244 kiss_fft_cpx *Fout0,*Fout1,*Fout2,*Fout3,*Fout4; 245 int i, u; 246 kiss_fft_cpx scratch[13]; 247 const kiss_twiddle_cpx *tw; 248 kiss_twiddle_cpx ya,yb; 249 kiss_fft_cpx * Fout_beg = Fout; 250 251 #ifdef FIXED_POINT 252 ya.r = 10126; 253 ya.i = -31164; 254 yb.r = -26510; 255 yb.i = -19261; 256 #else 257 ya = st->twiddles[fstride*m]; 258 yb = st->twiddles[fstride*2*m]; 259 #endif 260 tw=st->twiddles; 261 262 for (i=0;i<N;i++) 263 { 264 Fout = Fout_beg + i*mm; 265 Fout0=Fout; 266 Fout1=Fout0+m; 267 Fout2=Fout0+2*m; 268 Fout3=Fout0+3*m; 269 Fout4=Fout0+4*m; 270 271 /* For non-custom modes, m is guaranteed to be a multiple of 4. */ 272 for ( u=0; u<m; ++u ) { 273 scratch[0] = *Fout0; 274 275 C_MUL(scratch[1] ,*Fout1, tw[u*fstride]); 276 C_MUL(scratch[2] ,*Fout2, tw[2*u*fstride]); 277 C_MUL(scratch[3] ,*Fout3, tw[3*u*fstride]); 278 C_MUL(scratch[4] ,*Fout4, tw[4*u*fstride]); 279 280 C_ADD( scratch[7],scratch[1],scratch[4]); 281 C_SUB( scratch[10],scratch[1],scratch[4]); 282 C_ADD( scratch[8],scratch[2],scratch[3]); 283 C_SUB( scratch[9],scratch[2],scratch[3]); 284 285 Fout0->r = ADD32_ovflw(Fout0->r, ADD32_ovflw(scratch[7].r, scratch[8].r)); 286 Fout0->i = ADD32_ovflw(Fout0->i, ADD32_ovflw(scratch[7].i, scratch[8].i)); 287 288 scratch[5].r = ADD32_ovflw(scratch[0].r, ADD32_ovflw(S_MUL(scratch[7].r,ya.r), S_MUL(scratch[8].r,yb.r))); 289 scratch[5].i = ADD32_ovflw(scratch[0].i, ADD32_ovflw(S_MUL(scratch[7].i,ya.r), S_MUL(scratch[8].i,yb.r))); 290 291 scratch[6].r = ADD32_ovflw(S_MUL(scratch[10].i,ya.i), S_MUL(scratch[9].i,yb.i)); 292 scratch[6].i = NEG32_ovflw(ADD32_ovflw(S_MUL(scratch[10].r,ya.i), S_MUL(scratch[9].r,yb.i))); 293 294 C_SUB(*Fout1,scratch[5],scratch[6]); 295 C_ADD(*Fout4,scratch[5],scratch[6]); 296 297 scratch[11].r = ADD32_ovflw(scratch[0].r, ADD32_ovflw(S_MUL(scratch[7].r,yb.r), S_MUL(scratch[8].r,ya.r))); 298 scratch[11].i = ADD32_ovflw(scratch[0].i, ADD32_ovflw(S_MUL(scratch[7].i,yb.r), S_MUL(scratch[8].i,ya.r))); 299 scratch[12].r = SUB32_ovflw(S_MUL(scratch[9].i,ya.i), S_MUL(scratch[10].i,yb.i)); 300 scratch[12].i = SUB32_ovflw(S_MUL(scratch[10].r,yb.i), S_MUL(scratch[9].r,ya.i)); 301 302 C_ADD(*Fout2,scratch[11],scratch[12]); 303 C_SUB(*Fout3,scratch[11],scratch[12]); 304 305 ++Fout0;++Fout1;++Fout2;++Fout3;++Fout4; 306 } 307 } 308 } 309 #endif /* OVERRIDE_kf_bfly5 */ 310 311 312 #endif 313 314 315 #ifdef CUSTOM_MODES 316 317 static 318 void compute_bitrev_table( 319 int Fout, 320 opus_int16 *f, 321 const size_t fstride, 322 int in_stride, 323 opus_int16 * factors, 324 const kiss_fft_state *st 325 ) 326 { 327 const int p=*factors++; /* the radix */ 328 const int m=*factors++; /* stage's fft length/p */ 329 330 /*printf ("fft %d %d %d %d %d %d\n", p*m, m, p, s2, fstride*in_stride, N);*/ 331 if (m==1) 332 { 333 int j; 334 for (j=0;j<p;j++) 335 { 336 *f = Fout+j; 337 f += fstride*in_stride; 338 } 339 } else { 340 int j; 341 for (j=0;j<p;j++) 342 { 343 compute_bitrev_table( Fout , f, fstride*p, in_stride, factors,st); 344 f += fstride*in_stride; 345 Fout += m; 346 } 347 } 348 } 349 350 /* facbuf is populated by p1,m1,p2,m2, ... 351 where 352 p[i] * m[i] = m[i-1] 353 m0 = n */ 354 static 355 int kf_factor(int n,opus_int16 * facbuf) 356 { 357 int p=4; 358 int i; 359 int stages=0; 360 int nbak = n; 361 362 /*factor out powers of 4, powers of 2, then any remaining primes */ 363 do { 364 while (n % p) { 365 switch (p) { 366 case 4: p = 2; break; 367 case 2: p = 3; break; 368 default: p += 2; break; 369 } 370 if (p>32000 || (opus_int32)p*(opus_int32)p > n) 371 p = n; /* no more factors, skip to end */ 372 } 373 n /= p; 374 #ifdef RADIX_TWO_ONLY 375 if (p!=2 && p != 4) 376 #else 377 if (p>5) 378 #endif 379 { 380 return 0; 381 } 382 facbuf[2*stages] = p; 383 if (p==2 && stages > 1) 384 { 385 facbuf[2*stages] = 4; 386 facbuf[2] = 2; 387 } 388 stages++; 389 } while (n > 1); 390 n = nbak; 391 /* Reverse the order to get the radix 4 at the end, so we can use the 392 fast degenerate case. It turns out that reversing the order also 393 improves the noise behaviour. */ 394 for (i=0;i<stages/2;i++) 395 { 396 int tmp; 397 tmp = facbuf[2*i]; 398 facbuf[2*i] = facbuf[2*(stages-i-1)]; 399 facbuf[2*(stages-i-1)] = tmp; 400 } 401 for (i=0;i<stages;i++) 402 { 403 n /= facbuf[2*i]; 404 facbuf[2*i+1] = n; 405 } 406 return 1; 407 } 408 409 static void compute_twiddles(kiss_twiddle_cpx *twiddles, int nfft) 410 { 411 int i; 412 #ifdef FIXED_POINT 413 for (i=0;i<nfft;++i) { 414 opus_val32 phase = -i; 415 kf_cexp2(twiddles+i, DIV32(SHL32(phase,17),nfft)); 416 } 417 #else 418 for (i=0;i<nfft;++i) { 419 const double pi=3.14159265358979323846264338327; 420 double phase = ( -2*pi /nfft ) * i; 421 kf_cexp(twiddles+i, phase ); 422 } 423 #endif 424 } 425 426 int opus_fft_alloc_arch_c(kiss_fft_state *st) { 427 (void)st; 428 return 0; 429 } 430 431 /* 432 * 433 * Allocates all necessary storage space for the fft and ifft. 434 * The return value is a contiguous block of memory. As such, 435 * It can be freed with free(). 436 * */ 437 kiss_fft_state *opus_fft_alloc_twiddles(int nfft,void * mem,size_t * lenmem, 438 const kiss_fft_state *base, int arch) 439 { 440 kiss_fft_state *st=NULL; 441 size_t memneeded = sizeof(struct kiss_fft_state); /* twiddle factors*/ 442 443 if ( lenmem==NULL ) { 444 st = ( kiss_fft_state*)KISS_FFT_MALLOC( memneeded ); 445 }else{ 446 if (mem != NULL && *lenmem >= memneeded) 447 st = (kiss_fft_state*)mem; 448 *lenmem = memneeded; 449 } 450 if (st) { 451 opus_int16 *bitrev; 452 kiss_twiddle_cpx *twiddles; 453 454 st->nfft=nfft; 455 #ifdef FIXED_POINT 456 st->scale_shift = celt_ilog2(st->nfft); 457 if (st->nfft == 1<<st->scale_shift) 458 st->scale = Q15ONE; 459 else 460 st->scale = (1073741824+st->nfft/2)/st->nfft>>(15-st->scale_shift); 461 #else 462 st->scale = 1.f/nfft; 463 #endif 464 if (base != NULL) 465 { 466 st->twiddles = base->twiddles; 467 st->shift = 0; 468 while (st->shift < 32 && nfft<<st->shift != base->nfft) 469 st->shift++; 470 if (st->shift>=32) 471 goto fail; 472 } else { 473 st->twiddles = twiddles = (kiss_twiddle_cpx*)KISS_FFT_MALLOC(sizeof(kiss_twiddle_cpx)*nfft); 474 compute_twiddles(twiddles, nfft); 475 st->shift = -1; 476 } 477 if (!kf_factor(nfft,st->factors)) 478 { 479 goto fail; 480 } 481 482 /* bitrev */ 483 st->bitrev = bitrev = (opus_int16*)KISS_FFT_MALLOC(sizeof(opus_int16)*nfft); 484 if (st->bitrev==NULL) 485 goto fail; 486 compute_bitrev_table(0, bitrev, 1,1, st->factors,st); 487 488 /* Initialize architecture specific fft parameters */ 489 if (opus_fft_alloc_arch(st, arch)) 490 goto fail; 491 } 492 return st; 493 fail: 494 opus_fft_free(st, arch); 495 return NULL; 496 } 497 498 kiss_fft_state *opus_fft_alloc(int nfft,void * mem,size_t * lenmem, int arch) 499 { 500 return opus_fft_alloc_twiddles(nfft, mem, lenmem, NULL, arch); 501 } 502 503 void opus_fft_free_arch_c(kiss_fft_state *st) { 504 (void)st; 505 } 506 507 void opus_fft_free(const kiss_fft_state *cfg, int arch) 508 { 509 if (cfg) 510 { 511 opus_fft_free_arch((kiss_fft_state *)cfg, arch); 512 opus_free((opus_int16*)cfg->bitrev); 513 if (cfg->shift < 0) 514 opus_free((kiss_twiddle_cpx*)cfg->twiddles); 515 opus_free((kiss_fft_state*)cfg); 516 } 517 } 518 519 #endif /* CUSTOM_MODES */ 520 521 void opus_fft_impl(const kiss_fft_state *st,kiss_fft_cpx *fout) 522 { 523 int m2, m; 524 int p; 525 int L; 526 int fstride[MAXFACTORS]; 527 int i; 528 int shift; 529 530 /* st->shift can be -1 */ 531 shift = st->shift>0 ? st->shift : 0; 532 533 fstride[0] = 1; 534 L=0; 535 do { 536 p = st->factors[2*L]; 537 m = st->factors[2*L+1]; 538 fstride[L+1] = fstride[L]*p; 539 L++; 540 } while(m!=1); 541 m = st->factors[2*L-1]; 542 for (i=L-1;i>=0;i--) 543 { 544 if (i!=0) 545 m2 = st->factors[2*i-1]; 546 else 547 m2 = 1; 548 switch (st->factors[2*i]) 549 { 550 case 2: 551 kf_bfly2(fout, m, fstride[i]); 552 break; 553 case 4: 554 kf_bfly4(fout,fstride[i]<<shift,st,m, fstride[i], m2); 555 break; 556 #ifndef RADIX_TWO_ONLY 557 case 3: 558 kf_bfly3(fout,fstride[i]<<shift,st,m, fstride[i], m2); 559 break; 560 case 5: 561 kf_bfly5(fout,fstride[i]<<shift,st,m, fstride[i], m2); 562 break; 563 #endif 564 } 565 m = m2; 566 } 567 } 568 569 void opus_fft_c(const kiss_fft_state *st,const kiss_fft_cpx *fin,kiss_fft_cpx *fout) 570 { 571 int i; 572 opus_val16 scale; 573 #ifdef FIXED_POINT 574 /* Allows us to scale with MULT16_32_Q16(), which is faster than 575 MULT16_32_Q15() on ARM. */ 576 int scale_shift = st->scale_shift-1; 577 #endif 578 scale = st->scale; 579 580 celt_assert2 (fin != fout, "In-place FFT not supported"); 581 /* Bit-reverse the input */ 582 for (i=0;i<st->nfft;i++) 583 { 584 kiss_fft_cpx x = fin[i]; 585 fout[st->bitrev[i]].r = SHR32(MULT16_32_Q16(scale, x.r), scale_shift); 586 fout[st->bitrev[i]].i = SHR32(MULT16_32_Q16(scale, x.i), scale_shift); 587 } 588 opus_fft_impl(st, fout); 589 } 590 591 592 void opus_ifft_c(const kiss_fft_state *st,const kiss_fft_cpx *fin,kiss_fft_cpx *fout) 593 { 594 int i; 595 celt_assert2 (fin != fout, "In-place FFT not supported"); 596 /* Bit-reverse the input */ 597 for (i=0;i<st->nfft;i++) 598 fout[st->bitrev[i]] = fin[i]; 599 for (i=0;i<st->nfft;i++) 600 fout[i].i = -fout[i].i; 601 opus_fft_impl(st, fout); 602 for (i=0;i<st->nfft;i++) 603 fout[i].i = -fout[i].i; 604 } 605
