fft_grids.cpp

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

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

#include "stdafx.h"
#include "fft_grids.h"
#include <stdlib.h>

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

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

fft_grids::fft_grids()
{

}

fft_grids::~fft_grids()
{

}

double* fft_grids::precomp_fft(int size)
{
	double *wSave;

  /* allocate space */
  wSave = (double *) malloc( sizeof(double) * (2 * size + 15) );

  /* now precompute data */
  fft2.rffti_( &size, wSave );

  return wSave;
}


void fft_grids::grid_fourier_FFTPACK(double *grid, int size, double *wSave)
{
	int i;
  double tmp;

  /* normalizing factor */
  tmp = 1.0 / ((double) size);

  for (i = 0 ; i < size ; i ++)
    {
      fft2.rfftf_(&size, grid+(i*size), wSave);
    }

  /* normalize result */
  for ( i = 0 ; i < size * size ; i ++ )
    grid[i] *= tmp;

  /* now transpose the result */
  pris.transpose(grid, size);
}


void fft_grids::grid_invfourier_FFTPACK( double *grid, int size, double *wSave )
{
  int i;
  for (i=0; i<size; i++)
    {
      fft2.rfftb_( &size, grid+(i*size), wSave);
    }
}

void fft_grids::grid_fourier(double *realgrid, double *imaggrid, double *rmatrix, double *imatrix, int size, double *workspace)
{
	  double *rout, *iout, *scratchpad;
  int i;

  /* need to assign workspace - need 6*size total */

  rout = workspace; /* needs size space */
  iout = rout + size; /* needs size space */
  scratchpad = iout + size; /* needs 4 * size space */

  if(1)
    {
      for (i=0; i<size; i++) 
	{
//进行一个反的FFT变换，第一个size is the number of samples, 第二个size是 the number of coefficents

	  fftc.FFTInterp(realgrid+(i*size), imaggrid+(i*size),
		    rmatrix+(i*size),
		    imatrix+(i*size),
		    size, size, scratchpad, 1);
	}
    }
  else
    {
      for (i=0; i<size; i++) 
	{
	  fftc.FFTInterp(realgrid+(i*size), imaggrid+(i*size),
		    rout,
		    iout,
		    size, size, scratchpad, 1);

	  memcpy(rmatrix+(i*size),rout, sizeof(double) * size);
	  memcpy(imatrix+(i*size),iout, sizeof(double) * size);

	}
    }

  
  /* now transpose the results */
  pris.transpose(rmatrix,size);
  pris.transpose(imatrix,size);


}

void fft_grids::grid_invfourier(double *realgrid, double *imaggrid, double *rmatrix, double *imatrix, int size, double *workspace)
{
	  double *rout, *iout, *scratchpad;
  int i;

  /* need to assign workspace - need 6*size total */
  
  rout = workspace; /* needs size space */
  iout = rout + size; /* needs size space */
  scratchpad = iout + size; /* needs 4 * size space */
  
  for (i=0; i<size; i++) {
    fftc.FFTEval(realgrid+(i*size), imaggrid+(i*size),
	    rmatrix+(i*size),
	    imatrix+(i*size), size, size, scratchpad, 1);
  }


}