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