EMAN2
spidfft.cpp
Go to the documentation of this file.
00001 /*
00002  * Author: Steven Ludtke, 04/10/2003 (sludtke@bcm.edu)
00003  * Copyright (c) 2000-2006 Baylor College of Medicine
00004  *
00005  * This software is issued under a joint BSD/GNU license. You may use the
00006  * source code in this file under either license. However, note that the
00007  * complete EMAN2 and SPARX software packages have some GPL dependencies,
00008  * so you are responsible for compliance with the licenses of these packages
00009  * if you opt to use BSD licensing. The warranty disclaimer below holds
00010  * in either instance.
00011  *
00012  * This complete copyright notice must be included in any revised version of the
00013  * source code. Additional authorship citations may be added, but existing
00014  * author citations must be preserved.
00015  *
00016  * This program is free software; you can redistribute it and/or modify
00017  * it under the terms of the GNU General Public License as published by
00018  * the Free Software Foundation; either version 2 of the License, or
00019  * (at your option) any later version.
00020  *
00021  * This program is distributed in the hope that it will be useful,
00022  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00023  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00024  * GNU General Public License for more details.
00025  *
00026  * You should have received a copy of the GNU General Public License
00027  * along with this program; if not, write to the Free Software
00028  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
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    // 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 }
00133 
00134 // -----------------------------------------------------------------
00135 void fftc_q(float *br, float *bi, int ln, int ks)
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 }
00283 
00284 // -------------------------------------------
00285 void  fftr_d(double *xcmplx, int nv) 
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 } 
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