lib/smallft.c

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