annotate windowfunctions.cpp @ 2:fc19d45615d1

* Make all file names lower-case to avoid case ambiguity
(some includes differed in case from the filenames they were
trying to include).  Also replace MinGW-specific mem.h with
string.h

//---------------------------------------------------------------------------
#include <math.h>
#include <string.h>
#include "windowfunctions.h"
#include "align8.h"
//---------------------------------------------------------------------------

/*
  function I0: Bessel function of order zero

  Returns the Bessel function of x.
*/
double I0(double x)
{
  const double eps=1e-9;

  int n=1;
  double S=1, D=1, T;

  while (D>eps*S)
  {
    T=x/(2*n++);
    D*=T*T;
    S+=D;
  }
  return S;
}//I0

/*
  function FillWindow: fills a buffer $Result with a window function.

  In: wt:window type
      Count: window size
      ips & ups: extra window specificatin parameters
  Out: newwindow[Count+1]: the window function.

  No return value. This function is designed assuming Count being even, and Count/2 is the symmetry
  centre of newwindow. newwindow has physical size Count+1 for compatibility purpose. For all vanishing
  windows (e.g. Hann) newwindow[0]=newwindow[Count]=0.
*/
void FillWindow(double* Result, WindowType wt, int Count, int* ips, double* dps)
{
	switch(wt)
  {
    case wtRectangle:
      for (int i=0; i<=Count; i++)
        Result[i]=1;
      break;
    case wtTriangular:
      for (int i=0; i<=Count; i++)
				Result[i]=1-2*fabs(i-Count/2.0)/Count;
      break;
    case wtHamming:
      for (int i=0; i<=Count; i++)
        Result[i]=(0.54-0.46*cos(2*M_PI*i/Count));
      break;
    case wtBlackman:
      for (int i=0; i<=Count; i++)
        Result[i]=0.42-0.5*cos(2*M_PI*i/Count)+0.08*cos(4*M_PI*i/Count);
      break;
    case wtGaussian:
    {
      double num=*dps*M_PI*M_PI/2/Count/Count;
      for (int i=0; i<=Count; i++)
        Result[i]=exp(-(i-Count/2.0)*(i-Count/2.0)*num);
      break;
    }
    case wtKaiser:
    {
      double ldps=*dps*M_PI*M_PI/2/Count;
      double den=I0(0.5*Count*ldps);
      for (int i=0; i<=Count; i++)
        Result[i]=I0(ldps*sqrt(0.25*Count*Count-(i-Count*0.5)*(i-Count*0.5)))/den;
      break;
    }
    case wtHalfCosine:
    {
      for (int i=0; i<=Count; i++)
				Result[i]=sin(M_PI*i/Count);
      break;
    }
    case wtHann:
    {
      for (int i=0; i<=Count; i++)
        Result[i]=(0.5-cos(M_PI*2*i/Count)/2);
      break;
    }
    case wtHannSqr:
    {
      for (int i=0; i<=Count; i++)
        Result[i]=(3-4*cos(M_PI*2*i/Count)+cos(M_PI*4*i/Count))/8;
      break;
		}
		case wtHann3sqrt:
		{
			for (int i=0; i<=Count; i++)
			{
				double tmp=sin(M_PI*i/Count);
				Result[i]=tmp*tmp*tmp;
			}
			break;
		}
	}
}//FillWindow

/*
  function NewWindow: creates a window function.

  In: wt: window type
      Count: window size
      ips & ups: extra window specificatin parameters
  Out: newwindow[Count+1]: the window function.

  Returns pointer to newwindow. newwindow is created anew if newwindow=0 is specified on start.
*/
double* NewWindow(WindowType wt, int Count, int* ips, double* dps, double* newwindow)
{
	double* Result=newwindow; if (!Result) Result=new double[Count+1];
	FillWindow(Result, wt, Count, ips, dps);
	return Result;
}//NewWindow

/*
  function NewWindow8: 8-byte aligned version of NewWindow.

  In: wt: window type
      Count: window size
      ips & ups: extra window specificatin parameters
  Out: newwindow[Count+1]: the window function.

  Returns pointer to newwindow. newwindow is created anew and 8-byte aligned if newwindow=0 is
  specified on start. However, if newwindow is not 8-byte aligned on start, the unaligned buffer is
  returned.
*/
double* NewWindow8(WindowType wt, int Count, int* ips, double* dps, double* newwindow)
{
	double* Result=newwindow; if (!Result) Result=(double*)malloc8(sizeof(double)*(Count+1));
	FillWindow(Result, wt, Count, ips, dps);
	return Result;
}//NewWindow8

