apps/codecs/libspeex/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$
  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-speex.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 }