xue@11: /* xue@11: Harmonic sinusoidal modelling and tools xue@11: xue@11: C++ code package for harmonic sinusoidal modelling and relevant signal processing. xue@11: Centre for Digital Music, Queen Mary, University of London. xue@11: This file copyright 2011 Wen Xue. xue@11: xue@11: This program is free software; you can redistribute it and/or xue@11: modify it under the terms of the GNU General Public License as xue@11: published by the Free Software Foundation; either version 2 of the xue@11: License, or (at your option) any later version. xue@11: */ xue@1: //--------------------------------------------------------------------------- xue@1: #include Chris@2: #include Chris@2: #include "windowfunctions.h" xue@1: #include "align8.h" xue@1: //--------------------------------------------------------------------------- xue@1: Chris@5: /** xue@1: function I0: Bessel function of order zero xue@1: xue@1: Returns the Bessel function of x. xue@1: */ xue@1: double I0(double x) xue@1: { xue@1: const double eps=1e-9; xue@1: xue@1: int n=1; xue@1: double S=1, D=1, T; xue@1: xue@1: while (D>eps*S) xue@1: { xue@1: T=x/(2*n++); xue@1: D*=T*T; xue@1: S+=D; xue@1: } xue@1: return S; xue@1: }//I0 xue@1: Chris@5: /** xue@1: function FillWindow: fills a buffer $Result with a window function. xue@1: xue@1: In: wt:window type xue@1: Count: window size xue@1: ips & ups: extra window specificatin parameters xue@1: Out: newwindow[Count+1]: the window function. xue@1: xue@1: No return value. This function is designed assuming Count being even, and Count/2 is the symmetry xue@1: centre of newwindow. newwindow has physical size Count+1 for compatibility purpose. For all vanishing xue@1: windows (e.g. Hann) newwindow[0]=newwindow[Count]=0. xue@1: */ xue@1: void FillWindow(double* Result, WindowType wt, int Count, int* ips, double* dps) xue@1: { xue@1: switch(wt) xue@1: { xue@1: case wtRectangle: xue@1: for (int i=0; i<=Count; i++) xue@1: Result[i]=1; xue@1: break; xue@1: case wtTriangular: xue@1: for (int i=0; i<=Count; i++) xue@1: Result[i]=1-2*fabs(i-Count/2.0)/Count; xue@1: break; xue@1: case wtHamming: xue@1: for (int i=0; i<=Count; i++) xue@1: Result[i]=(0.54-0.46*cos(2*M_PI*i/Count)); xue@1: break; xue@1: case wtBlackman: xue@1: for (int i=0; i<=Count; i++) xue@1: Result[i]=0.42-0.5*cos(2*M_PI*i/Count)+0.08*cos(4*M_PI*i/Count); xue@1: break; xue@1: case wtGaussian: xue@1: { xue@1: double num=*dps*M_PI*M_PI/2/Count/Count; xue@1: for (int i=0; i<=Count; i++) xue@1: Result[i]=exp(-(i-Count/2.0)*(i-Count/2.0)*num); xue@1: break; xue@1: } xue@1: case wtKaiser: xue@1: { xue@1: double ldps=*dps*M_PI*M_PI/2/Count; xue@1: double den=I0(0.5*Count*ldps); xue@1: for (int i=0; i<=Count; i++) xue@1: Result[i]=I0(ldps*sqrt(0.25*Count*Count-(i-Count*0.5)*(i-Count*0.5)))/den; xue@1: break; xue@1: } xue@1: case wtHalfCosine: xue@1: { xue@1: for (int i=0; i<=Count; i++) xue@1: Result[i]=sin(M_PI*i/Count); xue@1: break; xue@1: } xue@1: case wtHann: xue@1: { xue@1: for (int i=0; i<=Count; i++) xue@1: Result[i]=(0.5-cos(M_PI*2*i/Count)/2); xue@1: break; xue@1: } xue@1: case wtHannSqr: xue@1: { xue@1: for (int i=0; i<=Count; i++) xue@1: Result[i]=(3-4*cos(M_PI*2*i/Count)+cos(M_PI*4*i/Count))/8; xue@1: break; xue@1: } xue@1: case wtHann3sqrt: xue@1: { xue@1: for (int i=0; i<=Count; i++) xue@1: { xue@1: double tmp=sin(M_PI*i/Count); xue@1: Result[i]=tmp*tmp*tmp; xue@1: } xue@1: break; xue@1: } xue@1: } xue@1: }//FillWindow xue@1: Chris@5: /** xue@1: function NewWindow: creates a window function. xue@1: xue@1: In: wt: window type xue@1: Count: window size xue@1: ips & ups: extra window specificatin parameters xue@1: Out: newwindow[Count+1]: the window function. xue@1: xue@1: Returns pointer to newwindow. newwindow is created anew if newwindow=0 is specified on start. xue@1: */ xue@1: double* NewWindow(WindowType wt, int Count, int* ips, double* dps, double* newwindow) xue@1: { xue@1: double* Result=newwindow; if (!Result) Result=new double[Count+1]; xue@1: FillWindow(Result, wt, Count, ips, dps); xue@1: return Result; xue@1: }//NewWindow xue@1: Chris@5: /** xue@1: function NewWindow8: 8-byte aligned version of NewWindow. xue@1: xue@1: In: wt: window type xue@1: Count: window size xue@1: ips & ups: extra window specificatin parameters xue@1: Out: newwindow[Count+1]: the window function. xue@1: xue@1: Returns pointer to newwindow. newwindow is created anew and 8-byte aligned if newwindow=0 is xue@1: specified on start. However, if newwindow is not 8-byte aligned on start, the unaligned buffer is xue@1: returned. xue@1: */ xue@1: double* NewWindow8(WindowType wt, int Count, int* ips, double* dps, double* newwindow) xue@1: { xue@1: double* Result=newwindow; if (!Result) Result=(double*)malloc8(sizeof(double)*(Count+1)); xue@1: FillWindow(Result, wt, Count, ips, dps); xue@1: return Result; xue@1: }//NewWindow8 xue@1: Chris@5: /** xue@1: function NewdWindow: computes the derivative of a window function. xue@1: xue@1: In: wt: window type xue@1: Count: window size xue@1: ips & ups: extra window specificatin parameters xue@1: Out: newdwindow[Count+1]: the derivative window function. xue@1: xue@1: Returns pointer to newdwindow. newdwindow is created anew if newdwindow=0 is specified on start. xue@1: */ xue@1: double* NewdWindow(WindowType wt, int Count, int* ips, double* dps, double* newdwindow) xue@1: { xue@1: double* Result=newdwindow; if (!Result) Result=new double[Count+1]; xue@1: double piiCount=M_PI/Count; xue@1: switch(wt) xue@1: { xue@1: case wtRectangle: xue@1: memset(Result, 0, sizeof(double)*(Count+1)); xue@1: break; xue@1: case wtTriangular: xue@1: throw("Trying to differentiate triangular window."); xue@1: break; xue@1: case wtHamming: xue@1: for (int i=0; i<=Count; i++) xue@1: Result[i]=0.92*piiCount*sin(2*piiCount*i); xue@1: break; xue@1: case wtBlackman: xue@1: for (int i=0; i<=Count; i++) xue@1: Result[i]=piiCount*sin(2*piiCount*i)-0.32*piiCount*sin(4*piiCount*i); xue@1: break; xue@1: case wtGaussian: xue@1: throw("Trying to differentiate Gaussian window."); xue@1: break; xue@1: case wtKaiser: xue@1: throw("Trying to differentiate Kaiser window."); xue@1: break; xue@1: case wtHalfCosine: xue@1: { xue@1: for (int i=0; i<=Count; i++) xue@1: Result[i]=piiCount*cos(piiCount*i); xue@1: break; xue@1: } xue@1: case wtHann: xue@1: { xue@1: for (int i=0; i<=Count; i++) xue@1: Result[i]=piiCount*sin(2*piiCount*i); xue@1: break; xue@1: } xue@1: case wtHannSqr: xue@1: { xue@1: for (int i=0; i<=Count; i++) xue@1: Result[i]=piiCount*sin(2*piiCount*i)-0.5*piiCount*sin(4*piiCount*i); xue@1: break; xue@1: } xue@1: case wtHann3sqrt: xue@1: { xue@1: for (int i=0; i<=Count; i++) xue@1: { xue@1: double s=sin(M_PI*i/Count), c=cos(M_PI*i/Count); xue@1: Result[i]=3*piiCount*s*s*c; xue@1: } xue@1: } xue@1: } xue@1: return Result; xue@1: }//NewdWindow xue@1: Chris@5: /** xue@1: function NewddWindow: computes the 2nd-order derivative of a window function. xue@1: xue@1: In: wt: window type xue@1: Count: window size xue@1: ips & ups: extra window specificatin parameters xue@1: Out: newddwindow[Count+1]: the 2nd-order derivative window function. xue@1: xue@1: Returns pointer to newddwindow. newddwindow is created anew if newddwindow=0 is specified on start. xue@1: */ xue@1: double* NewddWindow(WindowType wt, int Count, int* ips, double* dps, double* newddwindow) xue@1: { xue@1: double* Result=newddwindow; if (!Result) Result=new double[Count+1]; xue@1: double piiCount=M_PI/Count; xue@1: double piC2=piiCount*piiCount; xue@1: switch(wt) xue@1: { xue@1: case wtRectangle: xue@1: memset(Result, 0, sizeof(double)*(Count+1)); xue@1: break; xue@1: case wtTriangular: xue@1: throw("Trying to double-differentiate triangular window."); xue@1: break; xue@1: case wtHamming: xue@1: for (int i=0; i<=Count; i++) xue@1: Result[i]=1.84*piC2*cos(2*piiCount*i); xue@1: break; xue@1: case wtBlackman: xue@1: for (int i=0; i<=Count; i++) xue@1: Result[i]=2*piC2*cos(2*piiCount*i)-1.28*piC2*cos(4*piiCount*i); xue@1: break; xue@1: case wtGaussian: xue@1: throw("Trying to double-differentiate Gaussian window."); xue@1: break; xue@1: case wtKaiser: xue@1: throw("Trying to double-differentiate Kaiser window."); xue@1: break; xue@1: case wtHalfCosine: xue@1: { xue@1: for (int i=0; i<=Count; i++) xue@1: Result[i]=-piC2*sin(piiCount*i); xue@1: break; xue@1: } xue@1: case wtHann: xue@1: { xue@1: for (int i=0; i<=Count; i++) xue@1: Result[i]=2*piC2*cos(2*piiCount*i); xue@1: break; xue@1: } xue@1: case wtHannSqr: xue@1: { xue@1: for (int i=0; i<=Count; i++) xue@1: Result[i]=2*piC2*cos(2*piiCount*i)-2*piC2*cos(4*piiCount*i); xue@1: break; xue@1: } xue@1: case wtHann3sqrt: xue@1: { xue@1: for (int i=0; i<=Count; i++) xue@1: { xue@1: double s=sin(M_PI*i/Count), c=cos(M_PI*i/Count); xue@1: Result[i]=3*piC2*(2*c*c-s*s)*s; xue@1: } xue@1: break; xue@1: } xue@1: } xue@1: return Result; xue@1: }//NewddWindow xue@1: Chris@5: /** xue@1: function NewdddWindow: computes the 3rd-order derivative of a window function. xue@1: In: wt: window type xue@1: Count: window size xue@1: ips & ups: extra window specificatin parameters xue@1: Out: newdddwindow[Count+1]: the 3rd-order derivative window function. xue@1: xue@1: Returns pointer to newdddwindow. newdddwindow is created anew if newdddwindow=0 is specified on start. xue@1: */ xue@1: double* NewdddWindow(WindowType wt, int Count, int* ips, double* dps, double* newdddwindow) xue@1: { xue@1: double* Result=newdddwindow; if (!Result) Result=new double[Count+1]; xue@1: double piiCount=M_PI/Count; xue@1: double piC2=piiCount*piiCount; xue@1: double piC3=piiCount*piC2; xue@1: switch(wt) xue@1: { xue@1: case wtRectangle: xue@1: memset(Result, 0, sizeof(double)*(Count+1)); xue@1: break; xue@1: case wtTriangular: xue@1: throw("Trying to triple-differentiate triangular window."); xue@1: break; xue@1: case wtHamming: xue@1: for (int i=0; i<=Count; i++) xue@1: Result[i]=-3.68*piC3*sin(2*piiCount*i); xue@1: break; xue@1: case wtBlackman: xue@1: for (int i=0; i<=Count; i++) xue@1: Result[i]=-4*piC3*sin(2*piiCount*i)+5.12*piC3*sin(4*piiCount*i); xue@1: break; xue@1: case wtGaussian: xue@1: throw("Trying to triple-differentiate Gaussian window."); xue@1: break; xue@1: case wtKaiser: xue@1: throw("Trying to triple-differentiate Kaiser window."); xue@1: break; xue@1: case wtHalfCosine: xue@1: { xue@1: for (int i=0; i<=Count; i++) xue@1: Result[i]=-piC3*cos(piiCount*i); xue@1: break; xue@1: } xue@1: case wtHann: xue@1: { xue@1: for (int i=0; i<=Count; i++) xue@1: Result[i]=-4*piC3*sin(2*piiCount*i); xue@1: break; xue@1: } xue@1: case wtHannSqr: xue@1: { xue@1: for (int i=0; i<=Count; i++) xue@1: Result[i]=-4*piC3*sin(2*piiCount*i)+8*piC3*sin(4*piiCount*i); xue@1: break; xue@1: } xue@1: case wtHann3sqrt: xue@1: throw("Trying to triple-differentiate Hann^1.5 window."); xue@1: break; xue@1: } xue@1: return Result; xue@1: }//NewdddWindow xue@1: xue@1: //--------------------------------------------------------------------------- Chris@5: /** xue@1: function windowspec: computes a few descriptors of a cosine family window. A window function in the xue@1: cosine window family is the linear combination of a few cosine functions, therefore has a cosine xue@1: decomposition. xue@1: xue@1: In: wt: window type xue@1: N: window size xue@1: Out: c[M+1], coefficients of cosine decomposition xue@1: iH2: reciprocal of square of L2 norm, xue@1: d[2M+1]: self-convolution of c[], optional xue@1: xue@1: No return value. xue@1: */ xue@1: void windowspec(WindowType wt, int N, int* M, double* c, double* iH2, double* d) xue@1: { xue@1: switch(wt) xue@1: { xue@1: case wtRectangle: xue@1: M[0]=0, c[0]=1, iH2[0]=1.0/N/N; xue@1: if (d) d[0]=1; xue@1: break; xue@1: case wtHamming: xue@1: M[0]=1, c[0]=0.54, c[1]=0.23, iH2[0]=1.0/(0.3974*N*N); xue@1: if (d) d[0]=0.3974, d[1]=0.2484, d[2]=0.0529; xue@1: break; xue@1: case wtBlackman: xue@1: M[0]=2, c[0]=0.42, c[1]=0.25, c[2]=0.04, iH2[0]=1.0/(0.3046*N*N); xue@1: if (d) d[0]=0.3046, d[1]=0.23, d[2]=0.0961, d[4]=0.02, d[5]=0.0016; xue@1: break; xue@1: case wtHann: xue@1: M[0]=1, c[0]=0.5, c[1]=0.25, iH2[0]=1.0/(0.375*N*N); xue@1: if (d) d[0]=0.375, d[1]=0.25, d[2]=0.0625; xue@1: break; xue@1: default: xue@1: M[0]=1, c[0]=0.5, c[1]=0.25, iH2[0]=1.0/(0.375*N*N); xue@1: if (d) d[0]=0.375, d[1]=0.25, d[2]=0.0625; xue@1: break; xue@1: } xue@1: }//windowspec