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