/*
  function NewdWindow: computes the derivative of a window function.

  In: wt: window type
      Count: window size
      ips & ups: extra window specificatin parameters
  Out: newdwindow[Count+1]: the derivative window function.

  Returns pointer to newdwindow. newdwindow is created anew if newdwindow=0 is specified on start.
*/
double* NewdWindow(WindowType wt, int Count, int* ips, double* dps, double* newdwindow)
{
	double* Result=newdwindow; if (!Result) Result=new double[Count+1];
	double piiCount=M_PI/Count;
	switch(wt)
	{
		case wtRectangle:
			memset(Result, 0, sizeof(double)*(Count+1));
			break;
    case wtTriangular:
      throw("Trying to differentiate triangular window.");
      break;
		case wtHamming:
			for (int i=0; i<=Count; i++)
				Result[i]=0.92*piiCount*sin(2*piiCount*i);
			break;
		case wtBlackman:
			for (int i=0; i<=Count; i++)
				Result[i]=piiCount*sin(2*piiCount*i)-0.32*piiCount*sin(4*piiCount*i);
			break;
    case wtGaussian:
      throw("Trying to differentiate Gaussian window.");
      break;
    case wtKaiser:
      throw("Trying to differentiate Kaiser window.");
      break;
    case wtHalfCosine:
		{
			for (int i=0; i<=Count; i++)
				Result[i]=piiCount*cos(piiCount*i);
			break;
		}
		case wtHann:
		{
			for (int i=0; i<=Count; i++)
				Result[i]=piiCount*sin(2*piiCount*i);
			break;
		}
		case wtHannSqr:
		{
			for (int i=0; i<=Count; i++)
				Result[i]=piiCount*sin(2*piiCount*i)-0.5*piiCount*sin(4*piiCount*i);
			break;
		}
		case wtHann3sqrt:
		{
			for (int i=0; i<=Count; i++)
			{
				double s=sin(M_PI*i/Count), c=cos(M_PI*i/Count);
				Result[i]=3*piiCount*s*s*c;
			}
		}
	}
	return Result;
}//NewdWindow

/*
  function NewddWindow: computes the 2nd-order derivative of a window function.

  In: wt: window type
      Count: window size
      ips & ups: extra window specificatin parameters
  Out: newddwindow[Count+1]: the 2nd-order derivative window function.

  Returns pointer to newddwindow. newddwindow is created anew if newddwindow=0 is specified on start.
*/
double* NewddWindow(WindowType wt, int Count, int* ips, double* dps, double* newddwindow)
{
	double* Result=newddwindow; if (!Result) Result=new double[Count+1];
	double piiCount=M_PI/Count;
	double piC2=piiCount*piiCount;
	switch(wt)
	{
		case wtRectangle:
			memset(Result, 0, sizeof(double)*(Count+1));
      break;
    case wtTriangular:
      throw("Trying to double-differentiate triangular window.");
      break;
    case wtHamming:
			for (int i=0; i<=Count; i++)
				Result[i]=1.84*piC2*cos(2*piiCount*i);
			break;
		case wtBlackman:
			for (int i=0; i<=Count; i++)
				Result[i]=2*piC2*cos(2*piiCount*i)-1.28*piC2*cos(4*piiCount*i);
			break;
    case wtGaussian:
      throw("Trying to double-differentiate Gaussian window.");
      break;
    case wtKaiser:
      throw("Trying to double-differentiate Kaiser window.");
      break;
    case wtHalfCosine:
		{
			for (int i=0; i<=Count; i++)
				Result[i]=-piC2*sin(piiCount*i);
			break;
		}
		case wtHann:
		{
			for (int i=0; i<=Count; i++)
				Result[i]=2*piC2*cos(2*piiCount*i);
			break;
		}
		case wtHannSqr:
		{
			for (int i=0; i<=Count; i++)
				Result[i]=2*piC2*cos(2*piiCount*i)-2*piC2*cos(4*piiCount*i);
			break;
		}
		case wtHann3sqrt:
		{
			for (int i=0; i<=Count; i++)
			{
				double s=sin(M_PI*i/Count), c=cos(M_PI*i/Count);
				Result[i]=3*piC2*(2*c*c-s*s)*s;
			}
			break;
		}
	}
	return Result;
}//NewddWindow

