The Tor Browser / file comparison

comparison: media/libvorbis/lib/vorbis_smallft.c

media/libvorbis/lib/vorbis_smallft.c

branch
TOR_BUG_9701
changeset 8
97036ab72558
equal deleted inserted replaced
-1:000000000000 0:901c1a8b5b76
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
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
108 static void fdrffti(int n, float *wsave, int *ifac){
109
110 if (n == 1) return;
111 drfti1(n, wsave+n, ifac);
112 }
113
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
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
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
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
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
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
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
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
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
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
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
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
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 }

Mercurial Repository: The Tor Browser

Europalab SCM Repository Server

mercurial