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: 
xue@1: #include "align8.h"
Chris@2: #include "sinsyn.h"
xue@1: #include "splines.h"
xue@1: 
Chris@5: /** \file sinsyn.h */
Chris@5: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@7:   function CosSin: recursive cos-sin generator with trinomial frequency
xue@7: 
xue@7:   In: CountSt, CountEn
xue@7:       f3, f2, f1, f0: trinomial coefficients of frequency
xue@7:       ph: initial phase angle at 0 (NOT at CountSt)
xue@7:   Out: datar[CountSt:CountEn-1], datai[CountSt:CountEn-1]: synthesized pair of cosine and sine functions
xue@7:        ph: phase angle at CountEn
xue@7: 
xue@7:   No return value.
xue@7: */
xue@7: void CosSin(double* datar, double* datai, int CountSt, int CountEn, double f3, double f2, double f1, double f0, double &ph)
xue@7: {
xue@7:   int i;
xue@7:   double dph, ddph, dddph, ddddph,
xue@7:          sph, cph, sdph, cdph, sddph, cddph, sdddph, cdddph, sddddph, cddddph,
xue@7:          p0=ph, p1=2*M_PI*f0, p2=M_PI*f1, p3=2.0*M_PI*f2/3, p4=2.0*M_PI*f3/4, tmp;
xue@7:   if (CountSt==0)
xue@7:   {
xue@7:     dph=p1+p2+p3+p4; ddph=2*p2+6*p3+14*p4; dddph=6*p3+36*p4; ddddph=24*p4;
xue@7:   }
xue@7:   else
xue@7:   {
xue@7:     ph=p0+CountSt*(p1+CountSt*(p2+CountSt*(p3+CountSt*p4)));
xue@7:     dph=p1+p2+p3+p4+CountSt*(2*p2+3*p3+4*p4+CountSt*(3*p3+6*p4+CountSt*4*p4));
xue@7:     ddph=2*p2+6*p3+14*p4+CountSt*(6*p3+24*p4+CountSt*12*p4);
xue@7:     dddph=6*p3+36*p4+CountSt*24*p4; ddddph=24*p4;
xue@7:   }
xue@7:   sph=sin(ph), cph=cos(ph);
xue@7:   sdph=sin(dph), cdph=cos(dph);
xue@7:   sddph=sin(ddph), cddph=cos(ddph);
xue@7:   sdddph=sin(dddph), cdddph=cos(dddph);
xue@7:   sddddph=sin(ddddph), cddddph=cos(ddddph);
xue@7: 
xue@7:   for (i=CountSt; i<CountEn; i++)
xue@7:   {
xue@7:     datar[i]=cph; datai[i]=sph;
xue@7:     tmp=cph*cdph-sph*sdph; sph=sph*cdph+cph*sdph; cph=tmp;
xue@7:     tmp=cdph*cddph-sdph*sddph; sdph=sdph*cddph+cdph*sddph; cdph=tmp;
xue@7:     tmp=cddph*cdddph-sddph*sdddph; sddph=sddph*cdddph+cddph*sdddph; cddph=tmp;
xue@7:     tmp=cdddph*cddddph-sdddph*sddddph; sdddph=sdddph*cddddph+cdddph*sddddph; cdddph=tmp;
xue@7:   }
xue@7:   ph=p0+CountEn*(p1+CountEn*(p2+CountEn*(p3+CountEn*p4)));
xue@7: }//CosSin*/
xue@7: 
xue@7: /**
xue@1:   function Sinuoid:	original McAuley-Quatieri synthesizer interpolation between two measurement points.
xue@1: 
xue@1:   In: T: length from measurement point 1 to measurement point 2
xue@1:       a1, f1, p2: amplitude, frequency and phase angle at measurement point 1
xue@1:       a2, f2, p2: amplitude, frequency and phase angle at measurement point 2
xue@1:       ad: specifies if the resynthesized sinusoid is to be added to or to replace the contents of output buffer
xue@1:   Out: data[T]: output buffer
xue@1:        a[T], f[T], p[T]: resynthesized amplitude, frequency and phase
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void Sinusoid(double* data, int T, double a1, double a2, double f1, double f2, double p1, double p2, double* a, double* f, double* p, bool ad)
xue@1: {
xue@1:   int M=floor(((p1-p2)/M_PI+(f1+f2)*T)/2.0+0.5);
xue@1:   double b1=p2-p1-2*M_PI*(f1*T-M), b2=2*M_PI*(f2-f1);
xue@1:   double pa=(3*b1/T-b2)/T, pb=(-2*b1/T+b2)/T/T, pc=2*M_PI*f1, pd=p1;
xue@1:   double la=a1, da=(a2-a1)/T;
xue@1:   if (ad)
xue@1:     for (int t=0; t<T; t++)
xue@1:     {
xue@1:       double lp=pd+t*(pc+t*(pa+t*pb)), lf=(pc+2*pa*t+3*pb*t*t)/2/M_PI;
xue@1:       data[t]+=la*cos(lp);
xue@1:       a[t]=la;
xue@1:       p[t]=lp;
xue@1:       f[t]=lf;
xue@1:       la=la+da;
xue@1:     }
xue@1:   else
xue@1:     for (int t=0; t<T; t++)
xue@1:     {
xue@1:       double lp=pd+t*(pc+t*(pa+t*pb)), lf=(pc+2*pa*t+3*pb*t*t)/2/M_PI;
xue@1:       data[t]=la*cos(lp);
xue@1:       a[t]=la;
xue@1:       p[t]=lp;
xue@1:       f[t]=lf;
xue@1:       la=la+da;
xue@1:     }
xue@1: }//Sinusoid
xue@1: 
Chris@5: /**
xue@1:   function Sinuoid:	original McAuley-Quatieri synthesizer interpolation between two measurement points,
xue@1:   without returning interpolated sinusoid parameters.
xue@1: 
xue@1:   In: T: length from measurement point 1 to measurement point 2
xue@1:       a1, f1, p2: amplitude, frequency and phase angle at measurement point 1
xue@1:       a2, f2, p2: amplitude, frequency and phase angle at measurement point 2
xue@1:       ad: specifies if the resynthesized sinusoid is to be added to or to replace the contents of output buffer
xue@1:   Out: data[T]: output buffer
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void Sinusoid(double* data, int T, double a1, double a2, double f1, double f2, double p1, double p2, bool ad)
xue@1: {
xue@1:   int M=floor(((p1-p2)/M_PI+(f1+f2)*T)/2.0+0.5);
xue@1:   double b1=p2-p1-2*M_PI*(f1*T-M), b2=2*M_PI*(f2-f1);
xue@1:   double pa=(3*b1/T-b2)/T, pb=(-2*b1/T+b2)/T/T, pc=2*M_PI*f1, pd=p1;
xue@1:   double la=a1, da=(a2-a1)/T;
xue@1:   if (ad)
xue@1:     for (int t=0; t<T; t++)
xue@1:     {
xue@1:       data[t]+=la*cos(pd+t*(pc+t*(pa+t*pb)));
xue@1:       la=la+da;
xue@1:     }
xue@1:   else
xue@1:     for (int t=0; t<T; t++)
xue@1:     {
xue@1:       data[t]=la*cos(pd+t*(pc+t*(pa+t*pb)));
xue@1:       la=la+da;
xue@1:     }
xue@1: }//Sinusoid
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function Sinusoid_direct: synthesizes sinusoid over [CountSt, CountEn) from tronomial coefficients of
xue@1:   amplitude and frequency, direct implementation.
xue@1: 
xue@1:   In: CountSt, CountEn
xue@1:       aa, ab, ac, ad: trinomial coefficients of amplitude
xue@1:       fa, fb, fc, fd: trinomial coefficients of frequency
xue@1:       p1: initial phase angle at 0 (NOT at CountSt)
xue@1:       add: specifies if the resynthesized sinusoid is to be added to or to replace the content of output buffer
xue@1:   Out: data[CountSt:CountEn-1]: output buffer.
xue@1:        p1: phase angle at CountEn
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void Sinusoid_direct(double* data, int CountSt, int CountEn, double aa, double ab, double ac, double ad,
xue@1:   double fa, double fb, double fc, double fd, double &p1, bool add)
xue@1: {
xue@1:   int i; double a, ph;
xue@1:   for (i=CountSt; i<CountEn; i++)
xue@1:   {
xue@1:      a=ad+i*(ac+i*(ab+i*aa));
xue@1:      ph=p1+2*M_PI*i*(fd+i*((fc/2)+i*((fb/3)+i*fa/4)));
xue@1:     if (add) data[i]+=a*cos(ph);
xue@1:     else data[i]=a*cos(ph);
xue@1:   }
xue@1:   p1=p1+2*M_PI*i*(fd+i*((fc/2)+i*((fb/3)+i*fa/4)));
xue@1: }//Sinusoid
xue@1: 
Chris@5: /**
xue@7:   function Sinusoid_direct: synthesizes sinusoid over [CountSt, CountEn) from tronomial coefficients of
xue@7:   amplitude and frequency, direct implementation returning sinusoid coefficients rather than waveform.
xue@7: 
xue@7:   In: CountSt, CountEn
xue@7:       aa, ab, ac, ad: trinomial coefficients of amplitude
xue@7:       fa, fb, fc, fd: trinomial coefficients of frequency
xue@7:       p1: initial phase angle at 0 (NOT at CountSt)
xue@7:   Out: f[CountSt:CountEn-1], a[:], ph[:], da[:]: output buffers.
xue@7:        p1: phase angle at CountEn
xue@7: 
xue@7:   No return value.
xue@7: */
xue@7: void Sinusoid_direct(double* f, double* a, double* ph, double* da, int CountSt, int CountEn, double aa, double ab, double ac, double ad,
xue@7:   double fa, double fb, double fc, double fd, double &p1, bool LogA)
xue@7: {
xue@7:   int i;
xue@7:   if (LogA)
xue@7:   {
xue@7:     for (i=CountSt; i<CountEn; i++)
xue@7:     {
xue@7:        a[i]=exp(ad+i*(ac+i*(ab+i*aa)));
xue@7:        f[i]=fd+i*(fc+i*(fb+i*fa));
xue@7:        ph[i]=p1+2*M_PI*i*(fd+i*((fc/2)+i*((fb/3)+i*fa/4)));
xue@7:        da[i]=a[i]*(ac+i*(2*ab+3*i*aa));
xue@7:     }
xue@7:   }
xue@7:   else
xue@7:   {
xue@7:     for (i=CountSt; i<CountEn; i++)
xue@7:     {
xue@7:        a[i]=ad+i*(ac+i*(ab+i*aa));
xue@7:        f[i]=fd+i*(fc+i*(fb+i*fa));
xue@7:        ph[i]=p1+2*M_PI*i*(fd+i*((fc/2)+i*((fb/3)+i*fa/4)));
xue@7:        da[i]=ac+i*(2*ab+3*i*aa);
xue@7:     }
xue@7:   }
xue@7:   p1=p1+2*M_PI*i*(fd+i*((fc/2)+i*((fb/3)+i*fa/4)));
xue@7: }//Sinusoid
xue@7: 
xue@7: /**
xue@7:   function Sinusoid_direct: synthesizes sinusoid over [CountSt, CountEn) from tronomial coefficients of
xue@7:   frequency, direct implementation returning frequency and phase rather than waveform.
xue@7: 
xue@7:   In: CountSt, CountEn
xue@7:       fa, fb, fc, fd: trinomial coefficients of frequency
xue@7:       p1: initial phase angle at 0 (NOT at CountSt)
xue@7:   Out: f[CountSt:CountEn-1], ph[CountSt:CountEn-1]: output buffers.
xue@7:        p1: phase angle at CountEn
xue@7: 
xue@7:   No return value.
xue@7: */
xue@7: void Sinusoid_direct(double* f, double* ph, int CountSt, int CountEn, double fa, double fb, double fc,
xue@7:       double fd, double &p1)
xue@7: {
xue@7:   int i;
xue@7:   for (i=CountSt; i<CountEn; i++)
xue@7:   {
xue@7:      f[i]=fd+i*(fc+i*(fb+i*fa));
xue@7:      ph[i]=p1+2*M_PI*i*(fd+i*((fc/2)+i*((fb/3)+i*fa/4)));
xue@7:   }
xue@7:   p1=p1+2*M_PI*i*(fd+i*((fc/2)+i*((fb/3)+i*fa/4)));
xue@7: }//Sinusoid
xue@7: 
xue@7: /**
xue@1:   function Sinusoid: synthesizes sinusoid over [CountSt, CountEn) from tronomial coefficients of
xue@1:   amplitude and frequency, recursive implementation.
xue@1: 
xue@1:   In: CountSt, CountEn
xue@1:       a3, a2, a1, a0: trinomial coefficients of amplitude
xue@1:       f3, f2, f1, f0: trinomial coefficients of frequency
xue@1:       ph: initial phase angle at 0 (NOT at CountSt)
xue@1:       add: specifies if the resynthesized sinusoid is to be added to or to replace the content of output buffer
xue@1:   Out: data[CountSt:CountEn-1]: output buffer.
xue@1:        ph: phase angle at CountEn
xue@1: 
xue@1:   No return value. This function requires 8-byte stack alignment for optimal speed.
xue@1: */
xue@1: void Sinusoid(double* data, int CountSt, int CountEn, double a3, double a2, double a1, double a0,
xue@1:   double f3, double f2, double f1, double f0, double &ph, bool add)
xue@1: {
xue@1:   int i;
xue@1:   double a, da, dda, ddda, dph, ddph, dddph, ddddph,
xue@1:          sph, cph, sdph, cdph, sddph, cddph, sdddph, cdddph, sddddph, cddddph,
xue@1:          p0=ph, p1=2*M_PI*f0, p2=M_PI*f1, p3=2.0*M_PI*f2/3, p4=2.0*M_PI*f3/4, tmp;
xue@1:   if (CountSt==0)
xue@1:   {
xue@1:     a=a0; da=a1+a2+a3; dda=2*a2+6*a3; ddda=6*a3;
xue@1:     dph=p1+p2+p3+p4; ddph=2*p2+6*p3+14*p4; dddph=6*p3+36*p4; ddddph=24*p4;
xue@1:   }
xue@1:   else
xue@1:   {
xue@1:     a=a0+CountSt*(a1+CountSt*(a2+CountSt*a3));
xue@1:     da=a1+a2+a3+CountSt*(2*a2+3*a3+CountSt*3*a3);
xue@1:     dda=2*a2+6*a3+CountSt*6*a3; ddda=6*a3;
xue@1:     ph=p0+CountSt*(p1+CountSt*(p2+CountSt*(p3+CountSt*p4)));
xue@1:     dph=p1+p2+p3+p4+CountSt*(2*p2+3*p3+4*p4+CountSt*(3*p3+6*p4+CountSt*4*p4));
xue@1:     ddph=2*p2+6*p3+14*p4+CountSt*(6*p3+24*p4+CountSt*12*p4);
xue@1:     dddph=6*p3+36*p4+CountSt*24*p4; ddddph=24*p4;
xue@1:   }
xue@1:   sph=sin(ph), cph=cos(ph);
xue@1:   sdph=sin(dph), cdph=cos(dph);
xue@1:   sddph=sin(ddph), cddph=cos(ddph);
xue@1:   sdddph=sin(dddph), cdddph=cos(dddph);
xue@1:   sddddph=sin(ddddph), cddddph=cos(ddddph);
xue@1:   if (add)
xue@1:   {
xue@1:     for (i=CountSt; i<CountEn; i++)
xue@1:     {
xue@1:       data[i]+=a*cph;
xue@1:       a=a+da; da=da+dda; dda=dda+ddda;
xue@1:       tmp=cph*cdph-sph*sdph; sph=sph*cdph+cph*sdph; cph=tmp;
xue@1:       tmp=cdph*cddph-sdph*sddph; sdph=sdph*cddph+cdph*sddph; cdph=tmp;
xue@1:       tmp=cddph*cdddph-sddph*sdddph; sddph=sddph*cdddph+cddph*sdddph; cddph=tmp;
xue@1:       tmp=cdddph*cddddph-sdddph*sddddph; sdddph=sdddph*cddddph+cdddph*sddddph; cdddph=tmp;
xue@1:     }
xue@1:   }
xue@1:   else
xue@1:   {
xue@1:     for (i=CountSt; i<CountEn; i++)
xue@1:     {
xue@1:       data[i]=a*cph;
xue@1:       a=a+da; da=da+dda; dda=dda+ddda;
xue@1:       tmp=cph*cdph-sph*sdph; sph=sph*cdph+cph*sdph; cph=tmp;
xue@1:       tmp=cdph*cddph-sdph*sddph; sdph=sdph*cddph+cdph*sddph; cdph=tmp;
xue@1:       tmp=cddph*cdddph-sddph*sdddph; sddph=sddph*cdddph+cddph*sdddph; cddph=tmp;
xue@1:       tmp=cdddph*cddddph-sdddph*sddddph; sdddph=sdddph*cddddph+cdddph*sddddph; cdddph=tmp;
xue@1:     }
xue@1:   }
xue@1:   ph=p0+CountEn*(p1+CountEn*(p2+CountEn*(p3+CountEn*p4)));
xue@1: }
xue@1: 
Chris@5: /**
xue@1:   function SinusoidExp: synthesizes complex sinusoid whose derivative log amplitude and frequency are
xue@1:   trinomials
xue@1: 
xue@1:   In: CountSt, CountEn
xue@1:       a3, a2, a1, a0: trinomial coefficients for the derivative of log amplitude
xue@1:       omg3, omg2, omg1, omg0: trinomial coefficients for angular frequency
xue@1:       ea, ph: initial log amplitude and phase angle at 0
xue@1:       add: specifies if the resynthesized sinusoid is to be added to or to replace the content of output buffer
xue@1:   Out: data[CountSt:CountEn-1]: output buffer.
xue@1:        ea, ph: log amplitude and phase angle at CountEn.
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void SinusoidExp(cdouble* data, int CountSt, int CountEn, double a3, double a2, double a1, double a0,
xue@1: 	double omg3, double omg2, double omg1, double omg0, double &ea, double &ph, bool add)
xue@1: {
xue@1: 	double e0=ea, e1=a0, e2=0.5*a1, e3=a2/3, e4=a3/4,
xue@1:          p0=ph, p1=omg0, p2=0.5*omg1, p3=omg2/3, p4=omg3/4;
xue@1: 	if (add) for (int i=CountSt; i<CountEn; i++)
xue@1: 	{
xue@1: 		double lea=e0+i*(e1+i*(e2+i*(e3+i*e4)));
xue@1: 		double lph=p0+i*(p1+i*(p2+i*(p3+i*p4)));
xue@1: 		data[i]+=exp(cdouble(lea, lph));
xue@1: 	}
xue@1: 	else for (int i=CountSt; i<CountEn; i++)
xue@1: 	{
xue@1: 		double lea=e0+i*(e1+i*(e2+i*(e3+i*e4)));
xue@1: 		double lph=p0+i*(p1+i*(p2+i*(p3+i*p4)));
xue@1: 		data[i]=exp(cdouble(lea, lph));
xue@1: 	}
xue@1: 	ea=e0+CountEn*(e1+CountEn*(e2+CountEn*(e3+CountEn*e4)));
xue@1: 	ph=p0+CountEn*(p1+CountEn*(p2+CountEn*(p3+CountEn*p4)));
xue@1: }//SinusoidExp
xue@1: 
Chris@5: /**
xue@1:   function SinusoidExp: synthesizes complex sinusoid piece whose derivative logarithm is h[M]'lamda[M].
xue@1:   This version also synthesizes its derivative.
xue@1: 
xue@1:   In: h[M][T], dih[M][T]: basis functions and their difference-integrals
xue@1:       lamda[M]: coefficients of h[M]
xue@1:       tmpexp: inital logarithm at 0
xue@1:   Out: s[T], ds[T]: synthesized sinusoid and its derivative
xue@1:        tmpexp: logarithm at T
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void SinusoidExp(int T, cdouble* s, cdouble* ds, int M, cdouble* lamda, double** h, double** dih, cdouble& tmpexp)
xue@1: {
xue@1: 	for (int t=0; t<T; t++)
xue@1: 	{
xue@1: 		s[t]=exp(tmpexp);
xue@1: 		cdouble dexp=0, dR=0;
xue@1: 		for (int m=0; m<M; m++)	dexp+=lamda[m]*dih[m][t], dR+=lamda[m]*h[m][t];
xue@1: 		tmpexp+=dexp;
xue@1: 		ds[t]=s[t]*dR;
xue@1: 	}
xue@1: }//SinusoidExp
xue@1: 
Chris@5: /**
xue@1:   function SinusoidExp: synthesizes complex sinusoid piece whose derivative logarithm is h[M]'lamda[M].
xue@1:   This version does not synthesize its derivative.
xue@1: 
xue@1:   In: dih[M][T]: difference of integrals of basis functions h[M]
xue@1:       lamda[M]: coefficients of h[M]
xue@1:       tmpexp: inital logarithm at 0
xue@1:   Out: s[T]: synthesized sinusoid
xue@1:        tmpexp: logarithm at T
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void SinusoidExp(int T, cdouble* s, int M, cdouble* lamda, double** dih, cdouble& tmpexp)
xue@1: {
xue@1: 	for (int t=0; t<T; t++)
xue@1: 	{
xue@1: 		s[t]=exp(tmpexp);
xue@1: 		cdouble dexp=0;
xue@1: 		for (int m=0; m<M; m++)	dexp+=lamda[m]*dih[m][t];
xue@1: 		tmpexp+=dexp;
xue@1: 	}
xue@1: }//SinusoidExp
xue@1: 
Chris@5: /**
xue@1:   function SinusoidExpA: synthesizes complex sinusoid whose log amplitude and frequency are trinomials
xue@1: 
xue@1:   In: CountSt, CountEn
xue@1:       a3, a2, a1, a0: trinomial coefficients for log amplitude
xue@1:       omg3, omg2, omg1, omg0: trinomial coefficients for angular frequency
xue@1:       ph: initial phase angle at 0
xue@1:       add: specifies if the resynthesized sinusoid is to be added to or to replace the content of output buffer
xue@1:   Out: data[CountSt:CountEn-1]: output buffer.
xue@1:        ph: phase angle at CountEn.
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void SinusoidExpA(cdouble* data, int CountSt, int CountEn, double a3, double a2, double a1, double a0,
xue@1: 	double omg3, double omg2, double omg1, double omg0, double &ph, bool add)
xue@1: {
xue@1:   double p0=ph, p1=omg0, p2=0.5*omg1, p3=omg2/3, p4=omg3/4;
xue@1: 	if (add) for (int i=CountSt; i<CountEn; i++)
xue@1: 	{
xue@1: 		double lea=a0+i*(a1+i*(a2+i*a3));
xue@1: 		double lph=p0+i*(p1+i*(p2+i*(p3+i*p4)));
xue@1: 		data[i]+=exp(cdouble(lea, lph));
xue@1: 	}
xue@1: 	else for (int i=CountSt; i<CountEn; i++)
xue@1: 	{
xue@1: 		double lea=a0+i*(a1+i*(a2+i*a3));
xue@1: 		double lph=p0+i*(p1+i*(p2+i*(p3+i*p4)));
xue@1: 		data[i]=exp(cdouble(lea, lph));
xue@1: 	}
xue@1: 	ph=p0+CountEn*(p1+CountEn*(p2+CountEn*(p3+CountEn*p4)));
xue@1: }//SinusoidExpA
xue@1: 
Chris@5: /**
xue@1:   function SinusoidExpA: synthesizes complex sinusoid whose log amplitude and frequency are trinomials
xue@1:   with phase angle specified at both ends.
xue@1: 
xue@7:   In: T
xue@1:       a3, a2, a1, a0: trinomial coefficients for log amplitude
xue@1:       omg3, omg2, omg1, omg0: trinomial coefficients for angular frequency
xue@1:       ph0, ph2: phase angles at 0 and CountEn.
xue@1:       add: specifies if the resynthesized sinusoid is to be added to or to replace the content of output buffer
xue@7:   Out: data[T]: output buffer.
xue@1: 
xue@1:   No return value.
xue@1: */
xue@7: void SinusoidExpA(int T, cdouble* data, double a3, double a2, double a1, double a0,
xue@1: 	double omg3, double omg2, double omg1, double omg0, double ph0, double ph2, bool add)
xue@1: {
xue@1:   double p0=ph0, p1=omg0, p2=0.5*omg1, p3=omg2/3, p4=omg3/4;
xue@7:   double pend=p0+T*(p1+T*(p2+T*(p3+T*p4)));
xue@1: 
xue@1: 	int k=floor((pend-ph2)/2/M_PI+0.5);
xue@1: 	double d=ph2-pend+2*M_PI*k;
xue@7:   double _p=-2*d/T/T/T;
xue@7:   double _q=3*d/T/T;
xue@1: 
xue@7:   if (add) for (int i=0; i<T; i++)
xue@1: 	{
xue@1: 		double lea=a0+i*(a1+i*(a2+i*a3));
xue@1: 		double lph=p0+i*(p1+i*(p2+i*(p3+i*p4)));
xue@1:     data[i]+=exp(cdouble(lea, lph+(i*i*(_q+i*_p))));
xue@1: 	}
xue@7:   else for (int i=0; i<T; i++)
xue@1: 	{
xue@1: 		double lea=a0+i*(a1+i*(a2+i*a3));
xue@1: 		double lph=p0+i*(p1+i*(p2+i*(p3+i*p4)));
xue@1: 		data[i]=exp(cdouble(lea, lph+(i*i*(_q+i*_p))));
xue@1: 	}
xue@1: }//SinusoidExpA
xue@1: 
Chris@5: /**
xue@1:   function SinusoidExpA: synthesizes complex sinusoid piece whose log amplitude is h[M]'p[M] and
xue@1:   frequency is h[M]'q[M]. This version also synthesizes its derivative.
xue@1: 
xue@1:   In: h[M][T], dh[M][T], dih[M][T]: basis functions and their derivatives and difference-integrals
xue@1:       p[M], q[M]: real and imaginary parts of coefficients of h[M]
xue@1:       tmpph: inital phase angle at 0
xue@1:   Out: s[T], ds[T]: synthesized sinusoid and its derivative
xue@1:        tmpph: phase angle at T
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void SinusoidExpA(int T, cdouble* s, cdouble* ds, int M, double* p, double* q, double** h, double** dh, double** dih, double& tmpph)
xue@1: {
xue@1: 	for (int t=0; t<T; t++)
xue@1: 	{
xue@1: 		double e=0, dph=0, drr=0, dri=0;
xue@1: 		for (int m=0; m<M; m++) e+=p[m]*h[m][t], dph+=q[m]*dih[m][t], drr+=p[m]*dh[m][t], dri+=q[m]*h[m][t];
xue@1: 		s[t]=exp(cdouble(e, tmpph));
xue@1: 		ds[t]=s[t]*cdouble(drr, dri);
xue@1: 		tmpph+=dph;
xue@1: 	}
xue@1: }//SinusoidExpA
xue@1: 
Chris@5: /**
xue@1:   function SinusoidExpA: synthesizes complex sinusoid piece whose log amplitude is h[M]'p[M] and
xue@1:   frequency is h[M]'q[M]. This version does not synthesize its derivative.
xue@1: 
xue@1:   In: h[M][T], dih[M][T]: basis functions and their difference-integrals
xue@1:       p[M], q[M]: real and imaginary parts of coefficients of h[M]
xue@1:       tmpph: inital phase angle at 0
xue@1:   Out: s[T]: synthesized sinusoid
xue@1:        tmpph: phase angle at T
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void SinusoidExpA(int T, cdouble* s, int M, double* p, double* q, double** h, double** dih, double& tmpph)
xue@1: {
xue@1: 	for (int t=0; t<T; t++)
xue@1: 	{
xue@1: 		double e=0, dph=0;
xue@1: 		for (int m=0; m<M; m++) e+=p[m]*h[m][t], dph+=q[m]*dih[m][t];
xue@1: 		s[t]=exp(cdouble(e, tmpph));
xue@1: 		tmpph+=dph;
xue@1: 	}
xue@1: }//SinusoidExpA
xue@1: 
Chris@5: /**
xue@1:   function SinusoidExpA: synthesizes complex sinusoid piece whose log amplitude is h[M]'p[M] and
xue@1:   frequency is h[M]'q[M] with phase angle specified at both ends. This version does not synthesize its
xue@1:   derivative.
xue@1: 
xue@1:   In: h[M][T], dih[M][T]: basis functions and their difference-integrals
xue@1:       p[M], q[M]: real and imaginary parts of coefficients of h[M]
xue@1:       ph1, ph2: phase angles at 0 and T.
xue@1:   Out: s[T]: synthesized sinusoid
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void SinusoidExpA(int T, cdouble* s, int M, double* p, double* q, double** h, double** dih, double ph1, double ph2)
xue@1: {
xue@1: 	double pend=ph1;
xue@1: 	for (int t=0; t<T; t++)
xue@1: 	{
xue@1: 		double dph=0;
xue@1: 		for (int m=0; m<M; m++) dph+=q[m]*dih[m][t];
xue@1: 		pend+=dph;
xue@1: 	}
xue@1: 
xue@1: 	int k=floor((pend-ph2)/2/M_PI+0.5);
xue@1: 	double d=ph2-pend+2*M_PI*k;
xue@1: 	double _p=-2*d/T/T/T;
xue@1: 	double _q=3*d/T/T;
xue@1: 
xue@1: 	double ph=ph1;
xue@1: 	for (int t=0; t<T; t++)
xue@1: 	{
xue@1: 		double e=0, dph=0;
xue@1: 		for (int m=0; m<M; m++) e+=p[m]*h[m][t], dph+=q[m]*dih[m][t];
xue@1: 		if (e>300) e=300;
xue@1: 		if (e<-300) e=-300;
xue@1: 		s[t]=exp(cdouble(e, ph+(t*t*(_q+t*_p))));
xue@1: 		ph+=dph;
xue@1: 	}
xue@1: }//SinusoidExpA
xue@1: 
xue@1: /*
xue@1: //This is not used any longer as the recursion does not seem to help saving computation with all its overheads.
xue@1: void SinusoidExp(cdouble* data, int CountSt, int CountEn, double a3, double a2, double a1, double a0,
xue@1: 	double omg3, double omg2, double omg1, double omg0, double &ea, double &ph, bool add)
xue@1: {
xue@1:   int i;
xue@1: 	double dea, ddea, dddea, ddddea,
xue@1: 				 dph, ddph, dddph, ddddph,
xue@1: 				 sph, cph, sdph, cdph, sddph, cddph, sdddph, cdddph, sddddph, cddddph,
xue@1: 				 e0=ea, e1=a0, e2=0.5*a1, e3=a2/3, e4=a3/4,
xue@1: 				 p0=ph, p1=omg0, p2=0.5*omg1, p3=omg2/3, p4=omg3/4, tmp;
xue@1: 	if (CountSt==0)
xue@1:   {
xue@1: 		dea=e1+e2+e3+e4; ddea=2*e2+6*e3+14*e4; dddea=6*e3+36*e4; ddddea=24*e3;
xue@1: 		dph=p1+p2+p3+p4; ddph=2*p2+6*p3+14*p4; dddph=6*p3+36*p4; ddddph=24*p4;
xue@1:   }
xue@1:   else
xue@1:   {
xue@1: 		ea=e0+CountSt*(e1+CountSt*(e2+CountSt*(e3+CountSt*e4)));
xue@1: 		dea=e1+e2+e3+e4+CountSt*(2*e2+3*e3+4*e4+CountSt*(3*e3+6*e4+CountSt*4*e4));
xue@1: 		ddea=2*e2+6*e3+14*e4+CountSt*(6*e3+24*e4+CountSt*12*e4);
xue@1: 		dddea=6*e3+36*e4+CountSt*24*e4; ddddea=24*e4;
xue@1: 		ph=p0+CountSt*(p1+CountSt*(p2+CountSt*(p3+CountSt*p4)));
xue@1: 		dph=p1+p2+p3+p4+CountSt*(2*p2+3*p3+4*p4+CountSt*(3*p3+6*p4+CountSt*4*p4));
xue@1: 		ddph=2*p2+6*p3+14*p4+CountSt*(6*p3+24*p4+CountSt*12*p4);
xue@1: 		dddph=6*p3+36*p4+CountSt*24*p4; ddddph=24*p4;
xue@1: 	}
xue@1:   sph=sin(ph), cph=cos(ph);
xue@1:   sdph=sin(dph), cdph=cos(dph);
xue@1:   sddph=sin(ddph), cddph=cos(ddph);
xue@1:   sdddph=sin(dddph), cdddph=cos(dddph);
xue@1: 	sddddph=sin(ddddph), cddddph=cos(ddddph);
xue@1:   if (add)
xue@1:   {
xue@1:     for (i=CountSt; i<CountEn; i++)
xue@1:     {
xue@1: 			data[i]+=exp(ea)*cdouble(cph, sph);
xue@1: 			ea=ea+dea; dea=dea+ddea; ddea=ddea+dddea; dddea+dddea+ddddea;
xue@1: 			tmp=cph*cdph-sph*sdph; sph=sph*cdph+cph*sdph; cph=tmp;
xue@1: 			tmp=cdph*cddph-sdph*sddph; sdph=sdph*cddph+cdph*sddph; cdph=tmp;
xue@1: 			tmp=cddph*cdddph-sddph*sdddph; sddph=sddph*cdddph+cddph*sdddph; cddph=tmp;
xue@1: 			tmp=cdddph*cddddph-sdddph*sddddph; sdddph=sdddph*cddddph+cdddph*sddddph; cdddph=tmp;
xue@1: 		}
xue@1: 	}
xue@1: 	else
xue@1: 	{
xue@1: 		for (i=CountSt; i<CountEn; i++)
xue@1: 		{
xue@1: 			data[i]=exp(ea)*cdouble(cph, sph);
xue@1: 			ea=ea+dea; dea=dea+ddea; ddea=ddea+dddea; dddea+dddea+ddddea;
xue@1: 			tmp=cph*cdph-sph*sdph; sph=sph*cdph+cph*sdph; cph=tmp;
xue@1: 			tmp=cdph*cddph-sdph*sddph; sdph=sdph*cddph+cdph*sddph; cdph=tmp;
xue@1: 			tmp=cddph*cdddph-sddph*sdddph; sddph=sddph*cdddph+cddph*sdddph; cddph=tmp;
xue@1: 			tmp=cdddph*cddddph-sdddph*sddddph; sdddph=sdddph*cddddph+cdddph*sddddph; cdddph=tmp;
xue@1: 		}
xue@1: 	}
xue@1: 	ea=e0+CountEn*(e1+CountEn*(e2+CountEn*(e3+CountEn*e4)));
xue@1: 	ph=p0+CountEn*(p1+CountEn*(p2+CountEn*(p3+CountEn*p4)));
xue@1: }   //*/
xue@1: 
xue@1: 
Chris@5: /**
xue@1:   function Sinusoids: recursive harmonic multi-sinusoid generator
xue@1: 
xue@1:   In: st, en
xue@1:       M: number of partials
xue@1:       a3[M], a2[M], a1[M], a0[M]: trinomial coefficients for partial amplitudes
xue@1:       f3, f2, f1, f0: trinomial coefficients for fundamental frequency
xue@1:       ph[M]: partial phases at 0.
xue@1:       add: specifies if the resynthesized sinusoid is to be added to or to replace the content of output buffer
xue@1:   Out: data[st:en-1]: output buffer.
xue@1:        ph[M]: partial phases at en.
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void Sinusoids(int M, double* data, int st, int en, double* a3, double* a2, double* a1, double* a0, double f3, double f2, double f1, double f0, double* ph, bool add)
xue@1: {
xue@1:   double dph, ddph, dddph, ddddph;
xue@1:   double sdph, cdph, cdph2, sddph, cddph, sdddph, cdddph, sddddph, cddddph, sdmph, cdmph, sdm_1ph, cdm_1ph;
xue@1:   double p0, p1, p2, p3, p4, tmp, tmp2;
xue@1:   double *a=(double*)malloc8(sizeof(double)*M*6), *da=&a[M], *dda=&a[M*2], *ddda=&a[M*3],
xue@1:          *sph=&a[M*4], *cph=&a[M*5];
xue@1: 
xue@1:   for (int m=0; m<M; m++)
xue@1:   {
xue@1:     p0=ph[m], p1=2*M_PI*f0, p2=M_PI*f1, p3=2.0*M_PI*f2/3, p4=2.0*M_PI*f3/4;
xue@1:     if (st==0)
xue@1:     {
xue@1:       a[m]=a0[m]; da[m]=a1[m]+a2[m]+a3[m]; dda[m]=2*a2[m]+6*a3[m]; ddda[m]=6*a3[m];
xue@1:     }
xue@1:     else
xue@1:     {
xue@1:       a[m]=a0[m]+st*(a1[m]+st*(a2[m]+st*a3[m]));
xue@1:       da[m]=a1[m]+a2[m]+a3[m]+st*(2*a2[m]+3*a3[m]+st*3*a3[m]);
xue@1:       dda[m]=2*a2[m]+6*a3[m]+st*6*a3[m]; ddda[m]=6*a3[m];
xue@1:       ph[m]=p0+st*(p1+st*(p2+st*(p3+st*p4)));
xue@1:     }
xue@1:     sph[m]=sin(ph[m]), cph[m]=cos(ph[m]);
xue@1:     ph[m]=p0+(m+1)*en*(p1+en*(p2+en*(p3+en*p4)));
xue@1:   }
xue@1: 
xue@1:     if (st==0)
xue@1:     {
xue@1:       dph=p1+p2+p3+p4; ddph=2*p2+6*p3+14*p4; dddph=6*p3+36*p4; ddddph=24*p4;
xue@1:     }
xue@1:     else
xue@1:     {
xue@1:       dph=p1+p2+p3+p4+st*(2*p2+3*p3+4*p4+st*(3*p3+6*p4+st*4*p4));
xue@1:       ddph=2*p2+6*p3+14*p4+st*(6*p3+24*p4+st*12*p4);
xue@1:       dddph=6*p3+36*p4+st*24*p4; ddddph=24*p4;
xue@1:     }
xue@1:     sdph=sin(dph), cdph=cos(dph);
xue@1:     sddph=sin(ddph), cddph=cos(ddph);
xue@1:     sdddph=sin(dddph), cdddph=cos(dddph);
xue@1:     sddddph=sin(ddddph), cddddph=cos(ddddph);
xue@1: 
xue@1:   if (add)
xue@1:   {
xue@1:     for (int i=st; i<en; i++)
xue@1:     {
xue@1:       data[i]+=a[0]*cph[0]; a[0]+=da[0]; da[0]+=dda[0]; dda[0]+=ddda[0];
xue@1:       tmp=cph[0]*cdph-sph[0]*sdph; sph[0]=sph[0]*cdph+cph[0]*sdph; cph[0]=tmp;
xue@1:       cdm_1ph=1, sdm_1ph=0, cdmph=cdph, sdmph=sdph, cdph2=2*cdph;
xue@1: 
xue@1:       for (int m=1; m<M; m++)
xue@1:       {
xue@1:         data[i]+=a[m]*cph[m]; a[m]+=da[m]; da[m]+=dda[m]; dda[m]+=ddda[m];
xue@1: //        asm{mov ecx,m} asm{mov eax,a} asm{fld qword ptr [eax+ecx*8]} asm{mov edx,cph} asm{fld qword ptr [edx+ecx*8]} asm{fmul st(0),st(1)} asm{mov edx,data} asm{mov ebx,i} asm{fadd qword ptr [edx+ebx*8]} asm{fstp qword ptr [edx+ebx*8]} asm{mov edx,da} asm{fld qword ptr [edx+ecx*8]} asm{fadd st(1),st(0)} asm{mov ebx,dda} asm{fld qword ptr [ebx+ecx*8]} asm{fadd st(1),st(0)} asm{mov edi,ddda} asm{fadd qword ptr [edi+ecx*8]} asm{fstp qword ptr [ebx+ecx*8]} asm{fstp qword ptr [edx+ecx*8]} asm{fstp qword ptr [eax+ecx*8]}
xue@1:         tmp=cdmph, tmp2=sdmph;
xue@1:         cdmph=cdmph*cdph2-cdm_1ph; sdmph=sdmph*cdph2-sdm_1ph;
xue@1:         cdm_1ph=tmp, sdm_1ph=tmp2;
xue@1: 
xue@1:         tmp=cph[m]*cdmph-sph[m]*sdmph; sph[m]=sph[m]*cdmph+cph[m]*sdmph; cph[m]=tmp;
xue@1: //        asm{mov ecx,m} asm{mov eax,cph} asm{fld qword ptr [eax+ecx*8]} asm{mov edx,sph} asm{fld qword ptr [edx+ecx*8]} asm{fld st(1)} asm{fmul sdmph} asm{fld st(1)} asm{fmul sdmph} asm{fld cdmph} asm{fmul st(4),st(0)} asm{fmulp st(3),st(0)} asm{fsubp st(3),st(0)} asm{faddp} asm{fstp qword ptr [edx+ecx*8]} asm{fstp qword ptr [eax+ecx*8]}
xue@1:       }
xue@1: 
xue@1:       tmp=cdph*cddph-sdph*sddph; sdph=sdph*cddph+cdph*sddph; cdph=tmp;
xue@1:       tmp=cddph*cdddph-sddph*sdddph; sddph=sddph*cdddph+cddph*sdddph; cddph=tmp;
xue@1:       tmp=cdddph*cddddph-sdddph*sddddph; sdddph=sdddph*cddddph+cdddph*sddddph; cdddph=tmp;
xue@1:     }
xue@1:   }
xue@1:   else
xue@1:   {
xue@1:   }
xue@1:   free8(a);
xue@1: }//Sinusoids*/
xue@1: 
Chris@5: /**
xue@1:   function Sinusoid: synthesizes sinusoid piece from trinomial frequency and amplitude coefficients.
xue@1: 
xue@7:   In: T
xue@1:       aa, ab, ac, ad: trinomial coefficients of amplitude.
xue@1:       fa, fb, fc, fd: trinomial coefficients of frequency.
xue@1:       ph0, ph2: phase angles at 0 and CountEn.
xue@1:       add: specifies if the resynthesized sinusoid is to be added to or to replace the content of output buffer
xue@7:   Out: data[T]: output buffer.
xue@1: 
xue@1:   No return value.
xue@1: */
xue@7: void Sinusoid(int T, double* data, double aa, double ab, double ac, double ad,
xue@1:   double fa, double fb, double fc, double fd, double ph0, double ph2, bool add)
xue@1: {
xue@7:   double pend=ph0+2*M_PI*T*(fd+T*(fc/2+T*(fb/3+T*fa/4)));
xue@1:   int k=floor((pend-ph2)/2/M_PI+0.5);
xue@1:   double d=ph2-pend+2*M_PI*k;
xue@7:   double p=-2*d/T/T/T;
xue@7:   double q=3*d/T/T, a, ph;
xue@7:   for (int i=0; i<T; i++)
xue@1:   {
xue@1:     a=ad+i*(ac+i*(ab+i*aa)); if (a<0) a=0;
xue@1:     ph=ph0+2*M_PI*i*(fd+i*((fc/2)+i*((fb/3)+i*fa/4)))+i*i*(q+i*p);
xue@1:     if (add) data[i]+=a*cos(ph);
xue@1:     else data[i]=a*cos(ph);
xue@1:   }
xue@1: }//Sinusoid
xue@1: 
Chris@5: /**
xue@1:   function Sinusoid: synthesizes sinusoid piece from trinomial frequency and amplitude coefficients,
xue@1:   returning sinusoid coefficients instead of waveform.
xue@1: 
xue@7:   In: T
xue@1:       aa, ab, ac, ad: trinomial coefficients of amplitude (or log amplitude if LogA=true)
xue@1:       fa, fb, fc, fd: trinomial coefficients of frequency.
xue@1:       ph0, ph2: phase angles at 0 and CountEn.
xue@1:       LogA: specifies whether log amplitude or amplitude is a trinomial
xue@1:   Out: f[CountSt:CountEn-1], a[CountSt:CountEn-1], ph[CountSt:CountEn-1]: synthesized sinusoid parameters
xue@1:        da[CountSt:CountEn-1]: derivative of synthesized amplitude, optional
xue@1: 
xue@1:   No return value.
xue@1: */
xue@7: void Sinusoid(int T, double* f, double* a, double* ph, double* da,	double aa, double ab,
xue@1:   double ac, double ad,	double fa, double fb, double fc, double fd, double ph0, double ph2, bool LogA)
xue@1: {
xue@7:   double pend=ph0+2*M_PI*T*(fd+T*(fc/2+T*(fb/3+T*fa/4)));
xue@1:   int k=floor((pend-ph2)/2/M_PI+0.5);
xue@1:   double d=ph2-pend+2*M_PI*k;
xue@7:   double p=-2*d/T/T/T;
xue@7:   double q=3*d/T/T;
xue@7:   if (LogA) for (int i=0; i<T; i++)
xue@1: 	{
xue@1: 		a[i]=exp(ad+i*(ac+i*(ab+i*aa)));
xue@1: 		if (da) da[i]=a[i]*(ac+i*(2*ab+i*3*aa));
xue@1: 		f[i]=fd+i*(fc+i*(fb+i*fa))+i*(2*q+3*i*p)/(2*M_PI);
xue@1:     ph[i]=ph0+2*M_PI*i*(fd+i*((fc/2)+i*((fb/3)+i*fa/4)))+i*i*(q+i*p);
xue@1: 	}
xue@7:   else for (int i=0; i<T; i++)
xue@1: 	{
xue@1: 		a[i]=ad+i*(ac+i*(ab+i*aa));
xue@1: 		if (da) da[i]=ac+i*(2*ab+i*3*aa);
xue@1: 		f[i]=fd+i*(fc+i*(fb+i*fa))+i*(2*q+3*i*p)/(2*M_PI);
xue@1:     ph[i]=ph0+2*M_PI*i*(fd+i*((fc/2)+i*((fb/3)+i*fa/4)))+i*i*(q+i*p);
xue@1: 	}
xue@1: }//Sinusoid
xue@1: 
Chris@5: /**
xue@1:   function Sinusoid: generates trinomial frequency and phase with phase correction.
xue@1: 
xue@7:   In: T
xue@1:       fa, fb, fc, fd: trinomial coefficients of frequency.
xue@1:       ph0, ph2: phase angles at 0 and CountEn.
xue@7:   Out: f[T], ph[T]: output buffers holding frequency and phase.
xue@1: 
xue@1:   No return value.
xue@1: */
xue@7: void Sinusoid(int T, double* f, double* ph, double fa, double fb,
xue@1:   double fc, double fd, double ph0, double ph2)
xue@1: {
xue@7:   double pend=ph0+2*M_PI*T*(fd+T*(fc/2+T*(fb/3+T*fa/4)));
xue@1:   int k=floor((pend-ph2)/2/M_PI+0.5);
xue@1:   double d=ph2-pend+2*M_PI*k;
xue@7:   double p=-2*d/T/T/T;
xue@7:   double q=3*d/T/T;
xue@7:   for (int i=0; i<T; i++)
xue@1: 	{
xue@1: 		f[i]=fd+i*(fc+i*(fb+i*fa))+i*(2*q+3*i*p)/(2*M_PI);
xue@1:     ph[i]=ph0+2*M_PI*i*(fd+i*((fc/2)+i*((fb/3)+i*fa/4)))+i*i*(q+i*p);
xue@1: 	}
xue@1: }//Sinusoid
xue@1: 
Chris@5: /**
xue@1:   function SynthesizeSinusoid: synthesizes a time-varying sinusoid from a sequence of frequencies and amplitudes
xue@1: 
xue@1:   In: xs[Fr]: measurement points, should be integers although *xs has double type.
xue@1:       fs[Fr], as[Fr]: sequence of frequencies and amplitudes at xs[Fr]
xue@1:       phs[0]: initial phase angle at (int)xs[0].
xue@1:       dst, den: start and end time of synthesis, dst<=xs[0], den>=xs[Fr-1]
xue@1:       add: specifies if the resynthesized sinusoid is to be added to or to replace the content of output buffer
xue@1:   Out: xrec[0:den-dst-1]: output buffer hosting synthesized sinusoid from dst to den.
xue@1:        phs[Fr]: phase angles at xs[Fr]
xue@1: 
xue@1:   Returns pointer to xrec.
xue@1: */
xue@1: double* SynthesizeSinusoid(double* xrec, int dst, int den, double* phs, int Fr, double* xs, double* fs, double* as, bool add, bool* terminatetag)
xue@1: {
xue@1:   double *f3=new double[Fr*8], *f2=&f3[Fr], *f1=&f3[Fr*2], *f0=&f3[Fr*3],
xue@1:          *a3=&f3[Fr*4], *a2=&a3[Fr], *a1=&a3[Fr*2], *a0=&a3[Fr*3];
xue@1:   CubicSpline(Fr-1, f3, f2, f1, f0, xs, fs, 1, 1);
xue@1:   CubicSpline(Fr-1, a3, a2, a1, a0, xs, as, 1, 1);
xue@1:   double ph=phs[0];
xue@1:   for (int fr=0; fr<Fr-1; fr++)
xue@1:   {
xue@1:     phs[fr]=ph;
xue@1:     ALIGN8(Sinusoid(&xrec[(int)xs[fr]-dst], 0, xs[fr+1]-xs[fr], a3[fr], a2[fr], a1[fr], a0[fr], f3[fr], f2[fr], f1[fr], f0[fr], ph, add);)
xue@1:     if (terminatetag && *terminatetag) {delete[] f3; return 0;}
xue@1:   }
xue@11:   phs[Fr-1]=ph; ph=phs[Fr-2];
xue@1:   ALIGN8(Sinusoid(&xrec[(int)xs[Fr-2]-dst], xs[Fr-1]-xs[Fr-2], den-xs[Fr-2], a3[Fr-2], a2[Fr-2], a1[Fr-2], a0[Fr-2], f3[Fr-2], f2[Fr-2], f1[Fr-2], f0[Fr-2], ph, add);
xue@1:          Sinusoid(&xrec[(int)xs[0]-dst], dst-xs[0], 0, a3[0], a2[0], a1[0], a0[0], f3[0], f2[0], f1[0], f0[0], phs[0], add);)
xue@1:   delete[] f3;
xue@1:   return xrec;
xue@1: }//SynthesizeSinusoid
xue@1: 
Chris@5: /**
xue@1:   function ShiftTrinomial: shifts the origin of a trinomial from 0 to T
xue@1: 
xue@1:   In: a3, a2, a1, a0.
xue@1:   Out: b3, b2, b1, b0, so that a3*x^3+a2*x^2+a1*x+a0=b3(x-T)^3+b2(x-T)^2+b1(x-T)+b0
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void ShiftTrinomial(double T, double& b3, double& b2, double& b1, double& b0, double a3, double a2, double a1, double a0)
xue@1: {
xue@1:   b3=a3;
xue@1:   b2=a2+T*3*b3;
xue@1:   b1=a1+T*(2*b2-T*3*b3);
xue@1:   b0=a0+T*(b1-T*(b2-T*b3));
xue@1: }//ShiftTrinomial
xue@1: 
Chris@5: /**
xue@1:   function SynthesizeSinusoidP: synthesizes a time-varying sinusoid from a sequence of frequencies,
xue@1:   amplitudes and phase angles
xue@1: 
xue@1:   In: xs[Fr]: measurement points, should be integers although *xs has double type.
xue@1:       fs[Fr], as[Fr], phs[Fr]: sequence of frequencies, amplitudes and phase angles at xs[Fr]
xue@1:       dst, den: start and end time of synthesis, dst<=xs[0], den>=xs[Fr-1]
xue@1:       add: specifies if the resynthesized sinusoid is to be added to or to replace the content of output
xue@1:         buffer
xue@1:   Out: xrecm[0:den-dst-1]: output buffer hosting synthesized sinusoid from dst to den.
xue@1: 
xue@1:   Returns pointer to xrecm.
xue@1: */
xue@1: double* SynthesizeSinusoidP(double* xrecm, int dst, int den, double* phs, int Fr, double* xs, double* fs, double* as, bool add)
xue@1: {
xue@1:   double *f3=new double[Fr*8], *f2=&f3[Fr], *f1=&f3[Fr*2], *f0=&f3[Fr*3],
xue@1:          *a3=&f3[Fr*4], *a2=&a3[Fr], *a1=&a3[Fr*2], *a0=&a3[Fr*3];
xue@1:   CubicSpline(Fr-1, f3, f2, f1, f0, xs, fs, 1, 1);
xue@1:   CubicSpline(Fr-1, a3, a2, a1, a0, xs, as, 1, 1);
xue@7:   for (int fr=0; fr<Fr-1; fr++) Sinusoid(xs[fr+1]-xs[fr], &xrecm[(int)xs[fr]-dst], a3[fr], a2[fr], a1[fr], a0[fr], f3[fr], f2[fr], f1[fr], f0[fr], phs[fr], phs[fr+1], add);
xue@1:   double tmpph=phs[0]; Sinusoid(&xrecm[(int)xs[0]-dst], dst-xs[0], 0, 0, 0, 0, a0[0], f3[0], f2[0], f1[0], f0[0], tmpph, add);
xue@1:   //extend the trinomials on [xs[Fr-2], xs[Fr-1]) based at xs[Fr-2] to beyond xs[Fr-1] based at xs[Fr-1].
xue@1:   tmpph=phs[Fr-1];
xue@1:   ShiftTrinomial(xs[Fr-1]-xs[Fr-2], f3[Fr-1], f2[Fr-1], f1[Fr-1], f0[Fr-1], f3[Fr-2], f2[Fr-2], f1[Fr-2], f0[Fr-2]);
xue@1:   ShiftTrinomial(xs[Fr-1]-xs[Fr-2], a3[Fr-1], a2[Fr-1], a1[Fr-1], a0[Fr-1], a3[Fr-2], a2[Fr-2], a1[Fr-2], a0[Fr-2]);
xue@1:   Sinusoid(&xrecm[(int)xs[Fr-1]-dst], 0, den-xs[Fr-1], 0, 0, 0, a0[Fr-1], f3[Fr-1], f2[Fr-1], f1[Fr-1], f0[Fr-1], tmpph, add);
xue@1:   delete[] f3;
xue@1:   return xrecm;
xue@1: }//SynthesizeSinusoidP