lib/smallft.c

   1 /********************************************************************
   2  *                                                                  *
   3  * THIS FILE IS PART OF THE Ogg Vorbis SOFTWARE CODEC SOURCE CODE.  *
   4  * USE, DISTRIBUTION AND REPRODUCTION OF THIS SOURCE IS GOVERNED BY *
   5  * THE GNU PUBLIC LICENSE 2, WHICH IS INCLUDED WITH THIS SOURCE.    *
   6  * PLEASE READ THESE TERMS DISTRIBUTING.                            *
   7  *                                                                  *
   8  * THE OggSQUISH SOURCE CODE IS (C) COPYRIGHT 1994-2000             *
   9  * by Monty <monty@xiph.org> and The XIPHOPHORUS Company            *
  10  * http://www.xiph.org/                                             *
  11  *                                                                  *
  12  ********************************************************************
  13
  14  function: *unnormalized* fft transform
  15  last mod: $Id: smallft.c,v 1.8 2000/03/10 13:21:18 xiphmont Exp $
  16
  17 ********************************************************************/
  18
  19 /* FFT implementation from OggSquish, minus cosine transforms,
  20  * minus all but radix 2/4 case.  In Vorbis we only need this
  21  * cut-down version.
  22  *
  23  * To do more than just power-of-two sized vectors, see the full
  24  * version I wrote for NetLib.
  25  *
  26  * Note that the packing is a little strange; rather than the FFT r/i
  27  * packing following R_0, I_n, R_1, I_1, R_2, I_2 ... R_n-1, I_n-1,
  28  * it follows R_0, R_1, I_1, R_2, I_2 ... R_n-1, I_n-1, I_n like the
  29  * FORTRAN version
  30  */
  31
  32 #include <stdlib.h>
  33 #include <string.h>
  34 #include <math.h>
  35 #include "smallft.h"
  36 #include "misc.h"
  37
  38 static void drfti1(int n, double *wa, int *ifac){
  39   static int ntryh[4] = { 4,2,3,5 };
  40   static double tpi = 6.28318530717958647692528676655900577;
  41   double 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=(double)ld*argh;
  95       fi=0.;
  96       for (ii=2;ii<ido;ii+=2){
  97         fi+=1.;
  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, double *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,double *cc,double *ch,double *wa1){
 115   int i,k;
 116   double 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,double *cc,double *ch,double *wa1,
 170             double *wa2,double *wa3){
 171   static double hsqt2 = .70710678118654752440084436210485;
 172   int i,k,t0,t1,t2,t3,t4,t5,t6;
 173   double 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,double *cc,double *c1,
 272                           double *c2,double *ch,double *ch2,double *wa){
 273
 274   static double tpi=6.28318530717958647692528676655900577;
 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   double dc2,ai1,ai2,ar1,ar2,ds2;
 278   int nbd;
 279   double dcp,arg,dsp,ar1h,ar2h;
 280   int idp2,ipp2;
 281
 282   arg=tpi/(double)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.;
 409   ai1=0.;
 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,double *c,double *ch,double *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,double *cc,double *ch,double *wa1){
 635   int i,k,t0,t1,t2,t3,t4,t5,t6;
 636   double 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,double *cc,double *ch,double *wa1,
 687                           double *wa2){
 688   static double taur = -.5;
 689   static double taui = .86602540378443864676372317075293618;
 690   int i,k,t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10;
 691   double 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,double *cc,double *ch,double *wa1,
 750                           double *wa2,double *wa3){
 751   static double sqrt2=1.4142135623730950488016887242097;
 752   int i,k,t0,t1,t2,t3,t4,t5,t6,t7,t8;
 753   double 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,double *cc,double *c1,
 841             double *c2,double *ch,double *ch2,double *wa){
 842   static double tpi=6.28318530717958647692528676655900577;
 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   double dc2,ai1,ai2,ar1,ar2,ds2;
 846   int nbd;
 847   double dcp,arg,dsp,ar1h,ar2h;
 848   int ipp2;
 849
 850   t10=ip*ido;
 851   t0=l1*ido;
 852   arg=tpi/(double)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.;
 980   ai1=0.;
 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, double *c, double *ch, double *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,double *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,double *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=calloc(3*n,sizeof(double));
1245   l->splitcache=calloc(32,sizeof(int));
1246   fdrffti(n, l->trigcache, l->splitcache);
1247 }
1248
1249 void drft_clear(drft_lookup *l){
1250   if(l){
1251     if(l->trigcache)free(l->trigcache);
1252     if(l->splitcache)free(l->splitcache);
1253     memset(l,0,sizeof(drft_lookup));
1254   }
1255 }