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