comparison 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