Thu, 22 Jan 2015 13:21:57 +0100
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 | } |