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