📄 fst_semi_memo.cpp
字号:
// FST_semi_memo.cpp: implementation of the FST_semi_memo class.
//
//////////////////////////////////////////////////////////////////////
#include "stdafx.h"
#include "FST_semi_memo.h"
#include <math.h>
#include <stdio.h>
#define M_PI 3.14159265358979323846
#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif
//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
FST_semi_memo::FST_semi_memo()
{
}
FST_semi_memo::~FST_semi_memo()
{
}
void FST_semi_memo::FST_semi_memo1(double *rdata, double *idata, double *rcoeffs, double *icoeffs, int size, double **seminaive_naive_table, double *workspace, int dataformat, int cutoff)
{
int bw, m, i, j;
double *rres, *ires;
double *rdataptr, *idataptr;
double *fltres, *scratchpad;
double *eval_pts;
int tmpindex[2];
double pow_one;
double *cos_even;
int l, dummy ;
double tmpA, tmpB ;
bw = size/2;
/* assign space */
cos_even = (double *) malloc(sizeof(double) * bw);
rres = workspace; /* needs (size * size) = (4 * bw^2) */
ires = rres + (size * size); /* needs (size * size) = (4 * bw^2) */
fltres = ires + (size * size); /* needs bw */
eval_pts = fltres + bw; /* needs (2*bw) */
scratchpad = eval_pts + (2*bw); /* needs (24 * bw) */
/* total workspace is (8 * bw^2) + (29 * bw) */
/* do the FFTs along phi */
fftg1.grid_fourier(rdata, idata, rres, ires, size, scratchpad);
/* point to start of output data buffers */
rdataptr = rcoeffs;
idataptr = icoeffs;
for (m=0; m<bw; m++) {
/*
fprintf(stderr,"m = %d\n",m);
*/
/*** test to see if before cutoff or after ***/
if (m < cutoff){
/* do the real part */
semi1.SemiNaiveReduced3(rres+(m*size),
bw,
m,
fltres,
seminaive_naive_table[m],
scratchpad,
cos_even);
/* now load real part of coefficients into output space */
memcpy(rdataptr, fltres, sizeof(double) * (bw - m));
rdataptr += bw-m;
/* do imaginary part */
semi1.SemiNaiveReduced3(ires+(m*size),
bw,
m,
fltres,
seminaive_naive_table[m],
scratchpad,
cos_even);
/* now load imaginary part of coefficients into output space */
memcpy(idataptr, fltres, sizeof(double) * (bw - m));
idataptr += bw-m;
}
else{
/* do real part */
naive1.Naive_AnalysisX(rres+(m*size),
bw,
m,
fltres,
seminaive_naive_table[m],
workspace);
memcpy(rdataptr, fltres, sizeof(double) * (bw - m));
rdataptr += bw-m;
/* do imaginary part */
naive1.Naive_AnalysisX(ires+(m*size),
bw,
m,
fltres,
seminaive_naive_table[m],
workspace);
memcpy(idataptr, fltres, sizeof(double) * (bw - m));
idataptr += bw-m;
}
}
/*** now do upper coefficients ****/
/* now if the data is real, we don't have to compute the
coefficients whose order is less than 0, i.e. since
the data is real, we know that
f-hat(l,-m) = (-1)^m * conjugate(f-hat(l,m)),
so use that to get the rest of the coefficients
dataformat =0 -> samples are complex, =1 -> samples real
*/
if( dataformat == 0 ){
/* note that m is greater than bw here, but this is for
purposes of indexing the input data arrays.
The "true" value of m as a parameter for Pml is
size - m */
/* fprintf(stderr,"\n now the higher order terms \n\n"); */
for (m=bw+1; m<size; m++) {
/*
fprintf(stderr,"m = %d\n",-(size-m));
*/
/* do real part */
semi1.SemiNaiveReduced3(rres+(m*size),
bw,
size-m,
fltres,
seminaive_naive_table[size-m],
scratchpad,
cos_even);
/* now load real part of coefficients into output space */
if ((m % 2) != 0) {
for (i=0; i<m-bw; i++)
rdataptr[i] = -fltres[i];
}
else {
memcpy(rdataptr, fltres, sizeof(double) * (m - bw));
}
rdataptr += m-bw;
/* do imaginary part */
semi1.SemiNaiveReduced3(ires+(m*size),
bw,
size-m,
fltres,
seminaive_naive_table[size-m],
scratchpad,
cos_even);
/* now load real part of coefficients into output space */
if ((m % 2) != 0) {
for (i=0; i<m-bw; i++)
idataptr[i] = -fltres[i];
}
else {
memcpy(idataptr, fltres, sizeof(double) * (m - bw));
}
idataptr += m-bw;
}
}
else /**** if the data is real ****/
{
pow_one = 1.0;
for(i = 1; i < bw; i++){
pow_one *= -1.0;
for( j = i; j < bw; j++){
prif.seanindex2(i, j, bw, tmpindex);
rcoeffs[tmpindex[1]] =
pow_one * rcoeffs[tmpindex[0]];
icoeffs[tmpindex[1]] =
-1.0 * pow_one * icoeffs[tmpindex[0]];
}
}
}
free(cos_even);
/*****************
New in 2.6: Need to massage the coefficients one more time,
to get that right factor in there (how embarrassing)
******************/
tmpA = 2. * sqrt( PI );
tmpB = sqrt( 2. * PI );
for(m=0;m<bw;m++)
for(l=m;l<bw;l++)
{
dummy = prif.seanindex(m,l,bw);
if ( m == 0 )
{
rcoeffs[dummy] *= tmpA ;
icoeffs[dummy] *= tmpA ;
}
else
{
rcoeffs[dummy] *= tmpB ;
icoeffs[dummy] *= tmpB ;
}
/* now for the negative-order coefficients */
if ( m != 0 )
{
dummy = prif.seanindex(-m,l,bw);
rcoeffs[dummy] *= tmpB ;
icoeffs[dummy] *= tmpB ;
}
}
}
void FST_semi_memo::InvFST_semi_memo(double *rcoeffs, double *icoeffs, double *rdata, double *idata, int size, double **transpose_seminaive_naive_table, double *workspace, int dataformat, int cutoff)
{
int bw, m, i, n;
double *rdataptr, *idataptr;
double *rfourdata, *ifourdata;
double *rinvfltres, *iminvfltres, *scratchpad;
double *sin_values, *eval_pts;
int l, dummy ;
double tmpA, tmpB ;
FILE *fp;
bw = size/2;
/* allocate space */
rfourdata = workspace; /* needs (size * size) */
ifourdata = rfourdata + (size * size); /* needs (size * size) */
rinvfltres = ifourdata + (size * size); /* needs (2 * bw) */
iminvfltres = rinvfltres + (2 * bw); /* needs (2 * bw) */
sin_values = iminvfltres + (2 * bw); /* needs (2 * bw) */
eval_pts = sin_values + (2 * bw); /* needs (2 * bw) */
scratchpad = eval_pts + (2 * bw); /* needs (24 * bw) */
/* total workspace = (8 * bw^2) + (32 * bw) */
/* load up the sin_values array */
n = 2*bw;
prim1.ArcCosEvalPts(n, eval_pts);
for (i=0; i<n; i++)
sin_values[i] = sin(eval_pts[i]);
/**********************
New in 2.6: Need to massage the coefficients, to get that
right factor in there (how embarrassing)
***********************/
tmpA = 1./(2. * sqrt(PI) );
tmpB = 1./sqrt(2. * PI ) ;
for(m=0;m<bw;m++)
for(l=m;l<bw;l++)
{
dummy = prif.seanindex(m,l,bw);
if ( m == 0 )
{
rcoeffs[dummy] *= tmpA ;
icoeffs[dummy] *= tmpA ;
}
else
{
rcoeffs[dummy] *= tmpB ;
icoeffs[dummy] *= tmpB ;
}
/* now for the negative-order coefficients */
if ( m != 0 )
{
dummy = prif.seanindex(-m,l,bw);
rcoeffs[dummy] *= tmpB ;
icoeffs[dummy] *= tmpB ;
}
}
/* Now do all of the inverse Legendre transforms */
rdataptr = rcoeffs;
idataptr = icoeffs;
for (m=0; m<bw; m++)
{
/*
fprintf(stderr,"m = %d\n",m);
*/
if(m < cutoff)
{
/* do real part first */
semi1.InvSemiNaiveReduced(rdataptr,
bw,
m,
rinvfltres,
transpose_seminaive_naive_table[m],
sin_values,
scratchpad);
/* now do imaginary part */
semi1.InvSemiNaiveReduced(idataptr,
bw,
m,
iminvfltres,
transpose_seminaive_naive_table[m],
sin_values,
scratchpad);
/* will store normal, then tranpose before doing inverse fft */
memcpy(rfourdata+(m*size), rinvfltres, sizeof(double) * size);
memcpy(ifourdata+(m*size), iminvfltres, sizeof(double) * size);
/* move to next set of coeffs */
rdataptr += bw-m;
idataptr += bw-m;
}
else
{
/* first do the real part */
naive1.Naive_SynthesizeX(rdataptr,
bw,
m,
rinvfltres,
transpose_seminaive_naive_table[m]);
/* now do the imaginary */
naive1.Naive_SynthesizeX(idataptr,
bw,
m,
iminvfltres,
transpose_seminaive_naive_table[m]);
/* will store normal, then tranpose before doing inverse fft */
memcpy(rfourdata+(m*size), rinvfltres, sizeof(double) * size);
memcpy(ifourdata+(m*size), iminvfltres, sizeof(double) * size);
/* move to next set of coeffs */
rdataptr += bw-m;
idataptr += bw-m;
}
} /* closes m loop */
/* now fill in zero values where m = bw (from problem definition) */
memset(rfourdata + (bw * size), 0, sizeof(double) * size);
memset(ifourdata + (bw * size), 0, sizeof(double) * size);
/* now if the data is real, we don't have to compute the
coefficients whose order is less than 0, i.e. since
the data is real, we know that
invf-hat(l,-m) = conjugate(invf-hat(l,m)),
so use that to get the rest of the real data
dataformat =0 -> samples are complex, =1 -> samples real
*/
if(dataformat == 0){
/* now do negative m values */
for (m=bw+1; m<size; m++)
{
/*
fprintf(stderr,"m = %d\n",-(size-m));
*/
/* do real part first */
semi1.InvSemiNaiveReduced(rdataptr,
bw,
size - m,
rinvfltres,
transpose_seminaive_naive_table[size - m],
sin_values,
scratchpad);
/* now do imaginary part */
semi1.InvSemiNaiveReduced(idataptr,
bw,
size - m,
iminvfltres,
transpose_seminaive_naive_table[size - m],
sin_values,
scratchpad);
/* will store normal, then tranpose before doing inverse fft */
if ((m % 2) != 0)
for(i=0; i< size; i++){
rinvfltres[i] = -rinvfltres[i];
iminvfltres[i] = -iminvfltres[i];
}
memcpy(rfourdata + (m*size), rinvfltres, sizeof(double) * size);
memcpy(ifourdata + (m*size), iminvfltres, sizeof(double) * size);
/* move to next set of coeffs */
rdataptr += bw-(size-m);
idataptr += bw-(size-m);
} /* closes m loop */
}
else {
for(m = bw + 1; m < size; m++){
memcpy(rfourdata+(m*size), rfourdata+((size-m)*size),
sizeof(double) * size);
memcpy(ifourdata+(m*size), ifourdata+((size-m)*size),
sizeof(double) * size);
for(i = 0; i < size; i++)
ifourdata[(m*size)+i] *= -1.0;
}
}
/** now transpose **/
prif.transpose(rfourdata, size);
prif.transpose(ifourdata, size);
/* now do inverse fourier grid computation */
fftg1.grid_invfourier(rfourdata, ifourdata, rdata, idata, size, scratchpad);
/* amscray */
}
void FST_semi_memo::FZT_semi_memo(double *rdata, double *idata, double *rres, double *ires, int size, double *cos_pml_table, double *workspace, int dataformat)
{
int i, j, bw;
double *r0, *i0, dsize;
double tmpreal, tmpimag;
double *scratchpad;
double *cos_even;
int l, dummy ;
double tmpA ;
bw = size/2;
/* assign memory */
r0 = workspace; /* needs (2 * bw) */
i0 = r0 + (2 * bw); /* needs (2 * bw) */
cos_even = i0 + (2 * bw); /* needs (1 * bw) */
scratchpad = cos_even + (1 * bw); /* needs (8 * bw) */
/* total workspace = 13*bw */
dsize = (double) size;
/* compute the m = 0 components */
for (i=0; i<size; i++) {
tmpreal = 0.0;
tmpimag = 0.0;
for(j=0; j<size; j++) {
tmpreal += rdata[(i*size)+j];
tmpimag += idata[(i*size)+j];
}
/* normalize */
r0[i] = tmpreal/dsize;
i0[i] = tmpimag/dsize;
}
/* do the real part */
semi1.SemiNaiveReduced3(r0,
bw,
0,
rres,
cos_pml_table,
scratchpad,
cos_even);
if(dataformat == 0) /* do imaginary part */
semi1.SemiNaiveReduced3(i0,
bw,
0,
ires,
cos_pml_table,
scratchpad,
cos_even);
else /* otherwise set coefficients = 0 */
memset(ires, 0, sizeof(double) * size);
/*****************
New in 2.6: Need to massage the coefficients one more time,
to get that right factor in there (how embarrassing)
******************/
tmpA = 2. * sqrt( PI );
for(l=0;l<bw;l++)
{
dummy = prif.seanindex(0,l,bw);
rres[dummy] *= tmpA ;
ires[dummy] *= tmpA ;
}
}
void FST_semi_memo::TransMult(double *rdatacoeffs, double *idatacoeffs, double *rfiltercoeffs, double *ifiltercoeffs, double *rres, double *ires, int size)
{
int m, l, bw;
double *rdptr, *idptr, *rrptr, *irptr;
bw = size/2;
rdptr = rdatacoeffs;
idptr = idatacoeffs;
rrptr = rres;
irptr = ires;
for (m=0; m<bw; m++) {
for (l=m; l<bw; l++) {
compmult(rfiltercoeffs[l], ifiltercoeffs[l],
rdptr[l-m], idptr[l-m],
rrptr[l-m], irptr[l-m]);
rrptr[l-m] *= sqrt(4*M_PI/(2*l+1));
irptr[l-m] *= sqrt(4*M_PI/(2*l+1));
}
rdptr += bw-m; idptr += bw-m;
rrptr += bw-m; irptr += bw-m;
}
for (m=bw+1; m<size; m++) {
for (l=size-m; l<bw; l++){
compmult(rfiltercoeffs[l], ifiltercoeffs[l],
rdptr[l-size+m], idptr[l-size+m],
rrptr[l-size+m], irptr[l-size+m]);
rrptr[l-size+m] *= sqrt(4*M_PI/(2*l+1));
irptr[l-size+m] *= sqrt(4*M_PI/(2*l+1));
}
rdptr += m-bw; idptr += m-bw;
rrptr += m-bw; irptr += m-bw;
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -