naive_synthesis.cpp

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

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

#include "stdafx.h"
#include "naive_synthesis.h"

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#ifndef PI
#define PI 3.14159265358979
#endif

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

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

naive_synthesis::naive_synthesis()
{

}

naive_synthesis::~naive_synthesis()
{

}

double naive_synthesis::ns_L2_an(int m, int l)
{
	return (sqrt((((double) (2*l+3))/((double) (2*l+1))) *
	       (((double) (l-m+1))/((double) (l+m+1)))) *
	  (((double) (2*l+1))/((double) (l-m+1))));


}


double naive_synthesis::ns_L2_bn()
{
	return ((double) 0.0);


}


double naive_synthesis::ns_L2_cn(int m, int l)
{
	if (l != 0) {
    return (-1.0 *
	  sqrt((((double) (2*l+3))/((double) (2*l-1))) *
	       (((double) (l-m+1))/((double) (l+m+1))) *
	       (((double) (l-m))/((double) (l+m)))) *
	  (((double) (l+m))/((double) (l-m+1))));
  }
  else
    return 0.0;


}


void naive_synthesis::ns_vec_add(double *data1, double *data2, double *result, int n)
{
	 int k;

  for(k = 0; k < n % 4; ++k)
    {
      result[k] = data1[k] + data2[k];
    }

  for( ; k < n; k += 4)
    {
      result[k] = data1[k] + data2[k];
      result[k + 1] = data1[k + 1] + data2[k + 1];
      result[k + 2] = data1[k + 2] + data2[k + 2];
      result[k + 3] = data1[k + 3] + data2[k + 3];
    }


}


void naive_synthesis::ns_vec_mul(double scalar, double *data1, double *result, int n)
{
	   int k;


   for( k = 0; k < n % 4; ++k)
     result[k] = scalar * data1[k];

   for( ; k < n; k += 4)
     {
      result[k] = scalar * data1[k];
      result[k + 1] = scalar * data1[k + 1];
      result[k + 2] = scalar * data1[k + 2];
      result[k + 3] = scalar * data1[k + 3];
    }


}



void naive_synthesis::ns_vec_pt_mul(const double *data1, const double *data2, double *result, int n)
{
	  int k;

   for (k = 0; k < n % 4; ++k)
     result[k] = data1[k] * data2[k];

   for ( ; k < n; k += 4)
     {
       result[k] = data1[k] * data2[k];
       result[k + 1] = data1[k + 1] * data2[k + 1];
       result[k + 2] = data1[k + 2] * data2[k + 2];
       result[k + 3] = data1[k + 3] * data2[k + 3];
     }


}


void naive_synthesis::ns_ArcCosEvalPts(int n, double *eval_pts)
{
	int i;
    double twoN;

    twoN = (double) (2 * n);

   for (i=0; i<n; i++)
     eval_pts[i] = (( 2.0*((double)i)+1.0 ) * PI) / twoN;


}


void naive_synthesis::ns_EvalPts(int n, double *eval_pts)
{
	 int i;
    double twoN;

    twoN = (double) (2*n);

   for (i=0; i<n; i++)
     eval_pts[i] = cos((( 2.0*((double)i)+1.0 ) * PI) / twoN);


}

void naive_synthesis::ns_Pmm_L2(int m, double *eval_pts, int n, double *result)

{
  int i;
  double md, id, mcons;

  id = (double) 0.0;
  md = (double) m;
  mcons = sqrt(md + 0.5);

  for (i=0; i<m; i++) {
    mcons *= sqrt((md-(id/2.0))/(md-id));
    id += 1.0;
  }
  if (m != 0 )
    mcons *= pow(2.0,-md/2.0);
  if ((m % 2) != 0) mcons *= -1.0;

  for (i=0; i<n; i++) 
    result[i] = mcons * pow(sin(eval_pts[i]),((double) m));

}

