1 /******************************************************************** 2 * * 3 * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE. * 4 * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS * 5 * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE * 6 * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING. * 7 * * 8 * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2009 * 9 * by the Xiph.Org Foundation http://www.xiph.org/ * 10 * * 11 ******************************************************************** 12 13 function: *unnormalized* fft transform 14 last mod: $Id: smallft.c 16227 2009-07-08 06:58:46Z xiphmont $ 15 16 ********************************************************************/ 17 18 /* FFT implementation from OggSquish, minus cosine transforms, 19 * minus all but radix 2/4 case. In Vorbis we only need this 20 * cut-down version. 21 * 22 * To do more than just power-of-two sized vectors, see the full 23 * version I wrote for NetLib. 24 * 25 * Note that the packing is a little strange; rather than the FFT r/i 26 * packing following R_0, I_n, R_1, I_1, R_2, I_2 ... R_n-1, I_n-1, 27 * it follows R_0, R_1, I_1, R_2, I_2 ... R_n-1, I_n-1, I_n like the 28 * FORTRAN version 29 */ 30 31 #include <stdlib.h> 32 #include <string.h> 33 #include <math.h> 34 #include "smallft.h" 35 #include "os.h" 36 #include "misc.h" 37 38 static void drfti1(int n, float *wa, int *ifac){ 39 static int ntryh[4] = { 4,2,3,5 }; 40 static float tpi = 6.28318530717958648f; 41 float arg,argh,argld,fi; 42 int ntry=0,i,j=-1; 43 int k1, l1, l2, ib; 44 int ld, ii, ip, is, nq, nr; 45 int ido, ipm, nfm1; 46 int nl=n; 47 int nf=0; 48 49 L101: 50 j++; 51 if (j < 4) 52 ntry=ntryh[j]; 53 else 54 ntry+=2; 55 56 L104: 57 nq=nl/ntry; 58 nr=nl-ntry*nq; 59 if (nr!=0) goto L101; 60 61 nf++; 62 ifac[nf+1]=ntry; 63 nl=nq; 64 if(ntry!=2)goto L107; 65 if(nf==1)goto L107; 66 67 for (i=1;i<nf;i++){ 68 ib=nf-i+1; 69 ifac[ib+1]=ifac[ib]; 70 } 71 ifac[2] = 2; 72 73 L107: 74 if(nl!=1)goto L104; 75 ifac[0]=n; 76 ifac[1]=nf; 77 argh=tpi/n; 78 is=0; 79 nfm1=nf-1; 80 l1=1; 81 82 if(nfm1==0)return; 83 84 for (k1=0;k1<nfm1;k1++){ 85 ip=ifac[k1+2]; 86 ld=0; 87 l2=l1*ip; 88 ido=n/l2; 89 ipm=ip-1; 90 91 for (j=0;j<ipm;j++){ 92 ld+=l1; 93 i=is; 94 argld=(float)ld*argh; 95 fi=0.f; 96 for (ii=2;ii<ido;ii+=2){ 97 fi+=1.f; 98 arg=fi*argld; 99 wa[i++]=cos(arg); 100 wa[i++]=sin(arg); 101 } 102 is+=ido; 103 } 104 l1=l2; 105 } 106 } 107 108 static void fdrffti(int n, float *wsave, int *ifac){ 109 110 if (n == 1) return; 111 drfti1(n, wsave+n, ifac); 112 } 113 114 static void dradf2(int ido,int l1,float *cc,float *ch,float *wa1){ 115 int i,k; 116 float ti2,tr2; 117 int t0,t1,t2,t3,t4,t5,t6; 118 119 t1=0; 120 t0=(t2=l1*ido); 121 t3=ido<<1; 122 for(k=0;k<l1;k++){ 123 ch[t1<<1]=cc[t1]+cc[t2]; 124 ch[(t1<<1)+t3-1]=cc[t1]-cc[t2]; 125 t1+=ido; 126 t2+=ido; 127 } 128 129 if(ido<2)return; 130 if(ido==2)goto L105; 131 132 t1=0; 133 t2=t0; 134 for(k=0;k<l1;k++){ 135 t3=t2; 136 t4=(t1<<1)+(ido<<1); 137 t5=t1; 138 t6=t1+t1; 139 for(i=2;i<ido;i+=2){ 140 t3+=2; 141 t4-=2; 142 t5+=2; 143 t6+=2; 144 tr2=wa1[i-2]*cc[t3-1]+wa1[i-1]*cc[t3]; 145 ti2=wa1[i-2]*cc[t3]-wa1[i-1]*cc[t3-1]; 146 ch[t6]=cc[t5]+ti2; 147 ch[t4]=ti2-cc[t5]; 148 ch[t6-1]=cc[t5-1]+tr2; 149 ch[t4-1]=cc[t5-1]-tr2; 150 } 151 t1+=ido; 152 t2+=ido; 153 } 154 155 if(ido%2==1)return; 156 157 L105: 158 t3=(t2=(t1=ido)-1); 159 t2+=t0; 160 for(k=0;k<l1;k++){ 161 ch[t1]=-cc[t2]; 162 ch[t1-1]=cc[t3]; 163 t1+=ido<<1; 164 t2+=ido; 165 t3+=ido; 166 } 167 } 168 169 static void dradf4(int ido,int l1,float *cc,float *ch,float *wa1, 170 float *wa2,float *wa3){ 171 static float hsqt2 = .70710678118654752f; 172 int i,k,t0,t1,t2,t3,t4,t5,t6; 173 float ci2,ci3,ci4,cr2,cr3,cr4,ti1,ti2,ti3,ti4,tr1,tr2,tr3,tr4; 174 t0=l1*ido; 175 176 t1=t0; 177 t4=t1<<1; 178 t2=t1+(t1<<1); 179 t3=0; 180 181 for(k=0;k<l1;k++){ 182 tr1=cc[t1]+cc[t2]; 183 tr2=cc[t3]+cc[t4]; 184 185 ch[t5=t3<<2]=tr1+tr2; 186 ch[(ido<<2)+t5-1]=tr2-tr1; 187 ch[(t5+=(ido<<1))-1]=cc[t3]-cc[t4]; 188 ch[t5]=cc[t2]-cc[t1]; 189 190 t1+=ido; 191 t2+=ido; 192 t3+=ido; 193 t4+=ido; 194 } 195 196 if(ido<2)return; 197 if(ido==2)goto L105; 198 199 200 t1=0; 201 for(k=0;k<l1;k++){ 202 t2=t1; 203 t4=t1<<2; 204 t5=(t6=ido<<1)+t4; 205 for(i=2;i<ido;i+=2){ 206 t3=(t2+=2); 207 t4+=2; 208 t5-=2; 209 210 t3+=t0; 211 cr2=wa1[i-2]*cc[t3-1]+wa1[i-1]*cc[t3]; 212 ci2=wa1[i-2]*cc[t3]-wa1[i-1]*cc[t3-1]; 213 t3+=t0; 214 cr3=wa2[i-2]*cc[t3-1]+wa2[i-1]*cc[t3]; 215 ci3=wa2[i-2]*cc[t3]-wa2[i-1]*cc[t3-1]; 216 t3+=t0; 217 cr4=wa3[i-2]*cc[t3-1]+wa3[i-1]*cc[t3]; 218 ci4=wa3[i-2]*cc[t3]-wa3[i-1]*cc[t3-1]; 219 220 tr1=cr2+cr4; 221 tr4=cr4-cr2; 222 ti1=ci2+ci4; 223 ti4=ci2-ci4; 224 225 ti2=cc[t2]+ci3; 226 ti3=cc[t2]-ci3; 227 tr2=cc[t2-1]+cr3; 228 tr3=cc[t2-1]-cr3; 229 230 ch[t4-1]=tr1+tr2; 231 ch[t4]=ti1+ti2; 232 233 ch[t5-1]=tr3-ti4; 234 ch[t5]=tr4-ti3; 235 236 ch[t4+t6-1]=ti4+tr3; 237 ch[t4+t6]=tr4+ti3; 238 239 ch[t5+t6-1]=tr2-tr1; 240 ch[t5+t6]=ti1-ti2; 241 } 242 t1+=ido; 243 } 244 if(ido&1)return; 245 246 L105: 247 248 t2=(t1=t0+ido-1)+(t0<<1); 249 t3=ido<<2; 250 t4=ido; 251 t5=ido<<1; 252 t6=ido; 253 254 for(k=0;k<l1;k++){ 255 ti1=-hsqt2*(cc[t1]+cc[t2]); 256 tr1=hsqt2*(cc[t1]-cc[t2]); 257 258 ch[t4-1]=tr1+cc[t6-1]; 259 ch[t4+t5-1]=cc[t6-1]-tr1; 260 261 ch[t4]=ti1-cc[t1+t0]; 262 ch[t4+t5]=ti1+cc[t1+t0]; 263 264 t1+=ido; 265 t2+=ido; 266 t4+=t3; 267 t6+=ido; 268 } 269 } 270 271 static void dradfg(int ido,int ip,int l1,int idl1,float *cc,float *c1, 272 float *c2,float *ch,float *ch2,float *wa){ 273 274 static float tpi=6.283185307179586f; 275 int idij,ipph,i,j,k,l,ic,ik,is; 276 int t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10; 277 float dc2,ai1,ai2,ar1,ar2,ds2; 278 int nbd; 279 float dcp,arg,dsp,ar1h,ar2h; 280 int idp2,ipp2; 281 282 arg=tpi/(float)ip; 283 dcp=cos(arg); 284 dsp=sin(arg); 285 ipph=(ip+1)>>1; 286 ipp2=ip; 287 idp2=ido; 288 nbd=(ido-1)>>1; 289 t0=l1*ido; 290 t10=ip*ido; 291 292 if(ido==1)goto L119; 293 for(ik=0;ik<idl1;ik++)ch2[ik]=c2[ik]; 294 295 t1=0; 296 for(j=1;j<ip;j++){ 297 t1+=t0; 298 t2=t1; 299 for(k=0;k<l1;k++){ 300 ch[t2]=c1[t2]; 301 t2+=ido; 302 } 303 } 304 305 is=-ido; 306 t1=0; 307 if(nbd>l1){ 308 for(j=1;j<ip;j++){ 309 t1+=t0; 310 is+=ido; 311 t2= -ido+t1; 312 for(k=0;k<l1;k++){ 313 idij=is-1; 314 t2+=ido; 315 t3=t2; 316 for(i=2;i<ido;i+=2){ 317 idij+=2; 318 t3+=2; 319 ch[t3-1]=wa[idij-1]*c1[t3-1]+wa[idij]*c1[t3]; 320 ch[t3]=wa[idij-1]*c1[t3]-wa[idij]*c1[t3-1]; 321 } 322 } 323 } 324 }else{ 325 326 for(j=1;j<ip;j++){ 327 is+=ido; 328 idij=is-1; 329 t1+=t0; 330 t2=t1; 331 for(i=2;i<ido;i+=2){ 332 idij+=2; 333 t2+=2; 334 t3=t2; 335 for(k=0;k<l1;k++){ 336 ch[t3-1]=wa[idij-1]*c1[t3-1]+wa[idij]*c1[t3]; 337 ch[t3]=wa[idij-1]*c1[t3]-wa[idij]*c1[t3-1]; 338 t3+=ido; 339 } 340 } 341 } 342 } 343 344 t1=0; 345 t2=ipp2*t0; 346 if(nbd<l1){ 347 for(j=1;j<ipph;j++){ 348 t1+=t0; 349 t2-=t0; 350 t3=t1; 351 t4=t2; 352 for(i=2;i<ido;i+=2){ 353 t3+=2; 354 t4+=2; 355 t5=t3-ido; 356 t6=t4-ido; 357 for(k=0;k<l1;k++){ 358 t5+=ido; 359 t6+=ido; 360 c1[t5-1]=ch[t5-1]+ch[t6-1]; 361 c1[t6-1]=ch[t5]-ch[t6]; 362 c1[t5]=ch[t5]+ch[t6]; 363 c1[t6]=ch[t6-1]-ch[t5-1]; 364 } 365 } 366 } 367 }else{ 368 for(j=1;j<ipph;j++){ 369 t1+=t0; 370 t2-=t0; 371 t3=t1; 372 t4=t2; 373 for(k=0;k<l1;k++){ 374 t5=t3; 375 t6=t4; 376 for(i=2;i<ido;i+=2){ 377 t5+=2; 378 t6+=2; 379 c1[t5-1]=ch[t5-1]+ch[t6-1]; 380 c1[t6-1]=ch[t5]-ch[t6]; 381 c1[t5]=ch[t5]+ch[t6]; 382 c1[t6]=ch[t6-1]-ch[t5-1]; 383 } 384 t3+=ido; 385 t4+=ido; 386 } 387 } 388 } 389 390 L119: 391 for(ik=0;ik<idl1;ik++)c2[ik]=ch2[ik]; 392 393 t1=0; 394 t2=ipp2*idl1; 395 for(j=1;j<ipph;j++){ 396 t1+=t0; 397 t2-=t0; 398 t3=t1-ido; 399 t4=t2-ido; 400 for(k=0;k<l1;k++){ 401 t3+=ido; 402 t4+=ido; 403 c1[t3]=ch[t3]+ch[t4]; 404 c1[t4]=ch[t4]-ch[t3]; 405 } 406 } 407 408 ar1=1.f; 409 ai1=0.f; 410 t1=0; 411 t2=ipp2*idl1; 412 t3=(ip-1)*idl1; 413 for(l=1;l<ipph;l++){ 414 t1+=idl1; 415 t2-=idl1; 416 ar1h=dcp*ar1-dsp*ai1; 417 ai1=dcp*ai1+dsp*ar1; 418 ar1=ar1h; 419 t4=t1; 420 t5=t2; 421 t6=t3; 422 t7=idl1; 423 424 for(ik=0;ik<idl1;ik++){ 425 ch2[t4++]=c2[ik]+ar1*c2[t7++]; 426 ch2[t5++]=ai1*c2[t6++]; 427 } 428 429 dc2=ar1; 430 ds2=ai1; 431 ar2=ar1; 432 ai2=ai1; 433 434 t4=idl1; 435 t5=(ipp2-1)*idl1; 436 for(j=2;j<ipph;j++){ 437 t4+=idl1; 438 t5-=idl1; 439 440 ar2h=dc2*ar2-ds2*ai2; 441 ai2=dc2*ai2+ds2*ar2; 442 ar2=ar2h; 443 444 t6=t1; 445 t7=t2; 446 t8=t4; 447 t9=t5; 448 for(ik=0;ik<idl1;ik++){ 449 ch2[t6++]+=ar2*c2[t8++]; 450 ch2[t7++]+=ai2*c2[t9++]; 451 } 452 } 453 } 454 455 t1=0; 456 for(j=1;j<ipph;j++){ 457 t1+=idl1; 458 t2=t1; 459 for(ik=0;ik<idl1;ik++)ch2[ik]+=c2[t2++]; 460 } 461 462 if(ido<l1)goto L132; 463 464 t1=0; 465 t2=0; 466 for(k=0;k<l1;k++){ 467 t3=t1; 468 t4=t2; 469 for(i=0;i<ido;i++)cc[t4++]=ch[t3++]; 470 t1+=ido; 471 t2+=t10; 472 } 473 474 goto L135; 475 476 L132: 477 for(i=0;i<ido;i++){ 478 t1=i; 479 t2=i; 480 for(k=0;k<l1;k++){ 481 cc[t2]=ch[t1]; 482 t1+=ido; 483 t2+=t10; 484 } 485 } 486 487 L135: 488 t1=0; 489 t2=ido<<1; 490 t3=0; 491 t4=ipp2*t0; 492 for(j=1;j<ipph;j++){ 493 494 t1+=t2; 495 t3+=t0; 496 t4-=t0; 497 498 t5=t1; 499 t6=t3; 500 t7=t4; 501 502 for(k=0;k<l1;k++){ 503 cc[t5-1]=ch[t6]; 504 cc[t5]=ch[t7]; 505 t5+=t10; 506 t6+=ido; 507 t7+=ido; 508 } 509 } 510 511 if(ido==1)return; 512 if(nbd<l1)goto L141; 513 514 t1=-ido; 515 t3=0; 516 t4=0; 517 t5=ipp2*t0; 518 for(j=1;j<ipph;j++){ 519 t1+=t2; 520 t3+=t2; 521 t4+=t0; 522 t5-=t0; 523 t6=t1; 524 t7=t3; 525 t8=t4; 526 t9=t5; 527 for(k=0;k<l1;k++){ 528 for(i=2;i<ido;i+=2){ 529 ic=idp2-i; 530 cc[i+t7-1]=ch[i+t8-1]+ch[i+t9-1]; 531 cc[ic+t6-1]=ch[i+t8-1]-ch[i+t9-1]; 532 cc[i+t7]=ch[i+t8]+ch[i+t9]; 533 cc[ic+t6]=ch[i+t9]-ch[i+t8]; 534 } 535 t6+=t10; 536 t7+=t10; 537 t8+=ido; 538 t9+=ido; 539 } 540 } 541 return; 542 543 L141: 544 545 t1=-ido; 546 t3=0; 547 t4=0; 548 t5=ipp2*t0; 549 for(j=1;j<ipph;j++){ 550 t1+=t2; 551 t3+=t2; 552 t4+=t0; 553 t5-=t0; 554 for(i=2;i<ido;i+=2){ 555 t6=idp2+t1-i; 556 t7=i+t3; 557 t8=i+t4; 558 t9=i+t5; 559 for(k=0;k<l1;k++){ 560 cc[t7-1]=ch[t8-1]+ch[t9-1]; 561 cc[t6-1]=ch[t8-1]-ch[t9-1]; 562 cc[t7]=ch[t8]+ch[t9]; 563 cc[t6]=ch[t9]-ch[t8]; 564 t6+=t10; 565 t7+=t10; 566 t8+=ido; 567 t9+=ido; 568 } 569 } 570 } 571 } 572 573 static void drftf1(int n,float *c,float *ch,float *wa,int *ifac){ 574 int i,k1,l1,l2; 575 int na,kh,nf; 576 int ip,iw,ido,idl1,ix2,ix3; 577 578 nf=ifac[1]; 579 na=1; 580 l2=n; 581 iw=n; 582 583 for(k1=0;k1<nf;k1++){ 584 kh=nf-k1; 585 ip=ifac[kh+1]; 586 l1=l2/ip; 587 ido=n/l2; 588 idl1=ido*l1; 589 iw-=(ip-1)*ido; 590 na=1-na; 591 592 if(ip!=4)goto L102; 593 594 ix2=iw+ido; 595 ix3=ix2+ido; 596 if(na!=0) 597 dradf4(ido,l1,ch,c,wa+iw-1,wa+ix2-1,wa+ix3-1); 598 else 599 dradf4(ido,l1,c,ch,wa+iw-1,wa+ix2-1,wa+ix3-1); 600 goto L110; 601 602 L102: 603 if(ip!=2)goto L104; 604 if(na!=0)goto L103; 605 606 dradf2(ido,l1,c,ch,wa+iw-1); 607 goto L110; 608 609 L103: 610 dradf2(ido,l1,ch,c,wa+iw-1); 611 goto L110; 612 613 L104: 614 if(ido==1)na=1-na; 615 if(na!=0)goto L109; 616 617 dradfg(ido,ip,l1,idl1,c,c,c,ch,ch,wa+iw-1); 618 na=1; 619 goto L110; 620 621 L109: 622 dradfg(ido,ip,l1,idl1,ch,ch,ch,c,c,wa+iw-1); 623 na=0; 624 625 L110: 626 l2=l1; 627 } 628 629 if(na==1)return; 630 631 for(i=0;i<n;i++)c[i]=ch[i]; 632 } 633 634 static void dradb2(int ido,int l1,float *cc,float *ch,float *wa1){ 635 int i,k,t0,t1,t2,t3,t4,t5,t6; 636 float ti2,tr2; 637 638 t0=l1*ido; 639 640 t1=0; 641 t2=0; 642 t3=(ido<<1)-1; 643 for(k=0;k<l1;k++){ 644 ch[t1]=cc[t2]+cc[t3+t2]; 645 ch[t1+t0]=cc[t2]-cc[t3+t2]; 646 t2=(t1+=ido)<<1; 647 } 648 649 if(ido<2)return; 650 if(ido==2)goto L105; 651 652 t1=0; 653 t2=0; 654 for(k=0;k<l1;k++){ 655 t3=t1; 656 t5=(t4=t2)+(ido<<1); 657 t6=t0+t1; 658 for(i=2;i<ido;i+=2){ 659 t3+=2; 660 t4+=2; 661 t5-=2; 662 t6+=2; 663 ch[t3-1]=cc[t4-1]+cc[t5-1]; 664 tr2=cc[t4-1]-cc[t5-1]; 665 ch[t3]=cc[t4]-cc[t5]; 666 ti2=cc[t4]+cc[t5]; 667 ch[t6-1]=wa1[i-2]*tr2-wa1[i-1]*ti2; 668 ch[t6]=wa1[i-2]*ti2+wa1[i-1]*tr2; 669 } 670 t2=(t1+=ido)<<1; 671 } 672 673 if(ido%2==1)return; 674 675 L105: 676 t1=ido-1; 677 t2=ido-1; 678 for(k=0;k<l1;k++){ 679 ch[t1]=cc[t2]+cc[t2]; 680 ch[t1+t0]=-(cc[t2+1]+cc[t2+1]); 681 t1+=ido; 682 t2+=ido<<1; 683 } 684 } 685 686 static void dradb3(int ido,int l1,float *cc,float *ch,float *wa1, 687 float *wa2){ 688 static float taur = -.5f; 689 static float taui = .8660254037844386f; 690 int i,k,t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10; 691 float ci2,ci3,di2,di3,cr2,cr3,dr2,dr3,ti2,tr2; 692 t0=l1*ido; 693 694 t1=0; 695 t2=t0<<1; 696 t3=ido<<1; 697 t4=ido+(ido<<1); 698 t5=0; 699 for(k=0;k<l1;k++){ 700 tr2=cc[t3-1]+cc[t3-1]; 701 cr2=cc[t5]+(taur*tr2); 702 ch[t1]=cc[t5]+tr2; 703 ci3=taui*(cc[t3]+cc[t3]); 704 ch[t1+t0]=cr2-ci3; 705 ch[t1+t2]=cr2+ci3; 706 t1+=ido; 707 t3+=t4; 708 t5+=t4; 709 } 710 711 if(ido==1)return; 712 713 t1=0; 714 t3=ido<<1; 715 for(k=0;k<l1;k++){ 716 t7=t1+(t1<<1); 717 t6=(t5=t7+t3); 718 t8=t1; 719 t10=(t9=t1+t0)+t0; 720 721 for(i=2;i<ido;i+=2){ 722 t5+=2; 723 t6-=2; 724 t7+=2; 725 t8+=2; 726 t9+=2; 727 t10+=2; 728 tr2=cc[t5-1]+cc[t6-1]; 729 cr2=cc[t7-1]+(taur*tr2); 730 ch[t8-1]=cc[t7-1]+tr2; 731 ti2=cc[t5]-cc[t6]; 732 ci2=cc[t7]+(taur*ti2); 733 ch[t8]=cc[t7]+ti2; 734 cr3=taui*(cc[t5-1]-cc[t6-1]); 735 ci3=taui*(cc[t5]+cc[t6]); 736 dr2=cr2-ci3; 737 dr3=cr2+ci3; 738 di2=ci2+cr3; 739 di3=ci2-cr3; 740 ch[t9-1]=wa1[i-2]*dr2-wa1[i-1]*di2; 741 ch[t9]=wa1[i-2]*di2+wa1[i-1]*dr2; 742 ch[t10-1]=wa2[i-2]*dr3-wa2[i-1]*di3; 743 ch[t10]=wa2[i-2]*di3+wa2[i-1]*dr3; 744 } 745 t1+=ido; 746 } 747 } 748 749 static void dradb4(int ido,int l1,float *cc,float *ch,float *wa1, 750 float *wa2,float *wa3){ 751 static float sqrt2=1.414213562373095f; 752 int i,k,t0,t1,t2,t3,t4,t5,t6,t7,t8; 753 float ci2,ci3,ci4,cr2,cr3,cr4,ti1,ti2,ti3,ti4,tr1,tr2,tr3,tr4; 754 t0=l1*ido; 755 756 t1=0; 757 t2=ido<<2; 758 t3=0; 759 t6=ido<<1; 760 for(k=0;k<l1;k++){ 761 t4=t3+t6; 762 t5=t1; 763 tr3=cc[t4-1]+cc[t4-1]; 764 tr4=cc[t4]+cc[t4]; 765 tr1=cc[t3]-cc[(t4+=t6)-1]; 766 tr2=cc[t3]+cc[t4-1]; 767 ch[t5]=tr2+tr3; 768 ch[t5+=t0]=tr1-tr4; 769 ch[t5+=t0]=tr2-tr3; 770 ch[t5+=t0]=tr1+tr4; 771 t1+=ido; 772 t3+=t2; 773 } 774 775 if(ido<2)return; 776 if(ido==2)goto L105; 777 778 t1=0; 779 for(k=0;k<l1;k++){ 780 t5=(t4=(t3=(t2=t1<<2)+t6))+t6; 781 t7=t1; 782 for(i=2;i<ido;i+=2){ 783 t2+=2; 784 t3+=2; 785 t4-=2; 786 t5-=2; 787 t7+=2; 788 ti1=cc[t2]+cc[t5]; 789 ti2=cc[t2]-cc[t5]; 790 ti3=cc[t3]-cc[t4]; 791 tr4=cc[t3]+cc[t4]; 792 tr1=cc[t2-1]-cc[t5-1]; 793 tr2=cc[t2-1]+cc[t5-1]; 794 ti4=cc[t3-1]-cc[t4-1]; 795 tr3=cc[t3-1]+cc[t4-1]; 796 ch[t7-1]=tr2+tr3; 797 cr3=tr2-tr3; 798 ch[t7]=ti2+ti3; 799 ci3=ti2-ti3; 800 cr2=tr1-tr4; 801 cr4=tr1+tr4; 802 ci2=ti1+ti4; 803 ci4=ti1-ti4; 804 805 ch[(t8=t7+t0)-1]=wa1[i-2]*cr2-wa1[i-1]*ci2; 806 ch[t8]=wa1[i-2]*ci2+wa1[i-1]*cr2; 807 ch[(t8+=t0)-1]=wa2[i-2]*cr3-wa2[i-1]*ci3; 808 ch[t8]=wa2[i-2]*ci3+wa2[i-1]*cr3; 809 ch[(t8+=t0)-1]=wa3[i-2]*cr4-wa3[i-1]*ci4; 810 ch[t8]=wa3[i-2]*ci4+wa3[i-1]*cr4; 811 } 812 t1+=ido; 813 } 814 815 if(ido%2 == 1)return; 816 817 L105: 818 819 t1=ido; 820 t2=ido<<2; 821 t3=ido-1; 822 t4=ido+(ido<<1); 823 for(k=0;k<l1;k++){ 824 t5=t3; 825 ti1=cc[t1]+cc[t4]; 826 ti2=cc[t4]-cc[t1]; 827 tr1=cc[t1-1]-cc[t4-1]; 828 tr2=cc[t1-1]+cc[t4-1]; 829 ch[t5]=tr2+tr2; 830 ch[t5+=t0]=sqrt2*(tr1-ti1); 831 ch[t5+=t0]=ti2+ti2; 832 ch[t5+=t0]=-sqrt2*(tr1+ti1); 833 834 t3+=ido; 835 t1+=t2; 836 t4+=t2; 837 } 838 } 839 840 static void dradbg(int ido,int ip,int l1,int idl1,float *cc,float *c1, 841 float *c2,float *ch,float *ch2,float *wa){ 842 static float tpi=6.283185307179586f; 843 int idij,ipph,i,j,k,l,ik,is,t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10, 844 t11,t12; 845 float dc2,ai1,ai2,ar1,ar2,ds2; 846 int nbd; 847 float dcp,arg,dsp,ar1h,ar2h; 848 int ipp2; 849 850 t10=ip*ido; 851 t0=l1*ido; 852 arg=tpi/(float)ip; 853 dcp=cos(arg); 854 dsp=sin(arg); 855 nbd=(ido-1)>>1; 856 ipp2=ip; 857 ipph=(ip+1)>>1; 858 if(ido<l1)goto L103; 859 860 t1=0; 861 t2=0; 862 for(k=0;k<l1;k++){ 863 t3=t1; 864 t4=t2; 865 for(i=0;i<ido;i++){ 866 ch[t3]=cc[t4]; 867 t3++; 868 t4++; 869 } 870 t1+=ido; 871 t2+=t10; 872 } 873 goto L106; 874 875 L103: 876 t1=0; 877 for(i=0;i<ido;i++){ 878 t2=t1; 879 t3=t1; 880 for(k=0;k<l1;k++){ 881 ch[t2]=cc[t3]; 882 t2+=ido; 883 t3+=t10; 884 } 885 t1++; 886 } 887 888 L106: 889 t1=0; 890 t2=ipp2*t0; 891 t7=(t5=ido<<1); 892 for(j=1;j<ipph;j++){ 893 t1+=t0; 894 t2-=t0; 895 t3=t1; 896 t4=t2; 897 t6=t5; 898 for(k=0;k<l1;k++){ 899 ch[t3]=cc[t6-1]+cc[t6-1]; 900 ch[t4]=cc[t6]+cc[t6]; 901 t3+=ido; 902 t4+=ido; 903 t6+=t10; 904 } 905 t5+=t7; 906 } 907 908 if (ido == 1)goto L116; 909 if(nbd<l1)goto L112; 910 911 t1=0; 912 t2=ipp2*t0; 913 t7=0; 914 for(j=1;j<ipph;j++){ 915 t1+=t0; 916 t2-=t0; 917 t3=t1; 918 t4=t2; 919 920 t7+=(ido<<1); 921 t8=t7; 922 for(k=0;k<l1;k++){ 923 t5=t3; 924 t6=t4; 925 t9=t8; 926 t11=t8; 927 for(i=2;i<ido;i+=2){ 928 t5+=2; 929 t6+=2; 930 t9+=2; 931 t11-=2; 932 ch[t5-1]=cc[t9-1]+cc[t11-1]; 933 ch[t6-1]=cc[t9-1]-cc[t11-1]; 934 ch[t5]=cc[t9]-cc[t11]; 935 ch[t6]=cc[t9]+cc[t11]; 936 } 937 t3+=ido; 938 t4+=ido; 939 t8+=t10; 940 } 941 } 942 goto L116; 943 944 L112: 945 t1=0; 946 t2=ipp2*t0; 947 t7=0; 948 for(j=1;j<ipph;j++){ 949 t1+=t0; 950 t2-=t0; 951 t3=t1; 952 t4=t2; 953 t7+=(ido<<1); 954 t8=t7; 955 t9=t7; 956 for(i=2;i<ido;i+=2){ 957 t3+=2; 958 t4+=2; 959 t8+=2; 960 t9-=2; 961 t5=t3; 962 t6=t4; 963 t11=t8; 964 t12=t9; 965 for(k=0;k<l1;k++){ 966 ch[t5-1]=cc[t11-1]+cc[t12-1]; 967 ch[t6-1]=cc[t11-1]-cc[t12-1]; 968 ch[t5]=cc[t11]-cc[t12]; 969 ch[t6]=cc[t11]+cc[t12]; 970 t5+=ido; 971 t6+=ido; 972 t11+=t10; 973 t12+=t10; 974 } 975 } 976 } 977 978 L116: 979 ar1=1.f; 980 ai1=0.f; 981 t1=0; 982 t9=(t2=ipp2*idl1); 983 t3=(ip-1)*idl1; 984 for(l=1;l<ipph;l++){ 985 t1+=idl1; 986 t2-=idl1; 987 988 ar1h=dcp*ar1-dsp*ai1; 989 ai1=dcp*ai1+dsp*ar1; 990 ar1=ar1h; 991 t4=t1; 992 t5=t2; 993 t6=0; 994 t7=idl1; 995 t8=t3; 996 for(ik=0;ik<idl1;ik++){ 997 c2[t4++]=ch2[t6++]+ar1*ch2[t7++]; 998 c2[t5++]=ai1*ch2[t8++]; 999 } 1000 dc2=ar1; 1001 ds2=ai1; 1002 ar2=ar1; 1003 ai2=ai1; 1004 1005 t6=idl1; 1006 t7=t9-idl1; 1007 for(j=2;j<ipph;j++){ 1008 t6+=idl1; 1009 t7-=idl1; 1010 ar2h=dc2*ar2-ds2*ai2; 1011 ai2=dc2*ai2+ds2*ar2; 1012 ar2=ar2h; 1013 t4=t1; 1014 t5=t2; 1015 t11=t6; 1016 t12=t7; 1017 for(ik=0;ik<idl1;ik++){ 1018 c2[t4++]+=ar2*ch2[t11++]; 1019 c2[t5++]+=ai2*ch2[t12++]; 1020 } 1021 } 1022 } 1023 1024 t1=0; 1025 for(j=1;j<ipph;j++){ 1026 t1+=idl1; 1027 t2=t1; 1028 for(ik=0;ik<idl1;ik++)ch2[ik]+=ch2[t2++]; 1029 } 1030 1031 t1=0; 1032 t2=ipp2*t0; 1033 for(j=1;j<ipph;j++){ 1034 t1+=t0; 1035 t2-=t0; 1036 t3=t1; 1037 t4=t2; 1038 for(k=0;k<l1;k++){ 1039 ch[t3]=c1[t3]-c1[t4]; 1040 ch[t4]=c1[t3]+c1[t4]; 1041 t3+=ido; 1042 t4+=ido; 1043 } 1044 } 1045 1046 if(ido==1)goto L132; 1047 if(nbd<l1)goto L128; 1048 1049 t1=0; 1050 t2=ipp2*t0; 1051 for(j=1;j<ipph;j++){ 1052 t1+=t0; 1053 t2-=t0; 1054 t3=t1; 1055 t4=t2; 1056 for(k=0;k<l1;k++){ 1057 t5=t3; 1058 t6=t4; 1059 for(i=2;i<ido;i+=2){ 1060 t5+=2; 1061 t6+=2; 1062 ch[t5-1]=c1[t5-1]-c1[t6]; 1063 ch[t6-1]=c1[t5-1]+c1[t6]; 1064 ch[t5]=c1[t5]+c1[t6-1]; 1065 ch[t6]=c1[t5]-c1[t6-1]; 1066 } 1067 t3+=ido; 1068 t4+=ido; 1069 } 1070 } 1071 goto L132; 1072 1073 L128: 1074 t1=0; 1075 t2=ipp2*t0; 1076 for(j=1;j<ipph;j++){ 1077 t1+=t0; 1078 t2-=t0; 1079 t3=t1; 1080 t4=t2; 1081 for(i=2;i<ido;i+=2){ 1082 t3+=2; 1083 t4+=2; 1084 t5=t3; 1085 t6=t4; 1086 for(k=0;k<l1;k++){ 1087 ch[t5-1]=c1[t5-1]-c1[t6]; 1088 ch[t6-1]=c1[t5-1]+c1[t6]; 1089 ch[t5]=c1[t5]+c1[t6-1]; 1090 ch[t6]=c1[t5]-c1[t6-1]; 1091 t5+=ido; 1092 t6+=ido; 1093 } 1094 } 1095 } 1096 1097 L132: 1098 if(ido==1)return; 1099 1100 for(ik=0;ik<idl1;ik++)c2[ik]=ch2[ik]; 1101 1102 t1=0; 1103 for(j=1;j<ip;j++){ 1104 t2=(t1+=t0); 1105 for(k=0;k<l1;k++){ 1106 c1[t2]=ch[t2]; 1107 t2+=ido; 1108 } 1109 } 1110 1111 if(nbd>l1)goto L139; 1112 1113 is= -ido-1; 1114 t1=0; 1115 for(j=1;j<ip;j++){ 1116 is+=ido; 1117 t1+=t0; 1118 idij=is; 1119 t2=t1; 1120 for(i=2;i<ido;i+=2){ 1121 t2+=2; 1122 idij+=2; 1123 t3=t2; 1124 for(k=0;k<l1;k++){ 1125 c1[t3-1]=wa[idij-1]*ch[t3-1]-wa[idij]*ch[t3]; 1126 c1[t3]=wa[idij-1]*ch[t3]+wa[idij]*ch[t3-1]; 1127 t3+=ido; 1128 } 1129 } 1130 } 1131 return; 1132 1133 L139: 1134 is= -ido-1; 1135 t1=0; 1136 for(j=1;j<ip;j++){ 1137 is+=ido; 1138 t1+=t0; 1139 t2=t1; 1140 for(k=0;k<l1;k++){ 1141 idij=is; 1142 t3=t2; 1143 for(i=2;i<ido;i+=2){ 1144 idij+=2; 1145 t3+=2; 1146 c1[t3-1]=wa[idij-1]*ch[t3-1]-wa[idij]*ch[t3]; 1147 c1[t3]=wa[idij-1]*ch[t3]+wa[idij]*ch[t3-1]; 1148 } 1149 t2+=ido; 1150 } 1151 } 1152 } 1153 1154 static void drftb1(int n, float *c, float *ch, float *wa, int *ifac){ 1155 int i,k1,l1,l2; 1156 int na; 1157 int nf,ip,iw,ix2,ix3,ido,idl1; 1158 1159 nf=ifac[1]; 1160 na=0; 1161 l1=1; 1162 iw=1; 1163 1164 for(k1=0;k1<nf;k1++){ 1165 ip=ifac[k1 + 2]; 1166 l2=ip*l1; 1167 ido=n/l2; 1168 idl1=ido*l1; 1169 if(ip!=4)goto L103; 1170 ix2=iw+ido; 1171 ix3=ix2+ido; 1172 1173 if(na!=0) 1174 dradb4(ido,l1,ch,c,wa+iw-1,wa+ix2-1,wa+ix3-1); 1175 else 1176 dradb4(ido,l1,c,ch,wa+iw-1,wa+ix2-1,wa+ix3-1); 1177 na=1-na; 1178 goto L115; 1179 1180 L103: 1181 if(ip!=2)goto L106; 1182 1183 if(na!=0) 1184 dradb2(ido,l1,ch,c,wa+iw-1); 1185 else 1186 dradb2(ido,l1,c,ch,wa+iw-1); 1187 na=1-na; 1188 goto L115; 1189 1190 L106: 1191 if(ip!=3)goto L109; 1192 1193 ix2=iw+ido; 1194 if(na!=0) 1195 dradb3(ido,l1,ch,c,wa+iw-1,wa+ix2-1); 1196 else 1197 dradb3(ido,l1,c,ch,wa+iw-1,wa+ix2-1); 1198 na=1-na; 1199 goto L115; 1200 1201 L109: 1202 /* The radix five case can be translated later..... */ 1203 /* if(ip!=5)goto L112; 1204 1205 ix2=iw+ido; 1206 ix3=ix2+ido; 1207 ix4=ix3+ido; 1208 if(na!=0) 1209 dradb5(ido,l1,ch,c,wa+iw-1,wa+ix2-1,wa+ix3-1,wa+ix4-1); 1210 else 1211 dradb5(ido,l1,c,ch,wa+iw-1,wa+ix2-1,wa+ix3-1,wa+ix4-1); 1212 na=1-na; 1213 goto L115; 1214 1215 L112:*/ 1216 if(na!=0) 1217 dradbg(ido,ip,l1,idl1,ch,ch,ch,c,c,wa+iw-1); 1218 else 1219 dradbg(ido,ip,l1,idl1,c,c,c,ch,ch,wa+iw-1); 1220 if(ido==1)na=1-na; 1221 1222 L115: 1223 l1=l2; 1224 iw+=(ip-1)*ido; 1225 } 1226 1227 if(na==0)return; 1228 1229 for(i=0;i<n;i++)c[i]=ch[i]; 1230 } 1231 1232 void drft_forward(drft_lookup *l,float *data){ 1233 if(l->n==1)return; 1234 drftf1(l->n,data,l->trigcache,l->trigcache+l->n,l->splitcache); 1235 } 1236 1237 void drft_backward(drft_lookup *l,float *data){ 1238 if (l->n==1)return; 1239 drftb1(l->n,data,l->trigcache,l->trigcache+l->n,l->splitcache); 1240 } 1241 1242 void drft_init(drft_lookup *l,int n){ 1243 l->n=n; 1244 l->trigcache=_ogg_calloc(3*n,sizeof(*l->trigcache)); 1245 l->splitcache=_ogg_calloc(32,sizeof(*l->splitcache)); 1246 fdrffti(n, l->trigcache, l->splitcache); 1247 } 1248 1249 void drft_clear(drft_lookup *l){ 1250 if(l){ 1251 if(l->trigcache)_ogg_free(l->trigcache); 1252 if(l->splitcache)_ogg_free(l->splitcache); 1253 memset(l,0,sizeof(*l)); 1254 } 1255 } 1256