fftcode.cpp

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

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

#include "stdafx.h"
#include "FFTcode.h"

#include <stdio.h>
#include <string.h>  /* for memcpy */

#define compmult(a,b,c,d,e,f) (e) = ((a)*(c))-((b)*(d)); (f) = ((a)*(d))+((b)*(c))
#define compadd(a,b,c,d,e,f) (e) = (a) + (c); (f) = (b) + (d)



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

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

FFTcode::FFTcode()
{

}

FFTcode::~FFTcode()
{

}


int FFTcode::ilog2fast(int value)
{
	  switch (value)
    {
    case 1: return 0;
    case 2: return 1;
    case 4: return 2;
    case 8: return 3;
    case 16: return 4;
    case 32: return 5;
    case 64: return 6;
    case 128: return 7;
    case 256: return 8;
    case 512: return 9;
    case 1024: return 10;
    case 2048: return 11;
    case 4096: return 12;
    case 8192: return 13;
    case 16384: return 14;
    case 32768: return 15;
    case 65536: return 16;
    default: return -1;
    }


}



void FFTcode::FFTEval(double *reals, double *imags, double *r_vals, double *i_vals, int N, int K, double *workspace, int brflag)
{
	  const double *rootptr;
  double *rcptr, *icptr, *rrptr, *irptr;
  double rtemp, itemp;
  int rem2, j, k, l, toggle;

  double tmproot0, tmproot1;

  memcpy(workspace, reals, sizeof(double) * N);
  memcpy(workspace+K, imags, sizeof(double) * N);

  rootptr = r4096; /* point to roots of unity in BR order */
  rem2 = N/2; /* current size of remainders */
  rcptr = workspace;  /* points to real part of current poly coeff */
  icptr = workspace+K;  /* points to imag part of current poly coeff */
  rrptr = r_vals; /* points to real part of current remainder coeff */
  irptr = i_vals; /* points to imag part of current remainder coeff */
  toggle = 0;

  for (l=0; l < ilog2fast(N); l++)
    {
      for (k=0; k<K; k+=rem2) 
	{
	  tmproot0 = rootptr[0]; tmproot1 = rootptr[1];
	  for (j=0; j<rem2; j++) 
	    {
	      compmult(tmproot0,tmproot1,rcptr[rem2+j],icptr[rem2+j],
		       rtemp,itemp);
	      compadd(rtemp,itemp,rcptr[j],icptr[j],
		      rrptr[j],irptr[j]);
	    }
	  rrptr += rem2;
	  irptr += rem2;
	  rootptr += 2;
	  toggle++;
	  if ((toggle % 2) == 0) 
	    {
	      /* check l==0 here. If l==0, then we need to
		 fake it into thinking that there are multiple
		 copies of the input coeffs */
	      if (l != 0) 
		{ 
		  rcptr += 2*rem2;
		  icptr += 2*rem2;
		}
	      else 
		{
		  rcptr = workspace;
		  icptr = workspace+K;
		}
	    }
	} /* closes k loop - done with one level */
    
      if ((l % 2) == 0) 
	{
	  rcptr = r_vals;
	  icptr = i_vals;
	  rrptr = workspace;
	  irptr = workspace + K;
	}
      else 
	{
	  rcptr = workspace;
	  icptr = workspace + K;
	  rrptr = r_vals;
	  irptr = i_vals;
	}
    
      rem2 /= 2;
      rootptr = r4096;
    }
  if ((l % 2) == 0) 
    {
      memcpy(r_vals, workspace, sizeof(double) * K);
      memcpy(i_vals, workspace+K, sizeof(double) * K);
    }
  
  if (brflag == 1) { /* bitreverse the data */
    ind.bitreverse(r_vals,K, workspace);
    ind.bitreverse(i_vals,K, workspace);
  }


}



void FFTcode::FFTInterp(double *reals, double *imags, double *r_vals, double *i_vals, int N, int K, double *workspace, int brflag)
{
	  const double *rootptr;
  double *rhptr, *ihptr, *rlptr, *ilptr;
  double rtemp0, itemp0, rtemp1, itemp1, dn;
  int degree, j, k, l, toggle, upperlim;
  double rptr0, rptr1, rptr2, rptr3;
  double tmprval, tmpival;

  dn = (double) N;

  /* normalize and permute if necessary */
  for (j=0; j<N; j++)
    {
      workspace[j] = reals[j]/dn;
      workspace[j+N] = imags[j]/dn;
    }
  if (brflag == 1)
    {
      ind.bitreverse(workspace, N, workspace+(2*N));
      ind.bitreverse(workspace+N, N, workspace+(2*N));
    }

  rootptr = r4096; /* point to roots of unity in BR order */
  degree = 1; /* next order of polynomials - stride length */
  rlptr = workspace;  /* points to real part of lower order polys */
  ilptr = workspace+N;  /* points to imag part of lower order polys */
  rhptr = r_vals; /* points to real part of higher order polys */
  ihptr = i_vals; /* points to imag part of higher order polys */
  toggle = 0;

  for (l=1; l <= ilog2fast(K); l++)
    {
      upperlim = N/(2*degree);
      for (k=0; k < upperlim; k++)
	{
	  for (j=0; j < degree; j++)
	    {
	      compadd(rlptr[j],ilptr[j],rlptr[j+degree],ilptr[j+degree],
		      rhptr[j],ihptr[j]);
	    }
	  rhptr += degree;
	  ihptr += degree;

	  rptr0 = rootptr[0]; rptr1 = rootptr[1];
	  rptr2 = rootptr[2]; rptr3 = rootptr[3];

	  for (j=0; j < degree; j++)
	    {

	      compmult(rptr0,-rptr1,rlptr[j],ilptr[j],
		       rtemp0,itemp0);
	      compmult(rptr2,-rptr3,rlptr[j+degree],ilptr[j+degree],
		       rtemp1,itemp1);
	      compadd(rtemp0,itemp0,rtemp1,itemp1,
		      rhptr[j],ihptr[j]);

	    }

	  rhptr += degree;
	  ihptr += degree;
	  rlptr += (2*degree);
	  ilptr += (2*degree);
	  rootptr += 4;
	} /* end of k loop */
    
      if ((l % 2) != 0)
	{
	  rlptr = r_vals;
	  ilptr = i_vals;
	  rhptr = workspace;
	  ihptr = workspace + N;
	}
      else
	{
	  rhptr = r_vals;
	  ihptr = i_vals;
	  rlptr = workspace;
	  ilptr = workspace + N;
	}
      degree *= 2;
      rootptr = r4096;
      toggle++;
    
    } /* end of l loop */
  
  if ((toggle % 2) == 0)
    {
      memcpy(r_vals, workspace, sizeof(double) * N);
      memcpy(i_vals, workspace+N, sizeof(double) * N);
    }
  
  for (j=0; j<K; j++)
    {
      tmprval = r_vals[j]; tmpival = i_vals[j];
      for (k=K; k<N; k+=K)
	{
	  tmprval += r_vals[j+k];
	  tmpival += i_vals[j+k];
	}
      r_vals[j] = tmprval; i_vals[j] = tmpival;
    }
  


}