]> Joshua Wise's Git repositories - patchfork.git/blame - jorbis/src/com/jcraft/jorbis/Drft.java
initial import from pitchfork-0.5.5
[patchfork.git] / jorbis / src / com / jcraft / jorbis / Drft.java
CommitLineData
964dd0bc
JW
1/* JOrbis
2 * Copyright (C) 2000 ymnk, JCraft,Inc.
3 *
4 * Written by: 2000 ymnk<ymnk@jcraft.com>
5 *
6 * Many thanks to
7 * Monty <monty@xiph.org> and
8 * The XIPHOPHORUS Company http://www.xiph.org/ .
9 * JOrbis has been based on their awesome works, Vorbis codec.
10 *
11 * This program is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Library General Public License
13 * as published by the Free Software Foundation; either version 2 of
14 * the License, or (at your option) any later version.
15
16 * This program is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 * GNU Library General Public License for more details.
20 *
21 * You should have received a copy of the GNU Library General Public
22 * License along with this program; if not, write to the Free Software
23 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
24 */
25
26package com.jcraft.jorbis;
27
28class Drft{
29 int n;
30 float[] trigcache;
31 int[] splitcache;
32
33 void backward(float[] data){
34 //System.err.println("Drft.backward");
35 if(n==1)return;
36 drftb1(n,data,trigcache,trigcache,n,splitcache);
37 }
38
39 void init(int n){
40 //System.err.println("Drft.init");
41 this.n=n;
42 trigcache=new float[3*n];
43 splitcache=new int[32];
44 fdrffti(n, trigcache, splitcache);
45 }
46
47 void clear(){
48 //System.err.println("Drft.clear");
49 if(trigcache!=null)trigcache=null;
50 if(splitcache!=null)splitcache=null;
51// memset(l,0,sizeof(drft_lookup));
52 }
53
54 static int[] ntryh = { 4,2,3,5 };
55 static float tpi = 6.28318530717958647692528676655900577f;
56 static float hsqt2 = .70710678118654752440084436210485f;
57 static float taui = .86602540378443864676372317075293618f;
58 static float taur = -.5f;
59 static float sqrt2 = 1.4142135623730950488016887242097f;
60
61 static void drfti1(int n, float[] wa, int index, int[] ifac){
62 float arg,argh,argld,fi;
63 int ntry=0,i,j=-1;
64 int k1, l1, l2, ib;
65 int ld, ii, ip, is, nq, nr;
66 int ido, ipm, nfm1;
67 int nl=n;
68 int nf=0;
69
70 int state=101;
71
72 loop: while(true){
73 switch(state){
74 case 101:
75 j++;
76 if (j < 4)
77 ntry=ntryh[j];
78 else
79 ntry+=2;
80 case 104:
81 nq=nl/ntry;
82 nr=nl-ntry*nq;
83 if(nr!=0){
84 state=101;
85 break;
86 }
87 nf++;
88 ifac[nf+1]=ntry;
89 nl=nq;
90 if(ntry!=2){
91 state=107;
92 break;
93 }
94 if(nf==1){
95 state=107;
96 break;
97 }
98
99 for(i=1;i<nf;i++){
100 ib=nf-i+1;
101 ifac[ib+1]=ifac[ib];
102 }
103 ifac[2] = 2;
104 case 107:
105 if(nl!=1){
106 state=104;
107 break;
108 }
109 ifac[0]=n;
110 ifac[1]=nf;
111 argh=tpi/n;
112 is=0;
113 nfm1=nf-1;
114 l1=1;
115
116 if(nfm1==0)return;
117
118 for (k1=0;k1<nfm1;k1++){
119 ip=ifac[k1+2];
120 ld=0;
121 l2=l1*ip;
122 ido=n/l2;
123 ipm=ip-1;
124
125 for (j=0;j<ipm;j++){
126 ld+=l1;
127 i=is;
128 argld=(float)ld*argh;
129 fi=0.f;
130 for (ii=2;ii<ido;ii+=2){
131 fi+=1.f;
132 arg=fi*argld;
133 wa[index+i++]=(float)Math.cos(arg);
134 wa[index+i++]=(float)Math.sin(arg);
135 }
136 is+=ido;
137 }
138 l1=l2;
139 }
140 break loop;
141 }
142 }
143 }
144
145 static void fdrffti(int n, float[] wsave, int[] ifac){
146//System.err.println("fdrffti: n="+n);
147 if(n == 1) return;
148 drfti1(n, wsave, n, ifac);
149 }
150
151 static void dradf2(int ido,int l1,float[] cc, float[] ch, float[] wa1, int index){
152 int i,k;
153 float ti2,tr2;
154 int t0,t1,t2,t3,t4,t5,t6;
155
156 t1=0;
157 t0=(t2=l1*ido);
158 t3=ido<<1;
159 for(k=0;k<l1;k++){
160 ch[t1<<1]=cc[t1]+cc[t2];
161 ch[(t1<<1)+t3-1]=cc[t1]-cc[t2];
162 t1+=ido;
163 t2+=ido;
164 }
165
166 if(ido<2)return;
167
168 if(ido!=2){
169 t1=0;
170 t2=t0;
171 for(k=0;k<l1;k++){
172 t3=t2;
173 t4=(t1<<1)+(ido<<1);
174 t5=t1;
175 t6=t1+t1;
176 for(i=2;i<ido;i+=2){
177 t3+=2;
178 t4-=2;
179 t5+=2;
180 t6+=2;
181 tr2=wa1[index+i-2]*cc[t3-1]+wa1[index+i-1]*cc[t3];
182 ti2=wa1[index+i-2]*cc[t3]-wa1[index+i-1]*cc[t3-1];
183 ch[t6]=cc[t5]+ti2;
184 ch[t4]=ti2-cc[t5];
185 ch[t6-1]=cc[t5-1]+tr2;
186 ch[t4-1]=cc[t5-1]-tr2;
187 }
188 t1+=ido;
189 t2+=ido;
190 }
191 if(ido%2==1)return;
192 }
193
194 t3=(t2=(t1=ido)-1);
195 t2+=t0;
196 for(k=0;k<l1;k++){
197 ch[t1]=-cc[t2];
198 ch[t1-1]=cc[t3];
199 t1+=ido<<1;
200 t2+=ido;
201 t3+=ido;
202 }
203 }
204
205 static void dradf4(int ido,int l1,float[] cc, float[] ch,
206 float[] wa1, int index1,
207 float[] wa2, int index2,
208 float[] wa3, int index3){
209 int i,k,t0,t1,t2,t3,t4,t5,t6;
210 float ci2,ci3,ci4,cr2,cr3,cr4,ti1,ti2,ti3,ti4,tr1,tr2,tr3,tr4;
211 t0=l1*ido;
212
213 t1=t0;
214 t4=t1<<1;
215 t2=t1+(t1<<1);
216 t3=0;
217
218 for(k=0;k<l1;k++){
219 tr1=cc[t1]+cc[t2];
220 tr2=cc[t3]+cc[t4];
221
222 ch[t5=t3<<2]=tr1+tr2;
223 ch[(ido<<2)+t5-1]=tr2-tr1;
224 ch[(t5+=(ido<<1))-1]=cc[t3]-cc[t4];
225 ch[t5]=cc[t2]-cc[t1];
226
227 t1+=ido;
228 t2+=ido;
229 t3+=ido;
230 t4+=ido;
231 }
232 if(ido<2)return;
233
234 if(ido!=2){
235 t1=0;
236 for(k=0;k<l1;k++){
237 t2=t1;
238 t4=t1<<2;
239 t5=(t6=ido<<1)+t4;
240 for(i=2;i<ido;i+=2){
241 t3=(t2+=2);
242 t4+=2;
243 t5-=2;
244
245 t3+=t0;
246 cr2=wa1[index1+i-2]*cc[t3-1]+wa1[index1+i-1]*cc[t3];
247 ci2=wa1[index1+i-2]*cc[t3]-wa1[index1+i-1]*cc[t3-1];
248 t3+=t0;
249 cr3=wa2[index2+i-2]*cc[t3-1]+wa2[index2+i-1]*cc[t3];
250 ci3=wa2[index2+i-2]*cc[t3]-wa2[index2+i-1]*cc[t3-1];
251 t3+=t0;
252 cr4=wa3[index3+i-2]*cc[t3-1]+wa3[index3+i-1]*cc[t3];
253 ci4=wa3[index3+i-2]*cc[t3]-wa3[index3+i-1]*cc[t3-1];
254
255 tr1=cr2+cr4;
256 tr4=cr4-cr2;
257 ti1=ci2+ci4;
258 ti4=ci2-ci4;
259
260 ti2=cc[t2]+ci3;
261 ti3=cc[t2]-ci3;
262 tr2=cc[t2-1]+cr3;
263 tr3=cc[t2-1]-cr3;
264
265 ch[t4-1]=tr1+tr2;
266 ch[t4]=ti1+ti2;
267
268 ch[t5-1]=tr3-ti4;
269 ch[t5]=tr4-ti3;
270
271 ch[t4+t6-1]=ti4+tr3;
272 ch[t4+t6]=tr4+ti3;
273
274 ch[t5+t6-1]=tr2-tr1;
275 ch[t5+t6]=ti1-ti2;
276 }
277 t1+=ido;
278 }
279 if((ido&1)!=0)return;
280 }
281
282 t2=(t1=t0+ido-1)+(t0<<1);
283 t3=ido<<2;
284 t4=ido;
285 t5=ido<<1;
286 t6=ido;
287
288 for(k=0;k<l1;k++){
289 ti1=-hsqt2*(cc[t1]+cc[t2]);
290 tr1=hsqt2*(cc[t1]-cc[t2]);
291
292 ch[t4-1]=tr1+cc[t6-1];
293 ch[t4+t5-1]=cc[t6-1]-tr1;
294
295 ch[t4]=ti1-cc[t1+t0];
296 ch[t4+t5]=ti1+cc[t1+t0];
297
298 t1+=ido;
299 t2+=ido;
300 t4+=t3;
301 t6+=ido;
302 }
303 }
304
305 static void dradfg(int ido,int ip,int l1,int idl1,float[] cc,float[] c1,
306 float[] c2, float[] ch, float[] ch2, float[] wa, int index){
307 int idij,ipph,i,j,k,l,ic,ik,is;
308 int t0,t1,t2=0,t3,t4,t5,t6,t7,t8,t9,t10;
309 float dc2,ai1,ai2,ar1,ar2,ds2;
310 int nbd;
311 float dcp=0,arg,dsp=0,ar1h,ar2h;
312 int idp2,ipp2;
313
314 arg=tpi/(float)ip;
315 dcp=(float)Math.cos(arg);
316 dsp=(float)Math.sin(arg);
317 ipph=(ip+1)>>1;
318 ipp2=ip;
319 idp2=ido;
320 nbd=(ido-1)>>1;
321 t0=l1*ido;
322 t10=ip*ido;
323
324 int state=100;
325 loop: while(true){
326 switch(state){
327 case 101:
328 if(ido==1){
329 state=119;
330 break;
331 }
332 for(ik=0;ik<idl1;ik++)ch2[ik]=c2[ik];
333
334 t1=0;
335 for(j=1;j<ip;j++){
336 t1+=t0;
337 t2=t1;
338 for(k=0;k<l1;k++){
339 ch[t2]=c1[t2];
340 t2+=ido;
341 }
342 }
343
344 is=-ido;
345 t1=0;
346 if(nbd>l1){
347 for(j=1;j<ip;j++){
348 t1+=t0;
349 is+=ido;
350 t2= -ido+t1;
351 for(k=0;k<l1;k++){
352 idij=is-1;
353 t2+=ido;
354 t3=t2;
355 for(i=2;i<ido;i+=2){
356 idij+=2;
357 t3+=2;
358 ch[t3-1]=wa[index+idij-1]*c1[t3-1]+wa[index+idij]*c1[t3];
359 ch[t3]=wa[index+idij-1]*c1[t3]-wa[index+idij]*c1[t3-1];
360 }
361 }
362 }
363 }
364 else{
365
366 for(j=1;j<ip;j++){
367 is+=ido;
368 idij=is-1;
369 t1+=t0;
370 t2=t1;
371 for(i=2;i<ido;i+=2){
372 idij+=2;
373 t2+=2;
374 t3=t2;
375 for(k=0;k<l1;k++){
376 ch[t3-1]=wa[index+idij-1]*c1[t3-1]+wa[index+idij]*c1[t3];
377 ch[t3]=wa[index+idij-1]*c1[t3]-wa[index+idij]*c1[t3-1];
378 t3+=ido;
379 }
380 }
381 }
382 }
383
384 t1=0;
385 t2=ipp2*t0;
386 if(nbd<l1){
387 for(j=1;j<ipph;j++){
388 t1+=t0;
389 t2-=t0;
390 t3=t1;
391 t4=t2;
392 for(i=2;i<ido;i+=2){
393 t3+=2;
394 t4+=2;
395 t5=t3-ido;
396 t6=t4-ido;
397 for(k=0;k<l1;k++){
398 t5+=ido;
399 t6+=ido;
400 c1[t5-1]=ch[t5-1]+ch[t6-1];
401 c1[t6-1]=ch[t5]-ch[t6];
402 c1[t5]=ch[t5]+ch[t6];
403 c1[t6]=ch[t6-1]-ch[t5-1];
404 }
405 }
406 }
407 }
408 else{
409 for(j=1;j<ipph;j++){
410 t1+=t0;
411 t2-=t0;
412 t3=t1;
413 t4=t2;
414 for(k=0;k<l1;k++){
415 t5=t3;
416 t6=t4;
417 for(i=2;i<ido;i+=2){
418 t5+=2;
419 t6+=2;
420 c1[t5-1]=ch[t5-1]+ch[t6-1];
421 c1[t6-1]=ch[t5]-ch[t6];
422 c1[t5]=ch[t5]+ch[t6];
423 c1[t6]=ch[t6-1]-ch[t5-1];
424 }
425 t3+=ido;
426 t4+=ido;
427 }
428 }
429 }
430 case 119:
431 for(ik=0;ik<idl1;ik++)c2[ik]=ch2[ik];
432
433 t1=0;
434 t2=ipp2*idl1;
435 for(j=1;j<ipph;j++){
436 t1+=t0;
437 t2-=t0;
438 t3=t1-ido;
439 t4=t2-ido;
440 for(k=0;k<l1;k++){
441 t3+=ido;
442 t4+=ido;
443 c1[t3]=ch[t3]+ch[t4];
444 c1[t4]=ch[t4]-ch[t3];
445 }
446 }
447
448 ar1=1.f;
449 ai1=0.f;
450 t1=0;
451 t2=ipp2*idl1;
452 t3=(ip-1)*idl1;
453 for(l=1;l<ipph;l++){
454 t1+=idl1;
455 t2-=idl1;
456 ar1h=dcp*ar1-dsp*ai1;
457 ai1=dcp*ai1+dsp*ar1;
458 ar1=ar1h;
459 t4=t1;
460 t5=t2;
461 t6=t3;
462 t7=idl1;
463
464 for(ik=0;ik<idl1;ik++){
465 ch2[t4++]=c2[ik]+ar1*c2[t7++];
466 ch2[t5++]=ai1*c2[t6++];
467 }
468
469 dc2=ar1;
470 ds2=ai1;
471 ar2=ar1;
472 ai2=ai1;
473
474 t4=idl1;
475 t5=(ipp2-1)*idl1;
476 for(j=2;j<ipph;j++){
477 t4+=idl1;
478 t5-=idl1;
479
480 ar2h=dc2*ar2-ds2*ai2;
481 ai2=dc2*ai2+ds2*ar2;
482 ar2=ar2h;
483
484 t6=t1;
485 t7=t2;
486 t8=t4;
487 t9=t5;
488 for(ik=0;ik<idl1;ik++){
489 ch2[t6++]+=ar2*c2[t8++];
490 ch2[t7++]+=ai2*c2[t9++];
491 }
492 }
493 }
494 t1=0;
495 for(j=1;j<ipph;j++){
496 t1+=idl1;
497 t2=t1;
498 for(ik=0;ik<idl1;ik++)ch2[ik]+=c2[t2++];
499 }
500
501 if(ido<l1){
502 state=132;
503 break;
504 }
505
506 t1=0;
507 t2=0;
508 for(k=0;k<l1;k++){
509 t3=t1;
510 t4=t2;
511 for(i=0;i<ido;i++)cc[t4++]=ch[t3++];
512 t1+=ido;
513 t2+=t10;
514 }
515 state=135;
516 break;
517
518 case 132:
519 for(i=0;i<ido;i++){
520 t1=i;
521 t2=i;
522 for(k=0;k<l1;k++){
523 cc[t2]=ch[t1];
524 t1+=ido;
525 t2+=t10;
526 }
527 }
528 case 135:
529 t1=0;
530 t2=ido<<1;
531 t3=0;
532 t4=ipp2*t0;
533 for(j=1;j<ipph;j++){
534 t1+=t2;
535 t3+=t0;
536 t4-=t0;
537
538 t5=t1;
539 t6=t3;
540 t7=t4;
541
542 for(k=0;k<l1;k++){
543 cc[t5-1]=ch[t6];
544 cc[t5]=ch[t7];
545 t5+=t10;
546 t6+=ido;
547 t7+=ido;
548 }
549 }
550
551 if(ido==1)return;
552 if(nbd<l1){
553 state=141;
554 break;
555 }
556
557 t1=-ido;
558 t3=0;
559 t4=0;
560 t5=ipp2*t0;
561 for(j=1;j<ipph;j++){
562 t1+=t2;
563 t3+=t2;
564 t4+=t0;
565 t5-=t0;
566 t6=t1;
567 t7=t3;
568 t8=t4;
569 t9=t5;
570 for(k=0;k<l1;k++){
571 for(i=2;i<ido;i+=2){
572 ic=idp2-i;
573 cc[i+t7-1]=ch[i+t8-1]+ch[i+t9-1];
574 cc[ic+t6-1]=ch[i+t8-1]-ch[i+t9-1];
575 cc[i+t7]=ch[i+t8]+ch[i+t9];
576 cc[ic+t6]=ch[i+t9]-ch[i+t8];
577 }
578 t6+=t10;
579 t7+=t10;
580 t8+=ido;
581 t9+=ido;
582 }
583 }
584 return;
585 case 141:
586 t1=-ido;
587 t3=0;
588 t4=0;
589 t5=ipp2*t0;
590 for(j=1;j<ipph;j++){
591 t1+=t2;
592 t3+=t2;
593 t4+=t0;
594 t5-=t0;
595 for(i=2;i<ido;i+=2){
596 t6=idp2+t1-i;
597 t7=i+t3;
598 t8=i+t4;
599 t9=i+t5;
600 for(k=0;k<l1;k++){
601 cc[t7-1]=ch[t8-1]+ch[t9-1];
602 cc[t6-1]=ch[t8-1]-ch[t9-1];
603 cc[t7]=ch[t8]+ch[t9];
604 cc[t6]=ch[t9]-ch[t8];
605 t6+=t10;
606 t7+=t10;
607 t8+=ido;
608 t9+=ido;
609 }
610 }
611 }
612 break loop;
613 }
614 }
615 }
616
617 static void drftf1(int n,float[] c, float[] ch, float[] wa, int[] ifac){
618 int i,k1,l1,l2;
619 int na,kh,nf;
620 int ip,iw,ido,idl1,ix2,ix3;
621
622 nf=ifac[1];
623 na=1;
624 l2=n;
625 iw=n;
626
627 for(k1=0;k1<nf;k1++){
628 kh=nf-k1;
629 ip=ifac[kh+1];
630 l1=l2/ip;
631 ido=n/l2;
632 idl1=ido*l1;
633 iw-=(ip-1)*ido;
634 na=1-na;
635
636 int state=100;
637 loop: while(true){
638 switch(state){
639 case 100:
640 if(ip!=4){
641 state=102;
642 break;
643 }
644
645 ix2=iw+ido;
646 ix3=ix2+ido;
647 if(na!=0)
648 dradf4(ido,l1,ch,c,wa,iw-1,wa,ix2-1,wa,ix3-1);
649 else
650 dradf4(ido,l1,c,ch,wa,iw-1,wa,ix2-1,wa,ix3-1);
651 state=110;
652 break;
653 case 102:
654 if(ip!=2){
655 state=104;
656 break;
657 }
658 if(na!=0){
659 state=103;
660 break;
661 }
662 dradf2(ido,l1,c,ch,wa, iw-1);
663 state=110;
664 break;
665 case 103:
666 dradf2(ido,l1,ch,c,wa, iw-1);
667 case 104:
668 if(ido==1)na=1-na;
669 if(na!=0){
670 state=109;
671 break;
672 }
673 dradfg(ido,ip,l1,idl1,c,c,c,ch,ch,wa,iw-1);
674 na=1;
675 state=110;
676 break;
677 case 109:
678 dradfg(ido,ip,l1,idl1,ch,ch,ch,c,c,wa,iw-1);
679 na=0;
680 case 110:
681 l2=l1;
682 break loop;
683 }
684 }
685 }
686 if(na==1)return;
687 for(i=0;i<n;i++)c[i]=ch[i];
688 }
689
690 static void dradb2(int ido,int l1,float[] cc,float[] ch,float[] wa1, int index){
691 int i,k,t0,t1,t2,t3,t4,t5,t6;
692 float ti2,tr2;
693
694 t0=l1*ido;
695
696 t1=0;
697 t2=0;
698 t3=(ido<<1)-1;
699 for(k=0;k<l1;k++){
700 ch[t1]=cc[t2]+cc[t3+t2];
701 ch[t1+t0]=cc[t2]-cc[t3+t2];
702 t2=(t1+=ido)<<1;
703 }
704
705 if(ido<2)return;
706 if(ido!=2){
707 t1=0;
708 t2=0;
709 for(k=0;k<l1;k++){
710 t3=t1;
711 t5=(t4=t2)+(ido<<1);
712 t6=t0+t1;
713 for(i=2;i<ido;i+=2){
714 t3+=2;
715 t4+=2;
716 t5-=2;
717 t6+=2;
718 ch[t3-1]=cc[t4-1]+cc[t5-1];
719 tr2=cc[t4-1]-cc[t5-1];
720 ch[t3]=cc[t4]-cc[t5];
721 ti2=cc[t4]+cc[t5];
722 ch[t6-1]=wa1[index+i-2]*tr2-wa1[index+i-1]*ti2;
723 ch[t6]=wa1[index+i-2]*ti2+wa1[index+i-1]*tr2;
724 }
725 t2=(t1+=ido)<<1;
726 }
727 if((ido%2)==1)return;
728 }
729
730 t1=ido-1;
731 t2=ido-1;
732 for(k=0;k<l1;k++){
733 ch[t1]=cc[t2]+cc[t2];
734 ch[t1+t0]=-(cc[t2+1]+cc[t2+1]);
735 t1+=ido;
736 t2+=ido<<1;
737 }
738 }
739
740 static void dradb3(int ido,int l1,float[] cc,float[] ch,
741 float[] wa1, int index1,
742 float[] wa2, int index2){
743 int i,k,t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10;
744 float ci2,ci3,di2,di3,cr2,cr3,dr2,dr3,ti2,tr2;
745 t0=l1*ido;
746
747 t1=0;
748 t2=t0<<1;
749 t3=ido<<1;
750 t4=ido+(ido<<1);
751 t5=0;
752 for(k=0;k<l1;k++){
753 tr2=cc[t3-1]+cc[t3-1];
754 cr2=cc[t5]+(taur*tr2);
755 ch[t1]=cc[t5]+tr2;
756 ci3=taui*(cc[t3]+cc[t3]);
757 ch[t1+t0]=cr2-ci3;
758 ch[t1+t2]=cr2+ci3;
759 t1+=ido;
760 t3+=t4;
761 t5+=t4;
762 }
763
764 if(ido==1)return;
765
766 t1=0;
767 t3=ido<<1;
768 for(k=0;k<l1;k++){
769 t7=t1+(t1<<1);
770 t6=(t5=t7+t3);
771 t8=t1;
772 t10=(t9=t1+t0)+t0;
773
774 for(i=2;i<ido;i+=2){
775 t5+=2;
776 t6-=2;
777 t7+=2;
778 t8+=2;
779 t9+=2;
780 t10+=2;
781 tr2=cc[t5-1]+cc[t6-1];
782 cr2=cc[t7-1]+(taur*tr2);
783 ch[t8-1]=cc[t7-1]+tr2;
784 ti2=cc[t5]-cc[t6];
785 ci2=cc[t7]+(taur*ti2);
786 ch[t8]=cc[t7]+ti2;
787 cr3=taui*(cc[t5-1]-cc[t6-1]);
788 ci3=taui*(cc[t5]+cc[t6]);
789 dr2=cr2-ci3;
790 dr3=cr2+ci3;
791 di2=ci2+cr3;
792 di3=ci2-cr3;
793 ch[t9-1]=wa1[index1+i-2]*dr2-wa1[index1+i-1]*di2;
794 ch[t9]=wa1[index1+i-2]*di2+wa1[index1+i-1]*dr2;
795 ch[t10-1]=wa2[index2+i-2]*dr3-wa2[index2+i-1]*di3;
796 ch[t10]=wa2[index2+i-2]*di3+wa2[index2+i-1]*dr3;
797 }
798 t1+=ido;
799 }
800 }
801
802 static void dradb4(int ido,int l1,float[] cc,float[] ch,
803 float[] wa1, int index1,
804 float[] wa2, int index2,
805 float[] wa3, int index3){
806 int i,k,t0,t1,t2,t3,t4,t5,t6,t7,t8;
807 float ci2,ci3,ci4,cr2,cr3,cr4,ti1,ti2,ti3,ti4,tr1,tr2,tr3,tr4;
808 t0=l1*ido;
809
810 t1=0;
811 t2=ido<<2;
812 t3=0;
813 t6=ido<<1;
814 for(k=0;k<l1;k++){
815 t4=t3+t6;
816 t5=t1;
817 tr3=cc[t4-1]+cc[t4-1];
818 tr4=cc[t4]+cc[t4];
819 tr1=cc[t3]-cc[(t4+=t6)-1];
820 tr2=cc[t3]+cc[t4-1];
821 ch[t5]=tr2+tr3;
822 ch[t5+=t0]=tr1-tr4;
823 ch[t5+=t0]=tr2-tr3;
824 ch[t5+=t0]=tr1+tr4;
825 t1+=ido;
826 t3+=t2;
827 }
828
829 if(ido<2)return;
830 if(ido!=2){
831 t1=0;
832 for(k=0;k<l1;k++){
833 t5=(t4=(t3=(t2=t1<<2)+t6))+t6;
834 t7=t1;
835 for(i=2;i<ido;i+=2){
836 t2+=2;
837 t3+=2;
838 t4-=2;
839 t5-=2;
840 t7+=2;
841 ti1=cc[t2]+cc[t5];
842 ti2=cc[t2]-cc[t5];
843 ti3=cc[t3]-cc[t4];
844 tr4=cc[t3]+cc[t4];
845 tr1=cc[t2-1]-cc[t5-1];
846 tr2=cc[t2-1]+cc[t5-1];
847 ti4=cc[t3-1]-cc[t4-1];
848 tr3=cc[t3-1]+cc[t4-1];
849 ch[t7-1]=tr2+tr3;
850 cr3=tr2-tr3;
851 ch[t7]=ti2+ti3;
852 ci3=ti2-ti3;
853 cr2=tr1-tr4;
854 cr4=tr1+tr4;
855 ci2=ti1+ti4;
856 ci4=ti1-ti4;
857
858 ch[(t8=t7+t0)-1]=wa1[index1+i-2]*cr2-wa1[index1+i-1]*ci2;
859 ch[t8]=wa1[index1+i-2]*ci2+wa1[index1+i-1]*cr2;
860 ch[(t8+=t0)-1]=wa2[index2+i-2]*cr3-wa2[index2+i-1]*ci3;
861 ch[t8]=wa2[index2+i-2]*ci3+wa2[index2+i-1]*cr3;
862 ch[(t8+=t0)-1]=wa3[index3+i-2]*cr4-wa3[index3+i-1]*ci4;
863 ch[t8]=wa3[index3+i-2]*ci4+wa3[index3+i-1]*cr4;
864 }
865 t1+=ido;
866 }
867 if(ido%2 == 1)return;
868 }
869
870 t1=ido;
871 t2=ido<<2;
872 t3=ido-1;
873 t4=ido+(ido<<1);
874 for(k=0;k<l1;k++){
875 t5=t3;
876 ti1=cc[t1]+cc[t4];
877 ti2=cc[t4]-cc[t1];
878 tr1=cc[t1-1]-cc[t4-1];
879 tr2=cc[t1-1]+cc[t4-1];
880 ch[t5]=tr2+tr2;
881 ch[t5+=t0]=sqrt2*(tr1-ti1);
882 ch[t5+=t0]=ti2+ti2;
883 ch[t5+=t0]=-sqrt2*(tr1+ti1);
884
885 t3+=ido;
886 t1+=t2;
887 t4+=t2;
888 }
889 }
890
891 static void dradbg(int ido,int ip,int l1,int idl1,float[] cc,float[] c1,
892 float[] c2,float[] ch,float[] ch2,float[] wa, int index ){
893
894 int idij,ipph=0,i,j,k,l,ik,is,t0=0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10=0,
895 t11,t12;
896 float dc2,ai1,ai2,ar1,ar2,ds2;
897 int nbd=0;
898 float dcp=0,arg,dsp=0,ar1h,ar2h;
899 int ipp2=0;
900
901 int state=100;
902
903 loop: while(true){
904 switch(state){
905 case 100:
906 t10=ip*ido;
907 t0=l1*ido;
908 arg=tpi/(float)ip;
909 dcp=(float)Math.cos(arg);
910 dsp=(float)Math.sin(arg);
911 nbd=(ido-1)>>>1;
912 ipp2=ip;
913 ipph=(ip+1)>>>1;
914 if(ido<l1){
915 state=103;
916 break;
917 }
918 t1=0;
919 t2=0;
920 for(k=0;k<l1;k++){
921 t3=t1;
922 t4=t2;
923 for(i=0;i<ido;i++){
924 ch[t3]=cc[t4];
925 t3++;
926 t4++;
927 }
928 t1+=ido;
929 t2+=t10;
930 }
931 state=106;
932 break;
933 case 103:
934 t1=0;
935 for(i=0;i<ido;i++){
936 t2=t1;
937 t3=t1;
938 for(k=0;k<l1;k++){
939 ch[t2]=cc[t3];
940 t2+=ido;
941 t3+=t10;
942 }
943 t1++;
944 }
945 case 106:
946 t1=0;
947 t2=ipp2*t0;
948 t7=(t5=ido<<1);
949 for(j=1;j<ipph;j++){
950 t1+=t0;
951 t2-=t0;
952 t3=t1;
953 t4=t2;
954 t6=t5;
955 for(k=0;k<l1;k++){
956 ch[t3]=cc[t6-1]+cc[t6-1];
957 ch[t4]=cc[t6]+cc[t6];
958 t3+=ido;
959 t4+=ido;
960 t6+=t10;
961 }
962 t5+=t7;
963 }
964 if (ido == 1){
965 state=116;
966 break;
967 }
968 if(nbd<l1){
969 state=112;
970 break;
971 }
972
973 t1=0;
974 t2=ipp2*t0;
975 t7=0;
976 for(j=1;j<ipph;j++){
977 t1+=t0;
978 t2-=t0;
979 t3=t1;
980 t4=t2;
981
982 t7+=(ido<<1);
983 t8=t7;
984 for(k=0;k<l1;k++){
985 t5=t3;
986 t6=t4;
987 t9=t8;
988 t11=t8;
989 for(i=2;i<ido;i+=2){
990 t5+=2;
991 t6+=2;
992 t9+=2;
993 t11-=2;
994 ch[t5-1]=cc[t9-1]+cc[t11-1];
995 ch[t6-1]=cc[t9-1]-cc[t11-1];
996 ch[t5]=cc[t9]-cc[t11];
997 ch[t6]=cc[t9]+cc[t11];
998 }
999 t3+=ido;
1000 t4+=ido;
1001 t8+=t10;
1002 }
1003 }
1004 state=116;
1005 break;
1006 case 112:
1007 t1=0;
1008 t2=ipp2*t0;
1009 t7=0;
1010 for(j=1;j<ipph;j++){
1011 t1+=t0;
1012 t2-=t0;
1013 t3=t1;
1014 t4=t2;
1015 t7+=(ido<<1);
1016 t8=t7;
1017 t9=t7;
1018 for(i=2;i<ido;i+=2){
1019 t3+=2;
1020 t4+=2;
1021 t8+=2;
1022 t9-=2;
1023 t5=t3;
1024 t6=t4;
1025 t11=t8;
1026 t12=t9;
1027 for(k=0;k<l1;k++){
1028 ch[t5-1]=cc[t11-1]+cc[t12-1];
1029 ch[t6-1]=cc[t11-1]-cc[t12-1];
1030 ch[t5]=cc[t11]-cc[t12];
1031 ch[t6]=cc[t11]+cc[t12];
1032 t5+=ido;
1033 t6+=ido;
1034 t11+=t10;
1035 t12+=t10;
1036 }
1037 }
1038 }
1039 case 116:
1040 ar1=1.f;
1041 ai1=0.f;
1042 t1=0;
1043 t9=(t2=ipp2*idl1);
1044 t3=(ip-1)*idl1;
1045 for(l=1;l<ipph;l++){
1046 t1+=idl1;
1047 t2-=idl1;
1048
1049 ar1h=dcp*ar1-dsp*ai1;
1050 ai1=dcp*ai1+dsp*ar1;
1051 ar1=ar1h;
1052 t4=t1;
1053 t5=t2;
1054 t6=0;
1055 t7=idl1;
1056 t8=t3;
1057 for(ik=0;ik<idl1;ik++){
1058 c2[t4++]=ch2[t6++]+ar1*ch2[t7++];
1059 c2[t5++]=ai1*ch2[t8++];
1060 }
1061 dc2=ar1;
1062 ds2=ai1;
1063 ar2=ar1;
1064 ai2=ai1;
1065
1066 t6=idl1;
1067 t7=t9-idl1;
1068 for(j=2;j<ipph;j++){
1069 t6+=idl1;
1070 t7-=idl1;
1071 ar2h=dc2*ar2-ds2*ai2;
1072 ai2=dc2*ai2+ds2*ar2;
1073 ar2=ar2h;
1074 t4=t1;
1075 t5=t2;
1076 t11=t6;
1077 t12=t7;
1078 for(ik=0;ik<idl1;ik++){
1079 c2[t4++]+=ar2*ch2[t11++];
1080 c2[t5++]+=ai2*ch2[t12++];
1081 }
1082 }
1083 }
1084
1085 t1=0;
1086 for(j=1;j<ipph;j++){
1087 t1+=idl1;
1088 t2=t1;
1089 for(ik=0;ik<idl1;ik++)ch2[ik]+=ch2[t2++];
1090 }
1091
1092 t1=0;
1093 t2=ipp2*t0;
1094 for(j=1;j<ipph;j++){
1095 t1+=t0;
1096 t2-=t0;
1097 t3=t1;
1098 t4=t2;
1099 for(k=0;k<l1;k++){
1100 ch[t3]=c1[t3]-c1[t4];
1101 ch[t4]=c1[t3]+c1[t4];
1102 t3+=ido;
1103 t4+=ido;
1104 }
1105 }
1106
1107 if(ido==1){
1108 state=132;
1109 break;
1110 }
1111 if(nbd<l1){
1112 state=128;
1113 break;
1114 }
1115
1116 t1=0;
1117 t2=ipp2*t0;
1118 for(j=1;j<ipph;j++){
1119 t1+=t0;
1120 t2-=t0;
1121 t3=t1;
1122 t4=t2;
1123 for(k=0;k<l1;k++){
1124 t5=t3;
1125 t6=t4;
1126 for(i=2;i<ido;i+=2){
1127 t5+=2;
1128 t6+=2;
1129 ch[t5-1]=c1[t5-1]-c1[t6];
1130 ch[t6-1]=c1[t5-1]+c1[t6];
1131 ch[t5]=c1[t5]+c1[t6-1];
1132 ch[t6]=c1[t5]-c1[t6-1];
1133 }
1134 t3+=ido;
1135 t4+=ido;
1136 }
1137 }
1138 state=132;
1139 break;
1140 case 128:
1141 t1=0;
1142 t2=ipp2*t0;
1143 for(j=1;j<ipph;j++){
1144 t1+=t0;
1145 t2-=t0;
1146 t3=t1;
1147 t4=t2;
1148 for(i=2;i<ido;i+=2){
1149 t3+=2;
1150 t4+=2;
1151 t5=t3;
1152 t6=t4;
1153 for(k=0;k<l1;k++){
1154 ch[t5-1]=c1[t5-1]-c1[t6];
1155 ch[t6-1]=c1[t5-1]+c1[t6];
1156 ch[t5]=c1[t5]+c1[t6-1];
1157 ch[t6]=c1[t5]-c1[t6-1];
1158 t5+=ido;
1159 t6+=ido;
1160 }
1161 }
1162 }
1163 case 132:
1164 if(ido==1)return;
1165
1166 for(ik=0;ik<idl1;ik++)c2[ik]=ch2[ik];
1167
1168 t1=0;
1169 for(j=1;j<ip;j++){
1170 t2=(t1+=t0);
1171 for(k=0;k<l1;k++){
1172 c1[t2]=ch[t2];
1173 t2+=ido;
1174 }
1175 }
1176
1177 if(nbd>l1){
1178 state=139;
1179 break;
1180 }
1181
1182 is= -ido-1;
1183 t1=0;
1184 for(j=1;j<ip;j++){
1185 is+=ido;
1186 t1+=t0;
1187 idij=is;
1188 t2=t1;
1189 for(i=2;i<ido;i+=2){
1190 t2+=2;
1191 idij+=2;
1192 t3=t2;
1193 for(k=0;k<l1;k++){
1194 c1[t3-1]=wa[index+idij-1]*ch[t3-1]-wa[index+idij]*ch[t3];
1195 c1[t3]=wa[index+idij-1]*ch[t3]+wa[index+idij]*ch[t3-1];
1196 t3+=ido;
1197 }
1198 }
1199 }
1200 return;
1201
1202 case 139:
1203 is= -ido-1;
1204 t1=0;
1205 for(j=1;j<ip;j++){
1206 is+=ido;
1207 t1+=t0;
1208 t2=t1;
1209 for(k=0;k<l1;k++){
1210 idij=is;
1211 t3=t2;
1212 for(i=2;i<ido;i+=2){
1213 idij+=2;
1214 t3+=2;
1215 c1[t3-1]=wa[index+idij-1]*ch[t3-1]-wa[index+idij]*ch[t3];
1216 c1[t3]=wa[index+idij-1]*ch[t3]+wa[index+idij]*ch[t3-1];
1217 }
1218 t2+=ido;
1219 }
1220 }
1221 break loop;
1222 }
1223 }
1224 }
1225
1226 static void drftb1(int n, float[] c, float[] ch, float[] wa, int index, int[] ifac){
1227 int i,k1,l1,l2=0;
1228 int na;
1229 int nf,ip=0,iw,ix2,ix3,ido=0,idl1=0;
1230
1231 nf=ifac[1];
1232 na=0;
1233 l1=1;
1234 iw=1;
1235
1236 for(k1=0;k1<nf;k1++){
1237 int state=100;
1238 loop: while(true){
1239 switch(state){
1240 case 100:
1241 ip=ifac[k1 + 2];
1242 l2=ip*l1;
1243 ido=n/l2;
1244 idl1=ido*l1;
1245 if(ip!=4){
1246 state=103;
1247 break;
1248 }
1249 ix2=iw+ido;
1250 ix3=ix2+ido;
1251
1252 if(na!=0)
1253 dradb4(ido,l1,ch,c,wa,index+iw-1,wa,index+ix2-1,wa,index+ix3-1);
1254 else
1255 dradb4(ido,l1,c,ch,wa,index+iw-1,wa,index+ix2-1,wa,index+ix3-1);
1256 na=1-na;
1257 state=115;
1258 break;
1259 case 103:
1260 if(ip!=2){
1261 state=106;
1262 break;
1263 }
1264
1265 if(na!=0)
1266 dradb2(ido,l1,ch,c,wa,index+iw-1);
1267 else
1268 dradb2(ido,l1,c,ch,wa,index+iw-1);
1269 na=1-na;
1270 state=115;
1271 break;
1272
1273 case 106:
1274 if(ip!=3){
1275 state=109;
1276 break;
1277 }
1278
1279 ix2=iw+ido;
1280 if(na!=0)
1281 dradb3(ido,l1,ch,c,wa,index+iw-1,wa,index+ix2-1);
1282 else
1283 dradb3(ido,l1,c,ch,wa,index+iw-1,wa,index+ix2-1);
1284 na=1-na;
1285 state=115;
1286 break;
1287 case 109:
1288 // The radix five case can be translated later.....
1289 // if(ip!=5)goto L112;
1290 //
1291 //ix2=iw+ido;
1292 //ix3=ix2+ido;
1293 //ix4=ix3+ido;
1294 //if(na!=0)
1295 // dradb5(ido,l1,ch,c,wa+iw-1,wa+ix2-1,wa+ix3-1,wa+ix4-1);
1296 //else
1297 // dradb5(ido,l1,c,ch,wa+iw-1,wa+ix2-1,wa+ix3-1,wa+ix4-1);
1298 //na=1-na;
1299 //state=115;
1300 //break;
1301 if(na!=0)
1302 dradbg(ido,ip,l1,idl1,ch,ch,ch,c,c,wa,index+iw-1);
1303 else
1304 dradbg(ido,ip,l1,idl1,c,c,c,ch,ch,wa,index+iw-1);
1305 if(ido==1)na=1-na;
1306
1307 case 115:
1308 l1=l2;
1309 iw+=(ip-1)*ido;
1310 break loop;
1311 }
1312 }
1313 }
1314 if(na==0)return;
1315 for(i=0;i<n;i++)c[i]=ch[i];
1316 }
1317}
This page took 0.229743 seconds and 4 git commands to generate.