dbf-halloween2015

annotate libs/vorbis/smallft.c @ 3:c37fe5d8a4ed

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