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