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