vrshoot

annotate libs/vorbis/smallft.c @ 0:b2f14e535253

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