media/libvorbis/lib/vorbis_smallft.c

Thu, 22 Jan 2015 13:21:57 +0100

author
Michael Schloh von Bennewitz <michael@schloh.com>
date
Thu, 22 Jan 2015 13:21:57 +0100
branch
TOR_BUG_9701
changeset 15
b8a032363ba2
permissions
-rw-r--r--

Incorporate requested changes from Mozilla in review:
https://bugzilla.mozilla.org/show_bug.cgi?id=1123480#c6

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

mercurial