fft.cpp

普林斯顿开发的快速球面调和变换算法

// FFT.cpp: implementation of the FFT class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "FFT.h"
#include "math.h"

#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

FFT::FFT()
{

}

FFT::~FFT()
{

}
int FFT::cosqi_(int *n, double *wsave)
{
		int i__1;

    /* Builtin functions */
//    double cos(doublereal);

    /* Local variables */
    static int k;
    static double fk, dt, pih, dum;
//    extern /* Subroutine */ int rffti_(int *, double *);
//    extern doublereal pimach_(double *);


    /* Parameter adjustments */
    --wsave;

    /* Function Body */
    pih = pimach_(&dum) * .5f;
    dt = pih / (double) (*n);
    fk = 0.f;
    i__1 = *n;
    for (k = 1; k <= i__1; ++k) {
	fk += 1.f;
	wsave[k] = cos(fk * dt);
/* L101: */
    }
    rffti_(n, &wsave[*n + 1]);
    return 0;

}

double FFT::pimach_(double *dum)
{
		double ret_val;

    /* Builtin functions */
    //double atan;
/*     PI=3.1415926535897932384626433832795028841971693993751058209749446 */
    ret_val = atan(1.f) * 4.f;
    return ret_val;

}

int FFT::rffti_(int *n, double *wsave)
{
	    /* Parameter adjustments */
    --wsave;

    /* Function Body */
    if (*n == 1) {
	return 0;
    }
    this->rffti1_(n, &wsave[*n + 1], &wsave[(*n << 1) + 1]);
    return 0;


}

int FFT::rffti1_(int *n, double *wa, double *ifac)
{
		 static int ntryh[4] = { 4,2,3,5 };

    /* System generated locals */
    int i__1, i__2, i__3;

    /* Builtin functions */
//    double cos(doublereal), sin(doublereal);

    /* Local variables */
    static int i__, j, k1, l1, l2, ib;
    static double fi;
    static int ld, ii, nf, ip, nl, is, nq, nr;
    static double arg;
    static int ido;
    static double dum;
    static int ipm;
    static double tpi;
    static int nfm1;
    static double argh;
    static int ntry;
    static double argld;

    /* Parameter adjustments */
    --ifac;
    --wa;

    /* Function Body */
    nl = *n;
    nf = 0;
    j = 0;
L101:
    ++j;
    if (j - 4 <= 0) {
	goto L102;
    } else {
	goto L103;
    }
L102:
    ntry = ntryh[j - 1];
    goto L104;
L103:
    ntry += 2;
L104:
    nq = nl / ntry;
    nr = nl - ntry * nq;
    if (nr != 0) {
	goto L101;
    } else {
	goto L105;
    }
L105:
    ++nf;
    ifac[nf + 2] = ntry;
    nl = nq;
    if (ntry != 2) {
	goto L107;
    }
    if (nf == 1) {
	goto L107;
    }
    i__1 = nf;
    for (i__ = 2; i__ <= i__1; ++i__) {
	ib = nf - i__ + 2;
	ifac[ib + 2] = ifac[ib + 1];
/* L106: */
    }
    ifac[3] = 2;
L107:
    if (nl != 1) {
	goto L104;
    }
    ifac[1] = *n;
    ifac[2] = nf;
    tpi = pimach_(&dum) * 2.f;
    argh = tpi / (double) (*n);
    is = 0;
    nfm1 = nf - 1;
    l1 = 1;
    if (nfm1 == 0) {
	return 0;
    }
    i__1 = nfm1;
    for (k1 = 1; k1 <= i__1; ++k1) {
	ip = ifac[k1 + 2];
	ld = 0;
	l2 = l1 * ip;
	ido = *n / l2;
	ipm = ip - 1;
	i__2 = ipm;
	for (j = 1; j <= i__2; ++j) {
	    ld += l1;
	    i__ = is;
	    argld = (double) ld * argh;
	    fi = 0.f;
	    i__3 = ido;
	    for (ii = 3; ii <= i__3; ii += 2) {
		i__ += 2;
		fi += 1.f;
		arg = fi * argld;
		wa[i__ - 1] = cos(arg);
		wa[i__] = sin(arg);
/* L108: */
	    }
	    is += ido;
/* L109: */
	}
	l1 = l2;
/* L110: */
    }
    return 0;



}

int FFT::cosqb_(int *n, int *kk, double *x, double *wsave)
{
	    static double tsqrt2 = 2.82842712474619f;

    /* System generated locals */
    int i__1;

    /* Local variables */
    static double x1;
    
    /* Parameter adjustments */
    --wsave;
    --x;

    /* Function Body */

    if ((i__1 = *n - 2) < 0) {
	goto L101;
    } else if (i__1 == 0) {
	goto L102;
    } else {
	goto L103;
    }
L101:
    x[1] *= 4.f;
    return 0;
L102:
    x1 = (x[1] + x[2]) * 4.f;
    x[2] = tsqrt2 * (x[1] - x[2]);
    x[1] = x1;
    return 0;
L103:
    cosqb1_(n, kk, &x[1], &wsave[1], &wsave[*n + 1]);
    return 0;


}

