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