doom3-gpl
Doom 3 GPL source release
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
smallft.c
Go to the documentation of this file.
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-2002 *
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.17 2002/07/11 06:40:50 xiphmont 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 
1232  if(l->n==1)return;
1233  drftf1(l->n,data,l->trigcache,l->trigcache+l->n,l->splitcache);
1234 }
1235 
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(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 
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 }
int * splitcache
Definition: smallft.h:26
GLenum GLsizei n
Definition: glext.h:3705
case const float
Definition: Callbacks.cpp:62
#define _ogg_free
Definition: os_types.h:41
void drft_init(drft_lookup *l, int n)
Definition: smallft.c:1241
int i
Definition: process.py:33
list l
Definition: prepare.py:17
const GLubyte * c
Definition: glext.h:4677
void drft_backward(drft_lookup *l, float *data)
Definition: smallft.c:1236
GLsizei GLsizei GLenum GLenum const GLvoid * data
Definition: glext.h:2853
void drft_clear(drft_lookup *l)
Definition: smallft.c:1248
#define _ogg_calloc
Definition: os_types.h:39
GLint j
Definition: qgl.h:264
float * trigcache
Definition: smallft.h:25
void drft_forward(drft_lookup *l, float *data)
Definition: smallft.c:1231