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-2009             *
9  * by the Xiph.Org Foundation http://www.xiph.org/                  *
10  *                                                                  *
11  ********************************************************************
12 
13  function: *unnormalized* fft transform
14  last mod: $Id: smallft.c 16227 2009-07-08 06:58:46Z xiphmont $
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 "os.h"
36 #include "misc.h"
37 
drfti1(int n,float * wa,int * ifac)38 static void drfti1(int n, float *wa, int *ifac){
39   static int ntryh[4] = { 4,2,3,5 };
40   static float tpi = 6.28318530717958648f;
41   float 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=(float)ld*argh;
95       fi=0.f;
96       for (ii=2;ii<ido;ii+=2){
97         fi+=1.f;
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 
fdrffti(int n,float * wsave,int * ifac)108 static void fdrffti(int n, float *wsave, int *ifac){
109 
110   if (n == 1) return;
111   drfti1(n, wsave+n, ifac);
112 }
113 
dradf2(int ido,int l1,float * cc,float * ch,float * wa1)114 static void dradf2(int ido,int l1,float *cc,float *ch,float *wa1){
115   int i,k;
116   float 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 
dradf4(int ido,int l1,float * cc,float * ch,float * wa1,float * wa2,float * wa3)169 static void dradf4(int ido,int l1,float *cc,float *ch,float *wa1,
170             float *wa2,float *wa3){
171   static float hsqt2 = .70710678118654752f;
172   int i,k,t0,t1,t2,t3,t4,t5,t6;
173   float 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 
dradfg(int ido,int ip,int l1,int idl1,float * cc,float * c1,float * c2,float * ch,float * ch2,float * wa)271 static void dradfg(int ido,int ip,int l1,int idl1,float *cc,float *c1,
272                           float *c2,float *ch,float *ch2,float *wa){
273 
274   static float tpi=6.283185307179586f;
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   float dc2,ai1,ai2,ar1,ar2,ds2;
278   int nbd;
279   float dcp,arg,dsp,ar1h,ar2h;
280   int idp2,ipp2;
281 
282   arg=tpi/(float)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.f;
409   ai1=0.f;
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 
drftf1(int n,float * c,float * ch,float * wa,int * ifac)573 static void drftf1(int n,float *c,float *ch,float *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 
dradb2(int ido,int l1,float * cc,float * ch,float * wa1)634 static void dradb2(int ido,int l1,float *cc,float *ch,float *wa1){
635   int i,k,t0,t1,t2,t3,t4,t5,t6;
636   float 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 
dradb3(int ido,int l1,float * cc,float * ch,float * wa1,float * wa2)686 static void dradb3(int ido,int l1,float *cc,float *ch,float *wa1,
687                           float *wa2){
688   static float taur = -.5f;
689   static float taui = .8660254037844386f;
690   int i,k,t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10;
691   float 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 
dradb4(int ido,int l1,float * cc,float * ch,float * wa1,float * wa2,float * wa3)749 static void dradb4(int ido,int l1,float *cc,float *ch,float *wa1,
750                           float *wa2,float *wa3){
751   static float sqrt2=1.414213562373095f;
752   int i,k,t0,t1,t2,t3,t4,t5,t6,t7,t8;
753   float 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 
dradbg(int ido,int ip,int l1,int idl1,float * cc,float * c1,float * c2,float * ch,float * ch2,float * wa)840 static void dradbg(int ido,int ip,int l1,int idl1,float *cc,float *c1,
841             float *c2,float *ch,float *ch2,float *wa){
842   static float tpi=6.283185307179586f;
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   float dc2,ai1,ai2,ar1,ar2,ds2;
846   int nbd;
847   float dcp,arg,dsp,ar1h,ar2h;
848   int ipp2;
849 
850   t10=ip*ido;
851   t0=l1*ido;
852   arg=tpi/(float)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.f;
980   ai1=0.f;
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 
drftb1(int n,float * c,float * ch,float * wa,int * ifac)1154 static void drftb1(int n, float *c, float *ch, float *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 
drft_forward(drft_lookup * l,float * data)1232 void drft_forward(drft_lookup *l,float *data){
1233   if(l->n==1)return;
1234   drftf1(l->n,data,l->trigcache,l->trigcache+l->n,l->splitcache);
1235 }
1236 
drft_backward(drft_lookup * l,float * data)1237 void drft_backward(drft_lookup *l,float *data){
1238   if (l->n==1)return;
1239   drftb1(l->n,data,l->trigcache,l->trigcache+l->n,l->splitcache);
1240 }
1241 
drft_init(drft_lookup * l,int n)1242 void drft_init(drft_lookup *l,int n){
1243   l->n=n;
1244   l->trigcache=_ogg_calloc(3*n,sizeof(*l->trigcache));
1245   l->splitcache=_ogg_calloc(32,sizeof(*l->splitcache));
1246   fdrffti(n, l->trigcache, l->splitcache);
1247 }
1248 
drft_clear(drft_lookup * l)1249 void drft_clear(drft_lookup *l){
1250   if(l){
1251     if(l->trigcache)_ogg_free(l->trigcache);
1252     if(l->splitcache)_ogg_free(l->splitcache);
1253     memset(l,0,sizeof(*l));
1254   }
1255 }
1256