libvorbis-1.0.1/lib/smallft.c

   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,v 1.17 2002/07/11 06:40:50 xiphmont Exp $
  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 }