void naive_synthesis::Naive_Synthesize1(int bw, int m, double *coeffs, double *result, double *workspace)
{
    double *prev, *prevprev, *temp1, *temp2, *temp3, *temp4, *x_i, *eval_args;
    double t1;
    int i, j;

    prevprev = workspace;
    prev = prevprev + (2*bw);
    temp1 = prev + (2*bw);
    temp2 = temp1 + (2*bw);
    temp3 = temp2 + (2*bw);
    temp4 = temp3 + (2*bw);
    x_i = temp4 + (2*bw);
    eval_args = x_i + (2*bw);


    /* get the evaluation nodes */
    ns_EvalPts(2*bw,x_i);
    ns_ArcCosEvalPts(2*bw,eval_args);
    

    /* set initial values of first two Pmls */
    for (i=0; i<2*bw; i++) 
      prevprev[i] = 0.0;
    if (m == 0) {
	for (i=0; i<2*bw; i++) {
	    prev[i] = 1.0;
	}
    }
    else 
      ns_Pmm_L2(m, eval_args, 2*bw, prev);

    /* this may be useful for generating Gmls later but
       removed for now

    if ((m % 2) == 1) { 
	for (i=0; i<2*bw; i++)
	  prev[i] /= sin(eval_args[i]);
    }
    */

    /* make sure result is zeroed out */
    for (i=0; i<(2*bw); i++)
      result[i] = 0.0;

    /* add in Pmm contribution */

    if (coeffs[0] != 0.0)
      {
	t1 = coeffs[0];
	for (j=0; j<(2*bw); j++)
	  result[j] += t1 * prev[j];
      }

    /* now generate remaining pmls while synthesizing function */

    for (i=0; i<bw-m-1; i++) {
	ns_vec_mul(ns_L2_cn(m,m+i),prevprev,temp1,2*bw);
	ns_vec_pt_mul(prev, x_i, temp2, 2*bw);
	ns_vec_mul(ns_L2_an(m,m+i), temp2, temp3, 2*bw);
	ns_vec_add(temp3, temp1, temp4, 2*bw); /* temp4 now contains P(m,m+i+1) */

	/* add in contribution */
	if (coeffs[i+1] != 0.0)
	  {
	    t1 = coeffs[i+1];
	    for (j=0; j<(2*bw); j++)
	      result[j] += (t1 * temp4[j]);
	  }

	/* now update Pi and P(i+1) */
	memcpy(prevprev, prev, (size_t) sizeof(double) * 2 * bw);
	memcpy(prev, temp4, (size_t) sizeof(double) * 2 * bw);

    }

}






void naive_synthesis::Naive_Analysis(int bw, int m, double *data, double *result, double *workspace)
{
	   double *prev, *prevprev, *temp1, *temp2, *temp3, *temp4, *x_i, *eval_args;
    double *wdata;
    const double *weight_vec;
    double t1;
    int i, j;

    prevprev = workspace;
    prev = prevprev + (2*bw);
    temp1 = prev + (2*bw);
    temp2 = temp1 + (2*bw);
    temp3 = temp2 + (2*bw);
    temp4 = temp3 + (2*bw);
    x_i = temp4 + (2*bw);
    eval_args = x_i + (2*bw);
    wdata = eval_args + (2*bw);


    /* get the evaluation nodes */
    ns_EvalPts(2*bw,x_i);
    ns_ArcCosEvalPts(2*bw,eval_args);
    

    /* set initial values of first two Pmls */
    for (i=0; i<2*bw; i++) 
      prevprev[i] = 0.0;
    if (m == 0) {
	for (i=0; i<2*bw; i++) {
	    prev[i] = 0.5;
	}
    }
    else 
      ns_Pmm_L2(m, eval_args, 2*bw, prev);

    /* make sure result is zeroed out */
    for (i=0; i<(bw-m); i++)
      result[i] = 0.0;

    /* apply quadrature weights */
    weight_vec = wei.get_weights(bw);

    for (i=0; i<(2*bw); i++)
      wdata[i] = data[i] * weight_vec[i];


    /* compute Pmm coefficient */
    t1 = 0.0;
    for (j=0; j<(2*bw); j++)
      t1 += wdata[j] * prev[j];
    result[0] = t1;

    /* now generate remaining pmls while computing coefficients */

    for (i=0; i<bw-m-1; i++) {
	ns_vec_mul(ns_L2_cn(m,m+i),prevprev,temp1,2*bw);
	ns_vec_pt_mul(prev, x_i, temp2, 2*bw);
	ns_vec_mul(ns_L2_an(m,m+i), temp2, temp3, 2*bw);
	ns_vec_add(temp3, temp1, temp4, 2*bw); /* temp4 now contains P(m,m+i+1) */
	
	/* compute this coefficient */
	t1 = 0.0;
	for (j=0; j<(2*bw); j++)
	  t1 += wdata[j] * temp4[j];
	result[i+1] = t1;

	/* now update Pi and P(i+1) */
	/***
	for (j=0; j<2*bw; j++) {
	    prevprev[j] = prev[j];
	    prev[j] = temp4[j];
	}
	***/
	memcpy( prevprev, prev, sizeof(double) * 2 * bw );
	memcpy( prev, temp4, sizeof(double) * 2 * bw );
    }