This graph shows which files directly or indirectly include this file:
Go to the source code of this file.
Defines | |
#define | tab1(i) tab1 [(i)-1] |
#define | xcmplx(i, j) xcmplx [((j)-1)*2 + (i)-1] |
#define | br(i) br [(i)-1] |
#define | bi(i) bi [(i)-1] |
Functions | |
void | fftr_q (float *xcmplx, int nv) |
void | fftr_d (double *xcmplx, int nv) |
void | fftc_q (float *br, float *bi, int ln, int ks) |
void | fftc_d (double *br, double *bi, int ln, int ks) |
|
|
|
|
|
|
|
|
|
Definition at line 366 of file spidfft.cpp. References abs, bi, br, status, t, and tab1. Referenced by fftr_d(). 00367 { 00368 double rni,sgn,tr1,tr2,ti1,ti2; 00369 double cc,c,ss,s,t,x2,x3,x4,x5; 00370 int b3,b4,b5,b6,b7,b56; 00371 int n, k, l, j, i, ix0, ix1, status=0; 00372 double tab1[15]; 00373 00374 tab1(1)=9.58737990959775e-5; 00375 tab1(2)=1.91747597310703e-4; 00376 tab1(3)=3.83495187571395e-4; 00377 tab1(4)=7.66990318742704e-4; 00378 tab1(5)=1.53398018628476e-3; 00379 tab1(6)=3.06795676296598e-3; 00380 tab1(7)=6.13588464915449e-3; 00381 tab1(8)=1.22715382857199e-2; 00382 tab1(9)=2.45412285229123e-2; 00383 tab1(10)=4.90676743274181e-2; 00384 tab1(11)=9.80171403295604e-2; 00385 tab1(12)=1.95090322016128e-1; 00386 tab1(13)=3.82683432365090e-1; 00387 tab1(14)=7.07106781186546e-1; 00388 tab1(15)=1.00000000000000; 00389 00390 n=(int)pow(2,ln); 00391 00392 k=abs(ks); 00393 l=16-ln; 00394 b3=n*k; 00395 b6=b3; 00396 b7=k; 00397 if (ks > 0) { 00398 sgn=1.0; 00399 } 00400 else { 00401 sgn=-1.0; 00402 rni=1.0/(float)(n); 00403 j=1; 00404 for (i=1;i<=n;i++) { 00405 br(j)=br(j)*rni; 00406 bi(j)=bi(j)*rni; 00407 j=j+k; 00408 } 00409 } 00410 00411 L12: 00412 b6=b6/2; 00413 b5=b6; 00414 b4=2*b6; 00415 b56=b5-b6; 00416 00417 L14: 00418 tr1=br(b5+1); 00419 ti1=bi(b5+1); 00420 tr2=br(b56+1); 00421 ti2=bi(b56+1); 00422 00423 br(b5+1)=tr2-tr1; 00424 bi(b5+1)=ti2-ti1; 00425 br(b56+1)=tr1+tr2; 00426 bi(b56+1)=ti1+ti2; 00427 00428 b5=b5+b4; 00429 b56=b5-b6; 00430 if ( b5 <= b3 ) goto L14; 00431 if ( b6 == b7 ) goto L20; 00432 00433 b4=b7; 00434 cc=2.0*pow(tab1(l),2); 00435 c=1.0-cc; 00436 l++; 00437 ss=sgn*tab1(l); 00438 s=ss; 00439 00440 L16: 00441 b5=b6+b4; 00442 b4=2*b6; 00443 b56=b5-b6; 00444 00445 L18: 00446 tr1=br(b5+1); 00447 ti1=bi(b5+1); 00448 tr2=br(b56+1); 00449 ti2=bi(b56+1); 00450 br(b5+1)=c*(tr2-tr1)-s*(ti2-ti1); 00451 bi(b5+1)=s*(tr2-tr1)+c*(ti2-ti1); 00452 br(b56+1)=tr1+tr2; 00453 bi(b56+1)=ti1+ti2; 00454 00455 b5=b5+b4; 00456 b56=b5-b6; 00457 if ( b5 <= b3 ) goto L18; 00458 b4=b5-b6; 00459 b5=b4-b3; 00460 c=-c; 00461 b4=b6-b5; 00462 if ( b5 < b4 ) goto L16; 00463 b4=b4+b7; 00464 if ( b4 >= b5 ) goto L12; 00465 00466 t=c-cc*c-ss*s; 00467 s=s+ss*c-cc*s; 00468 c=t; 00469 goto L16; 00470 00471 L20: 00472 ix0=b3/2; 00473 b3=b3-b7; 00474 b4=0; 00475 b5=0; 00476 b6=ix0; 00477 ix1=0; 00478 if (b6 == b7) goto EXIT; 00479 00480 L22: 00481 b4=b3-b4; 00482 b5=b3-b5; 00483 x2=br(b4+1); 00484 x3=br(b5+1); 00485 x4=bi(b4+1); 00486 x5=bi(b5+1); 00487 br(b4+1)=x3; 00488 br(b5+1)=x2; 00489 bi(b4+1)=x5; 00490 bi(b5+1)=x4; 00491 if(b6 < b4) goto L22; 00492 00493 L24: 00494 b4=b4+b7; 00495 b5=b6+b5; 00496 x2=br(b4+1); 00497 x3=br(b5+1); 00498 x4=bi(b4+1); 00499 x5=bi(b5+1); 00500 br(b4+1)=x3; 00501 br(b5+1)=x2; 00502 bi(b4+1)=x5; 00503 bi(b5+1)=x4; 00504 ix0=b6; 00505 00506 L26: 00507 ix0=ix0/2; 00508 ix1=ix1-ix0; 00509 if( ix1 >= 0) goto L26; 00510 00511 ix0=2*ix0; 00512 b4=b4+b7; 00513 ix1=ix1+ix0; 00514 b5=ix1; 00515 if ( b5 >= b4) goto L22; 00516 if ( b4 < b6) goto L24; 00517 00518 EXIT: 00519 status = 0; 00520 }
|
|
Definition at line 135 of file spidfft.cpp. 00136 { 00137 // dimension br(1),bi(1) 00138 00139 int b3,b4,b5,b6,b7,b56; 00140 int n, k, l, j, i, ix0, ix1; 00141 float rni, tr1, ti1, tr2, ti2, cc, c, ss, s, t, x2, x3, x4, x5, sgn; 00142 float tab1[15]; 00143 int status=0; 00144 00145 tab1(1)=9.58737990959775e-5; 00146 tab1(2)=1.91747597310703e-4; 00147 tab1(3)=3.83495187571395e-4; 00148 tab1(4)=7.66990318742704e-4; 00149 tab1(5)=1.53398018628476e-3; 00150 tab1(6)=3.06795676296598e-3; 00151 tab1(7)=6.13588464915449e-3; 00152 tab1(8)=1.22715382857199e-2; 00153 tab1(9)=2.45412285229123e-2; 00154 tab1(10)=4.90676743274181e-2; 00155 tab1(11)=9.80171403295604e-2; 00156 tab1(12)=1.95090322016128e-1; 00157 tab1(13)=3.82683432365090e-1; 00158 tab1(14)=7.07106781186546e-1; 00159 tab1(15)=1.00000000000000; 00160 00161 n=(int)pow(2,ln); 00162 00163 k=abs(ks); 00164 l=16-ln; 00165 b3=n*k; 00166 b6=b3; 00167 b7=k; 00168 if( ks > 0 ) { 00169 sgn=1.0; 00170 } 00171 else { 00172 sgn=-1.0; 00173 rni=1.0/(float)n; 00174 j=1; 00175 for (i=1; i<=n;i++) { 00176 br(j)=br(j)*rni; 00177 bi(j)=bi(j)*rni; 00178 j=j+k; 00179 } 00180 } 00181 L12: 00182 b6=b6/2; 00183 b5=b6; 00184 b4=2*b6; 00185 b56=b5-b6; 00186 L14: 00187 tr1=br(b5+1); 00188 ti1=bi(b5+1); 00189 tr2=br(b56+1); 00190 ti2=bi(b56+1); 00191 00192 br(b5+1)=tr2-tr1; 00193 bi(b5+1)=ti2-ti1; 00194 br(b56+1)=tr1+tr2; 00195 bi(b56+1)=ti1+ti2; 00196 00197 b5=b5+b4; 00198 b56=b5-b6; 00199 if (b5 <= b3) goto L14; 00200 if (b6 == b7) goto L20; 00201 00202 b4=b7; 00203 cc=2.0*pow(tab1(l),2); 00204 c=1.0-cc; 00205 l=l+1; 00206 ss=sgn*tab1(l); 00207 s=ss; 00208 L16: 00209 b5=b6+b4; 00210 b4=2*b6; 00211 b56=b5-b6; 00212 L18: 00213 tr1=br(b5+1); 00214 ti1=bi(b5+1); 00215 tr2=br(b56+1); 00216 ti2=bi(b56+1); 00217 br(b5+1)=c*(tr2-tr1)-s*(ti2-ti1); 00218 bi(b5+1)=s*(tr2-tr1)+c*(ti2-ti1); 00219 br(b56+1)=tr1+tr2; 00220 bi(b56+1)=ti1+ti2; 00221 00222 b5=b5+b4; 00223 b56=b5-b6; 00224 if(b5 <= b3) goto L18; 00225 b4=b5-b6; 00226 b5=b4-b3; 00227 c=-c; 00228 b4=b6-b5; 00229 if(b5 < b4) goto L16; 00230 b4=b4+b7; 00231 if(b4 >= b5) goto L12; 00232 00233 t=c-cc*c-ss*s; 00234 s=s+ss*c-cc*s; 00235 c=t; 00236 goto L16; 00237 L20: 00238 ix0=b3/2; 00239 b3=b3-b7; 00240 b4=0; 00241 b5=0; 00242 b6=ix0; 00243 ix1=0; 00244 if ( b6 == b7) goto EXIT; 00245 L22: 00246 b4=b3-b4; 00247 b5=b3-b5; 00248 x2=br(b4+1); 00249 x3=br(b5+1); 00250 x4=bi(b4+1); 00251 x5=bi(b5+1); 00252 br(b4+1)=x3; 00253 br(b5+1)=x2; 00254 bi(b4+1)=x5; 00255 bi(b5+1)=x4; 00256 if (b6 < b4) goto L22; 00257 L24: 00258 b4=b4+b7; 00259 b5=b6+b5; 00260 x2=br(b4+1); 00261 x3=br(b5+1); 00262 x4=bi(b4+1); 00263 x5=bi(b5+1); 00264 br(b4+1)=x3; 00265 br(b5+1)=x2; 00266 bi(b4+1)=x5; 00267 bi(b5+1)=x4; 00268 ix0=b6; 00269 L26: 00270 ix0=ix0/2; 00271 ix1=ix1-ix0; 00272 if(ix1 >= 0) goto L26; 00273 00274 ix0=2*ix0; 00275 b4=b4+b7; 00276 ix1=ix1+ix0; 00277 b5=ix1; 00278 if (b5 >= b4) goto L22; 00279 if (b4 < b6) goto L24; 00280 EXIT: 00281 status = 0; 00282 }
|
|
Definition at line 285 of file spidfft.cpp. References abs, fftc_d(), t, tab1, and xcmplx. Referenced by crosrng_e(), and crosrng_ms(). 00286 { 00287 // double precision x(2,1) 00288 int i1, i2, nu, inv, nu1, n, isub, n2, i; 00289 double tr1,tr2,ti1,ti2,tr,ti; 00290 double cc,c,ss,s,t; 00291 double tab1[15]; 00292 00293 tab1(1)=9.58737990959775e-5; 00294 tab1(2)=1.91747597310703e-4; 00295 tab1(3)=3.83495187571395e-4; 00296 tab1(4)=7.66990318742704e-4; 00297 tab1(5)=1.53398018628476e-3; 00298 tab1(6)=3.06795676296598e-3; 00299 tab1(7)=6.13588464915449e-3; 00300 tab1(8)=1.22715382857199e-2; 00301 tab1(9)=2.45412285229123e-2; 00302 tab1(10)=4.90676743274181e-2; 00303 tab1(11)=9.80171403295604e-2; 00304 tab1(12)=1.95090322016128e-1; 00305 tab1(13)=3.82683432365090e-1; 00306 tab1(14)=7.07106781186546e-1; 00307 tab1(15)=1.00000000000000; 00308 00309 nu=abs(nv); 00310 inv=nv/nu; 00311 nu1=nu-1; 00312 n=(int)pow(2,nu1); 00313 isub=16-nu1; 00314 ss=-tab1(isub); 00315 cc=-2.0*pow(tab1(isub-1),2); 00316 c=1.0; 00317 s=0.0; 00318 n2=n/2; 00319 00320 if ( inv > 0 ) { 00321 fftc_d(&xcmplx(1,1),&xcmplx(2,1),nu1,2); 00322 tr=xcmplx(1,1); 00323 ti=xcmplx(2,1); 00324 xcmplx(1,1)=tr+ti; 00325 xcmplx(2,1)=tr-ti; 00326 for (i=1;i<=n2;i++) { 00327 i1=i+1; 00328 i2=n-i+1; 00329 tr1=xcmplx(1,i1); 00330 tr2=xcmplx(1,i2); 00331 ti1=xcmplx(2,i1); 00332 ti2=xcmplx(2,i2); 00333 t=(cc*c-ss*s)+c; 00334 s=(cc*s+ss*c)+s; 00335 c=t; 00336 xcmplx(1,i1)=0.5*((tr1+tr2)+(ti1+ti2)*c-(tr1-tr2)*s); 00337 xcmplx(1,i2)=0.5*((tr1+tr2)-(ti1+ti2)*c+(tr1-tr2)*s); 00338 xcmplx(2,i1)=0.5*((ti1-ti2)-(ti1+ti2)*s-(tr1-tr2)*c); 00339 xcmplx(2,i2)=0.5*(-(ti1-ti2)-(ti1+ti2)*s-(tr1-tr2)*c); 00340 } 00341 } 00342 else { 00343 tr=xcmplx(1,1); 00344 ti=xcmplx(2,1); 00345 xcmplx(1,1)=0.5*(tr+ti); 00346 xcmplx(2,1)=0.5*(tr-ti); 00347 for (i=1;i<=n2;i++) { 00348 i1=i+1; 00349 i2=n-i+1; 00350 tr1=xcmplx(1,i1); 00351 tr2=xcmplx(1,i2); 00352 ti1=xcmplx(2,i1); 00353 ti2=xcmplx(2,i2); 00354 t=(cc*c-ss*s)+c; 00355 s=(cc*s+ss*c)+s; 00356 c=t; 00357 xcmplx(1,i1)=0.5*((tr1+tr2)-(tr1-tr2)*s-(ti1+ti2)*c); 00358 xcmplx(1,i2)=0.5*((tr1+tr2)+(tr1-tr2)*s+(ti1+ti2)*c); 00359 xcmplx(2,i1)=0.5*((ti1-ti2)+(tr1-tr2)*c-(ti1+ti2)*s); 00360 xcmplx(2,i2)=0.5*(-(ti1-ti2)+(tr1-tr2)*c-(ti1+ti2)*s); 00361 } 00362 fftc_d(&xcmplx(1,1),&xcmplx(2,1),nu1,-2); 00363 } 00364 }
|
|
Definition at line 53 of file spidfft.cpp. 00054 { 00055 // dimension xcmplx(2,1); xcmplx(1,i) --- real, xcmplx(2,i) --- imaginary 00056 00057 float tab1[15]; 00058 int nu, inv, nu1, n, isub, n2, i1, i2, i; 00059 float ss, cc, c, s, tr, ti, tr1, tr2, ti1, ti2, t; 00060 00061 tab1(1)=9.58737990959775e-5; 00062 tab1(2)=1.91747597310703e-4; 00063 tab1(3)=3.83495187571395e-4; 00064 tab1(4)=7.66990318742704e-4; 00065 tab1(5)=1.53398018628476e-3; 00066 tab1(6)=3.06795676296598e-3; 00067 tab1(7)=6.13588464915449e-3; 00068 tab1(8)=1.22715382857199e-2; 00069 tab1(9)=2.45412285229123e-2; 00070 tab1(10)=4.90676743274181e-2; 00071 tab1(11)=9.80171403295604e-2; 00072 tab1(12)=1.95090322016128e-1; 00073 tab1(13)=3.82683432365090e-1; 00074 tab1(14)=7.07106781186546e-1; 00075 tab1(15)=1.00000000000000; 00076 00077 nu=abs(nv); 00078 inv=nv/nu; 00079 nu1=nu-1; 00080 n=(int)pow(2,nu1); 00081 isub=16-nu1; 00082 00083 ss=-tab1(isub); 00084 cc=-2.0*pow(tab1(isub-1),2); 00085 c=1.0; 00086 s=0.0; 00087 n2=n/2; 00088 if ( inv > 0) { 00089 fftc_q(&xcmplx(1,1),&xcmplx(2,1),nu1,2); 00090 tr=xcmplx(1,1); 00091 ti=xcmplx(2,1); 00092 xcmplx(1,1)=tr+ti; 00093 xcmplx(2,1)=tr-ti; 00094 for (i=1;i<=n2;i++) { 00095 i1=i+1; 00096 i2=n-i+1; 00097 tr1=xcmplx(1,i1); 00098 tr2=xcmplx(1,i2); 00099 ti1=xcmplx(2,i1); 00100 ti2=xcmplx(2,i2); 00101 t=(cc*c-ss*s)+c; 00102 s=(cc*s+ss*c)+s; 00103 c=t; 00104 xcmplx(1,i1)=0.5*((tr1+tr2)+(ti1+ti2)*c-(tr1-tr2)*s); 00105 xcmplx(1,i2)=0.5*((tr1+tr2)-(ti1+ti2)*c+(tr1-tr2)*s); 00106 xcmplx(2,i1)=0.5*((ti1-ti2)-(ti1+ti2)*s-(tr1-tr2)*c); 00107 xcmplx(2,i2)=0.5*(-(ti1-ti2)-(ti1+ti2)*s-(tr1-tr2)*c); 00108 } 00109 } 00110 else { 00111 tr=xcmplx(1,1); 00112 ti=xcmplx(2,1); 00113 xcmplx(1,1)=0.5*(tr+ti); 00114 xcmplx(2,1)=0.5*(tr-ti); 00115 for (i=1; i<=n2; i++) { 00116 i1=i+1; 00117 i2=n-i+1; 00118 tr1=xcmplx(1,i1); 00119 tr2=xcmplx(1,i2); 00120 ti1=xcmplx(2,i1); 00121 ti2=xcmplx(2,i2); 00122 t=(cc*c-ss*s)+c; 00123 s=(cc*s+ss*c)+s; 00124 c=t; 00125 xcmplx(1,i1)=0.5*((tr1+tr2)-(tr1-tr2)*s-(ti1+ti2)*c); 00126 xcmplx(1,i2)=0.5*((tr1+tr2)+(tr1-tr2)*s+(ti1+ti2)*c); 00127 xcmplx(2,i1)=0.5*((ti1-ti2)+(tr1-tr2)*c-(ti1+ti2)*s); 00128 xcmplx(2,i2)=0.5*(-(ti1-ti2)+(tr1-tr2)*c-(ti1+ti2)*s); 00129 } 00130 fftc_q(&xcmplx(1,1),&xcmplx(2,1),nu1,-2); 00131 } 00132 }
|