00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032 #include "ctf.h"
00033 #include "log.h"
00034 #include "emdata.h"
00035 #include "xydata.h"
00036 #include "Assert.h"
00037 #include "projector.h"
00038
00039 #include <sys/types.h>
00040 #include <sys/times.h>
00041 #include <time.h>
00042 #include <sys/time.h>
00043
00044 using namespace EMAN;
00045
00046 using std::cout;
00047 using std::cin;
00048 using std::endl;
00049 using std::string;
00050
00051 #include "spidfft.h"
00052
00053 void fftr_q(float *xcmplx, int nv)
00054 {
00055
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 }
00133
00134
00135 void fftc_q(float *br, float *bi, int ln, int ks)
00136 {
00137
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 }
00283
00284
00285 void fftr_d(double *xcmplx, int nv)
00286 {
00287
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 }
00365
00366 void fftc_d(double *br, double *bi, int ln, int ks)
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 }
00521