int FFT::cosqb1_(int *n, int *kk, double *x, double *w, double *xh)
{
	    int i__1;

    /* Local variables */
    static int i__, k, kc, np2, ns2;
    static double xim1;
    static int modn;
    static double scale;
//    extern /* Subroutine */ int rfftb_(integer *, real *, real *);

    /* Parameter adjustments */
    --xh;
    --w;
    --x;

    /* Function Body */
    scale = 1.f / (*n * 2.f);
    ns2 = (*n + 1) / 2;
    np2 = *n + 2;
    i__1 = *n;
    for (i__ = 3; i__ <= i__1; i__ += 2) {
	xim1 = x[i__ - 1] + x[i__];
	x[i__] -= x[i__ - 1];
	x[i__ - 1] = xim1;
/* L101: */
    }
    x[1] += x[1];
    modn = *n % 2;
    if (modn == 0) {
	x[*n] += x[*n];
    }
    rfftb_(n, &x[1], &xh[1]);
/*     if asking for all the coefficients */
    if (*n == *kk) {
	i__1 = ns2;
	for (k = 2; k <= i__1; ++k) {
	    kc = np2 - k;
	    xh[k] = w[k - 1] * x[kc] + w[kc - 1] * x[k];
	    xh[kc] = w[k - 1] * x[k] - w[kc - 1] * x[kc];
/* L102: */
	}
	if (modn == 0) {
	    x[ns2 + 1] = scale * w[ns2] * (x[ns2 + 1] + x[ns2 + 1]);
	}
	i__1 = ns2;
	for (k = 2; k <= i__1; ++k) {
	    kc = np2 - k;
	    x[k] = scale * (xh[k] + xh[kc]);
	    x[kc] = scale * (xh[k] - xh[kc]);
/* L103: */
	}
	x[1] = scale * .5f * (x[1] + x[1]);
/*     otherwise */
    } else {
	i__1 = *kk;
	for (k = 2; k <= i__1; ++k) {
	    kc = np2 - k;
	    xh[k] = w[k - 1] * x[kc] + w[kc - 1] * x[k];
	    xh[kc] = w[k - 1] * x[k] - w[kc - 1] * x[kc];
/* L202: */
	}
	i__1 = *kk;
	for (k = 2; k <= i__1; ++k) {
	    kc = np2 - k;
	    x[k] = scale * (xh[k] + xh[kc]);
/* L303: */
	}
	x[1] = scale * .5f * (x[1] + x[1]);
    }
    return 0;
}


int FFT::rfftb_(int *n, double *r__, double *wsave)
{
//	extern /* Subroutine */ int rfftb1_(integer *, real *, real *, real *, real *);


    /* Parameter adjustments */
    --wsave;
    --r__;

    /* Function Body */
    if (*n == 1) {
	return 0;
    }
    rfftb1_(n, &r__[1], &wsave[1], &wsave[*n + 1], &wsave[(*n << 1) + 1]);
    return 0;
}



int FFT::rfftb1_(int *n, double *c__, double *ch, double *wa, double *ifac)
{
	    /* System generated locals */
    int i__1;

    /* Local variables */
    static int i__, k1, l1, l2, na, nf, ip, iw, ix2, ix3, ix4, ido, idl1;
/*    extern  int radb2_(integer *, integer *, real *, real *, 
	    real *), radb3_(integer *, integer *, real *, real *, real *, 
	    real *), radb4_(integer *, integer *, real *, real *, real *, 
	    real *, real *), radb5_(integer *, integer *, real *, real *, 
	    real *, real *, real *, real *), radbg_(integer *, integer *, 
	    integer *, integer *, real *, real *, real *, real *, real *, 
	    real *);
*/
    /* Parameter adjustments */
    --ifac;
    --wa;
    --ch;
    --c__;

    /* Function Body */
    nf = ifac[2];
    na = 0;
    l1 = 1;
    iw = 1;
    i__1 = nf;
    for (k1 = 1; k1 <= i__1; ++k1) {
	ip = ifac[k1 + 2];
	l2 = ip * l1;
	ido = *n / l2;
	idl1 = ido * l1;
	if (ip != 4) {
	    goto L103;
	}
	ix2 = iw + ido;
	ix3 = ix2 + ido;
	if (na != 0) {
	    goto L101;
	}
	radb4_(&ido, &l1, &c__[1], &ch[1], &wa[iw], &wa[ix2], &wa[ix3]);
	goto L102;
L101:
	radb4_(&ido, &l1, &ch[1], &c__[1], &wa[iw], &wa[ix2], &wa[ix3]);
L102:
	na = 1 - na;
	goto L115;
L103:
	if (ip != 2) {
	    goto L106;
	}
	if (na != 0) {
	    goto L104;
	}
	radb2_(&ido, &l1, &c__[1], &ch[1], &wa[iw]);
	goto L105;
L104:
	radb2_(&ido, &l1, &ch[1], &c__[1], &wa[iw]);
L105:
	na = 1 - na;
	goto L115;
L106:
	if (ip != 3) {
	    goto L109;
	}
	ix2 = iw + ido;
	if (na != 0) {
	    goto L107;
	}
	radb3_(&ido, &l1, &c__[1], &ch[1], &wa[iw], &wa[ix2]);
	goto L108;
L107:
	radb3_(&ido, &l1, &ch[1], &c__[1], &wa[iw], &wa[ix2]);
L108:
	na = 1 - na;
	goto L115;
L109:
	if (ip != 5) {
	    goto L112;
	}
	ix2 = iw + ido;
	ix3 = ix2 + ido;
	ix4 = ix3 + ido;
	if (na != 0) {
	    goto L110;
	}
	radb5_(&ido, &l1, &c__[1], &ch[1], &wa[iw], &wa[ix2], &wa[ix3], &wa[
		ix4]);
	goto L111;
L110:
	radb5_(&ido, &l1, &ch[1], &c__[1], &wa[iw], &wa[ix2], &wa[ix3], &wa[
		ix4]);
L111:
	na = 1 - na;
	goto L115;
L112:
	if (na != 0) {
	    goto L113;
	}
	radbg_(&ido, &ip, &l1, &idl1, &c__[1], &c__[1], &c__[1], &ch[1], &ch[
		1], &wa[iw]);
	goto L114;
L113:
	radbg_(&ido, &ip, &l1, &idl1, &ch[1], &ch[1], &ch[1], &c__[1], &c__[1]
		, &wa[iw]);
L114:
	if (ido == 1) {
	    na = 1 - na;
	}
L115:
	l1 = l2;
	iw += (ip - 1) * ido;
/* L116: */
    }
    if (na == 0) {
	return 0;
    }
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
	c__[i__] = ch[i__];
/* L117: */
    }
    return 0;


}

