Home | History | Annotate | Download | only in lib 
      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