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 https://xiph.org/ *
10 * *
11 ********************************************************************
12
13 function: *unnormalized* fft transform
14
15 ********************************************************************/
16
17/* FFT implementation from OggSquish, minus cosine transforms,
18 * minus all but radix 2/4 case. In Vorbis we only need this
19 * cut-down version.
20 *
21 * To do more than just power-of-two sized vectors, see the full
22 * version I wrote for NetLib.
23 *
24 * Note that the packing is a little strange; rather than the FFT r/i
25 * packing following R_0, I_n, R_1, I_1, R_2, I_2 ... R_n-1, I_n-1,
26 * it follows R_0, R_1, I_1, R_2, I_2 ... R_n-1, I_n-1, I_n like the
27 * FORTRAN version
28 */
29
30#include <stdlib.h>
31#include <string.h>
32#include <math.h>
33#include "smallft.h"
34#include "os.h"
35#include "misc.h"
36
37static 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
107static void fdrffti(int n, float *wsave, int *ifac){
108
109 if (n == 1) return;
110 drfti1(n, wsave+n, ifac);
111}
112
113static 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
168static 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
270static 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
389L119:
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
572static 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
633static 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
674L105:
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
685static 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
748static 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
839static 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
977L116:
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
1096L132:
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
1153static 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
1231void drft_forward(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
1236void drft_backward(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
1241void drft_init(drft_lookup *l,int n){
1242 l->n=n;
1243 l->trigcache=_ogg_calloc(3*n,sizeof(*l->trigcache));
1244 l->splitcache=_ogg_calloc(32,sizeof(*l->splitcache));
1245 fdrffti(n, l->trigcache, l->splitcache);
1246}
1247
1248void drft_clear(drft_lookup *l){
1249 if(l){
1250 if(l->trigcache)_ogg_free(l->trigcache);
1251 if(l->splitcache)_ogg_free(l->splitcache);
1252 memset(l,0,sizeof(*l));
1253 }
1254}
1255