int FFT::radb2_(int *ido, int *l1, double *cc, double *ch, double *wa1)
{
	    /* System generated locals */
    int cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i__1, i__2;

    /* Local variables */
    static int i__, k, ic;
    static double ti2, tr2;
    static double idp2;

    /* Parameter adjustments */
    ch_dim1 = *ido;
    ch_dim2 = *l1;
    ch_offset = 1 + ch_dim1 * (1 + ch_dim2);
    ch -= ch_offset;
    cc_dim1 = *ido;
    cc_offset = 1 + cc_dim1 * 3;
    cc -= cc_offset;
    --wa1;

    /* Function Body */
    i__1 = *l1;
    for (k = 1; k <= i__1; ++k) {
	ch[(k + ch_dim2) * ch_dim1 + 1] = cc[((k << 1) + 1) * cc_dim1 + 1] + 
		cc[*ido + ((k << 1) + 2) * cc_dim1];
	ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1] = cc[((k << 1) + 1) * cc_dim1 
		+ 1] - cc[*ido + ((k << 1) + 2) * cc_dim1];
/* L101: */
    }
    if ((i__1 = *ido - 2) < 0) {
	goto L107;
    } else if (i__1 == 0) {
	goto L105;
    } else {
	goto L102;
    }
L102:
    idp2 = *ido + 2;
    i__1 = *l1;
    for (k = 1; k <= i__1; ++k) {
	i__2 = *ido;
	for (i__ = 3; i__ <= i__2; i__ += 2) {
	    ic = idp2 - i__;
	    ch[i__ - 1 + (k + ch_dim2) * ch_dim1] = cc[i__ - 1 + ((k << 1) + 
		    1) * cc_dim1] + cc[ic - 1 + ((k << 1) + 2) * cc_dim1];
	    tr2 = cc[i__ - 1 + ((k << 1) + 1) * cc_dim1] - cc[ic - 1 + ((k << 
		    1) + 2) * cc_dim1];
	    ch[i__ + (k + ch_dim2) * ch_dim1] = cc[i__ + ((k << 1) + 1) * 
		    cc_dim1] - cc[ic + ((k << 1) + 2) * cc_dim1];
	    ti2 = cc[i__ + ((k << 1) + 1) * cc_dim1] + cc[ic + ((k << 1) + 2) 
		    * cc_dim1];
	    ch[i__ - 1 + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i__ - 2] * tr2 
		    - wa1[i__ - 1] * ti2;
	    ch[i__ + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i__ - 2] * ti2 + 
		    wa1[i__ - 1] * tr2;
/* L103: */
	}
/* L104: */
    }
    if (*ido % 2 == 1) {
	return 0;
    }
L105:
    i__1 = *l1;
    for (k = 1; k <= i__1; ++k) {
	ch[*ido + (k + ch_dim2) * ch_dim1] = cc[*ido + ((k << 1) + 1) * 
		cc_dim1] + cc[*ido + ((k << 1) + 1) * cc_dim1];
	ch[*ido + (k + (ch_dim2 << 1)) * ch_dim1] = -(cc[((k << 1) + 2) * 
		cc_dim1 + 1] + cc[((k << 1) + 2) * cc_dim1 + 1]);
/* L106: */
    }
L107:
    return 0;


}

int FFT::radb3_(int *ido, int *l1, double *cc, double *ch, double *wa1, double *wa2)
{
	    /* Initialized data */

    static double taur = -.5f;
    static double taui = .866025403784439f;

    /* System generated locals */
    int cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i__1, i__2;