/*
  function NewdddWindow: computes the 3rd-order derivative of a window function.
  In: wt: window type
      Count: window size
      ips & ups: extra window specificatin parameters
  Out: newdddwindow[Count+1]: the 3rd-order derivative window function.

  Returns pointer to newdddwindow. newdddwindow is created anew if newdddwindow=0 is specified on start.
*/
double* NewdddWindow(WindowType wt, int Count, int* ips, double* dps, double* newdddwindow)
{
	double* Result=newdddwindow; if (!Result) Result=new double[Count+1];
	double piiCount=M_PI/Count;
	double piC2=piiCount*piiCount;
	double piC3=piiCount*piC2;
	switch(wt)
	{
    case wtRectangle:
      memset(Result, 0, sizeof(double)*(Count+1));
      break;
    case wtTriangular:
      throw("Trying to triple-differentiate triangular window.");
      break;
    case wtHamming:
			for (int i=0; i<=Count; i++)
				Result[i]=-3.68*piC3*sin(2*piiCount*i);
			break;
		case wtBlackman:
			for (int i=0; i<=Count; i++)
				Result[i]=-4*piC3*sin(2*piiCount*i)+5.12*piC3*sin(4*piiCount*i);
			break;
    case wtGaussian:
      throw("Trying to triple-differentiate Gaussian window.");
      break;
    case wtKaiser:
      throw("Trying to triple-differentiate Kaiser window.");
      break;
    case wtHalfCosine:
		{
			for (int i=0; i<=Count; i++)
				Result[i]=-piC3*cos(piiCount*i);
			break;
		}
		case wtHann:
		{
			for (int i=0; i<=Count; i++)
				Result[i]=-4*piC3*sin(2*piiCount*i);
			break;
    }
		case wtHannSqr:
    {
      for (int i=0; i<=Count; i++)
				Result[i]=-4*piC3*sin(2*piiCount*i)+8*piC3*sin(4*piiCount*i);
			break;
    }
    case wtHann3sqrt:
      throw("Trying to triple-differentiate Hann^1.5 window.");
      break;
	}
  return Result;
}//NewdddWindow

//---------------------------------------------------------------------------
/*
  function windowspec: computes a few descriptors of a cosine family window. A window function in the
  cosine window family is the linear combination of a few cosine functions, therefore has a cosine
  decomposition.

  In: wt: window type
      N: window size
  Out: c[M+1], coefficients of cosine decomposition
       iH2: reciprocal of square of L2 norm,
       d[2M+1]: self-convolution of c[], optional

  No return value.
*/
void windowspec(WindowType wt, int N, int* M, double* c, double* iH2, double* d)
{
  switch(wt)
  {
    case wtRectangle:
      M[0]=0, c[0]=1, iH2[0]=1.0/N/N;
      if (d) d[0]=1;
      break;
    case wtHamming:
      M[0]=1, c[0]=0.54, c[1]=0.23, iH2[0]=1.0/(0.3974*N*N);
      if (d) d[0]=0.3974, d[1]=0.2484, d[2]=0.0529;
      break;
    case wtBlackman:
      M[0]=2, c[0]=0.42, c[1]=0.25, c[2]=0.04, iH2[0]=1.0/(0.3046*N*N);
      if (d) d[0]=0.3046, d[1]=0.23, d[2]=0.0961, d[4]=0.02, d[5]=0.0016;
      break;
    case wtHann:
      M[0]=1, c[0]=0.5, c[1]=0.25, iH2[0]=1.0/(0.375*N*N);
      if (d) d[0]=0.375, d[1]=0.25, d[2]=0.0625;
      break;
    default:
      M[0]=1, c[0]=0.5, c[1]=0.25, iH2[0]=1.0/(0.375*N*N);
      if (d) d[0]=0.375, d[1]=0.25, d[2]=0.0625;
      break;
  }
}//windowspec