comparison sinsyn.cpp @ 7:9b1c0825cc77

TFileStream; redefining some function arguments

 #include "splines.h"
 
 /** \file sinsyn.h */
 
 //---------------------------------------------------------------------------
+/**
+  function CosSin: recursive cos-sin generator with trinomial frequency
+
+  In: CountSt, CountEn
+      f3, f2, f1, f0: trinomial coefficients of frequency
+      ph: initial phase angle at 0 (NOT at CountSt)
+  Out: datar[CountSt:CountEn-1], datai[CountSt:CountEn-1]: synthesized pair of cosine and sine functions
+       ph: phase angle at CountEn
+
+  No return value.
+*/
+void CosSin(double* datar, double* datai, int CountSt, int CountEn, double f3, double f2, double f1, double f0, double &ph)
+{
+  int i;
+  double dph, ddph, dddph, ddddph,
+         sph, cph, sdph, cdph, sddph, cddph, sdddph, cdddph, sddddph, cddddph,
+         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;
+  if (CountSt==0)
+  {
+    dph=p1+p2+p3+p4; ddph=2*p2+6*p3+14*p4; dddph=6*p3+36*p4; ddddph=24*p4;
+  }
+  else
+  {
+    ph=p0+CountSt*(p1+CountSt*(p2+CountSt*(p3+CountSt*p4)));
+    dph=p1+p2+p3+p4+CountSt*(2*p2+3*p3+4*p4+CountSt*(3*p3+6*p4+CountSt*4*p4));
+    ddph=2*p2+6*p3+14*p4+CountSt*(6*p3+24*p4+CountSt*12*p4);
+    dddph=6*p3+36*p4+CountSt*24*p4; ddddph=24*p4;
+  }
+  sph=sin(ph), cph=cos(ph);
+  sdph=sin(dph), cdph=cos(dph);
+  sddph=sin(ddph), cddph=cos(ddph);
+  sdddph=sin(dddph), cdddph=cos(dddph);
+  sddddph=sin(ddddph), cddddph=cos(ddddph);
+
+  for (i=CountSt; i<CountEn; i++)
+  {
+    datar[i]=cph; datai[i]=sph;
+    tmp=cph*cdph-sph*sdph; sph=sph*cdph+cph*sdph; cph=tmp;
+    tmp=cdph*cddph-sdph*sddph; sdph=sdph*cddph+cdph*sddph; cdph=tmp;
+    tmp=cddph*cdddph-sddph*sdddph; sddph=sddph*cdddph+cddph*sdddph; cddph=tmp;
+    tmp=cdddph*cddddph-sdddph*sddddph; sdddph=sdddph*cddddph+cdddph*sddddph; cdddph=tmp;
+  }
+  ph=p0+CountEn*(p1+CountEn*(p2+CountEn*(p3+CountEn*p4)));
+}//CosSin*/
+
 /**
   function Sinuoid:	original McAuley-Quatieri synthesizer interpolation between two measurement points.
 
   In: T: length from measurement point 1 to measurement point 2
       a1, f1, p2: amplitude, frequency and phase angle at measurement point 1

   }
   p1=p1+2*M_PI*i*(fd+i*((fc/2)+i*((fb/3)+i*fa/4)));
 }//Sinusoid
 
 /**
+  function Sinusoid_direct: synthesizes sinusoid over [CountSt, CountEn) from tronomial coefficients of
+  amplitude and frequency, direct implementation returning sinusoid coefficients rather than waveform.
+
+  In: CountSt, CountEn
+      aa, ab, ac, ad: trinomial coefficients of amplitude
+      fa, fb, fc, fd: trinomial coefficients of frequency
+      p1: initial phase angle at 0 (NOT at CountSt)
+  Out: f[CountSt:CountEn-1], a[:], ph[:], da[:]: output buffers.
+       p1: phase angle at CountEn
+
+  No return value.
+*/
+void Sinusoid_direct(double* f, double* a, double* ph, double* da, int CountSt, int CountEn, double aa, double ab, double ac, double ad,
+  double fa, double fb, double fc, double fd, double &p1, bool LogA)
+{
+  int i;
+  if (LogA)
+  {
+    for (i=CountSt; i<CountEn; i++)
+    {
+       a[i]=exp(ad+i*(ac+i*(ab+i*aa)));
+       f[i]=fd+i*(fc+i*(fb+i*fa));
+       ph[i]=p1+2*M_PI*i*(fd+i*((fc/2)+i*((fb/3)+i*fa/4)));
+       da[i]=a[i]*(ac+i*(2*ab+3*i*aa));
+    }
+  }
+  else
+  {
+    for (i=CountSt; i<CountEn; i++)
+    {
+       a[i]=ad+i*(ac+i*(ab+i*aa));
+       f[i]=fd+i*(fc+i*(fb+i*fa));
+       ph[i]=p1+2*M_PI*i*(fd+i*((fc/2)+i*((fb/3)+i*fa/4)));
+       da[i]=ac+i*(2*ab+3*i*aa);
+    }
+  }
+  p1=p1+2*M_PI*i*(fd+i*((fc/2)+i*((fb/3)+i*fa/4)));
+}//Sinusoid
+
+/**
+  function Sinusoid_direct: synthesizes sinusoid over [CountSt, CountEn) from tronomial coefficients of
+  frequency, direct implementation returning frequency and phase rather than waveform.
+
+  In: CountSt, CountEn
+      fa, fb, fc, fd: trinomial coefficients of frequency
+      p1: initial phase angle at 0 (NOT at CountSt)
+  Out: f[CountSt:CountEn-1], ph[CountSt:CountEn-1]: output buffers.
+       p1: phase angle at CountEn
+
+  No return value.
+*/
+void Sinusoid_direct(double* f, double* ph, int CountSt, int CountEn, double fa, double fb, double fc,
+      double fd, double &p1)
+{
+  int i;
+  for (i=CountSt; i<CountEn; i++)
+  {
+     f[i]=fd+i*(fc+i*(fb+i*fa));
+     ph[i]=p1+2*M_PI*i*(fd+i*((fc/2)+i*((fb/3)+i*fa/4)));
+  }
+  p1=p1+2*M_PI*i*(fd+i*((fc/2)+i*((fb/3)+i*fa/4)));
+}//Sinusoid
+
+/**
   function Sinusoid: synthesizes sinusoid over [CountSt, CountEn) from tronomial coefficients of
   amplitude and frequency, recursive implementation.
 
   In: CountSt, CountEn
       a3, a2, a1, a0: trinomial coefficients of amplitude

 
 /**
   function SinusoidExpA: synthesizes complex sinusoid whose log amplitude and frequency are trinomials
   with phase angle specified at both ends.
 
-  In: CountSt, CountEn
+  In: T
       a3, a2, a1, a0: trinomial coefficients for log amplitude
       omg3, omg2, omg1, omg0: trinomial coefficients for angular frequency
       ph0, ph2: phase angles at 0 and CountEn.
       add: specifies if the resynthesized sinusoid is to be added to or to replace the content of output buffer
-  Out: data[CountSt:CountEn-1]: output buffer.
-
-  No return value.
-*/
-void SinusoidExpA(cdouble* data, int CountSt, int CountEn, double a3, double a2, double a1, double a0,
+  Out: data[T]: output buffer.
+
+  No return value.
+*/
+void SinusoidExpA(int T, cdouble* data, double a3, double a2, double a1, double a0,
 	double omg3, double omg2, double omg1, double omg0, double ph0, double ph2, bool add)
 {
   double p0=ph0, p1=omg0, p2=0.5*omg1, p3=omg2/3, p4=omg3/4;
-	double pend=p0+CountEn*(p1+CountEn*(p2+CountEn*(p3+CountEn*p4)));
+  double pend=p0+T*(p1+T*(p2+T*(p3+T*p4)));
 
 	int k=floor((pend-ph2)/2/M_PI+0.5);
 	double d=ph2-pend+2*M_PI*k;
-	double _p=-2*d/CountEn/CountEn/CountEn;
-	double _q=3*d/CountEn/CountEn;
-
-	if (add) for (int i=CountSt; i<CountEn; i++)
+  double _p=-2*d/T/T/T;
+  double _q=3*d/T/T;
+
+  if (add) for (int i=0; i<T; i++)
 	{
 		double lea=a0+i*(a1+i*(a2+i*a3));
 		double lph=p0+i*(p1+i*(p2+i*(p3+i*p4)));
     data[i]+=exp(cdouble(lea, lph+(i*i*(_q+i*_p))));
 	}
-	else for (int i=CountSt; i<CountEn; i++)
+  else for (int i=0; i<T; i++)
 	{
 		double lea=a0+i*(a1+i*(a2+i*a3));
 		double lph=p0+i*(p1+i*(p2+i*(p3+i*p4)));
 		data[i]=exp(cdouble(lea, lph+(i*i*(_q+i*_p))));
 	}

 	}
 	ea=e0+CountEn*(e1+CountEn*(e2+CountEn*(e3+CountEn*e4)));
 	ph=p0+CountEn*(p1+CountEn*(p2+CountEn*(p3+CountEn*p4)));
 }   //*/
 
-/**
-  function Sinusoid: recursive cos-sin generator with trinomial frequency
-
-  In: CountSt, CountEn
-      f3, f2, f1, f0: trinomial coefficients of frequency
-      ph: initial phase angle at 0 (NOT at CountSt)
-  Out: datar[CountSt:CountEn-1], datai[CountSt:CountEn-1]: synthesized pair of cosine and sine functions
-       ph: phase angle at CountEn
-
-  No return value.
-*/
-void Sinusoid(double* datar, double* datai, int CountSt, int CountEn, double f3, double f2, double f1, double f0, double &ph)
-{
-  int i;
-  double dph, ddph, dddph, ddddph,
-         sph, cph, sdph, cdph, sddph, cddph, sdddph, cdddph, sddddph, cddddph,
-         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;
-  if (CountSt==0)
-  {
-    dph=p1+p2+p3+p4; ddph=2*p2+6*p3+14*p4; dddph=6*p3+36*p4; ddddph=24*p4;
-  }
-  else
-  {
-    ph=p0+CountSt*(p1+CountSt*(p2+CountSt*(p3+CountSt*p4)));
-    dph=p1+p2+p3+p4+CountSt*(2*p2+3*p3+4*p4+CountSt*(3*p3+6*p4+CountSt*4*p4));
-    ddph=2*p2+6*p3+14*p4+CountSt*(6*p3+24*p4+CountSt*12*p4);
-    dddph=6*p3+36*p4+CountSt*24*p4; ddddph=24*p4;
-  }
-  sph=sin(ph), cph=cos(ph);
-  sdph=sin(dph), cdph=cos(dph);
-  sddph=sin(ddph), cddph=cos(ddph);
-  sdddph=sin(dddph), cdddph=cos(dddph);
-  sddddph=sin(ddddph), cddddph=cos(ddddph);
-
-  for (i=CountSt; i<CountEn; i++)
-  {
-    datar[i]=cph; datai[i]=sph;
-    tmp=cph*cdph-sph*sdph; sph=sph*cdph+cph*sdph; cph=tmp;
-    tmp=cdph*cddph-sdph*sddph; sdph=sdph*cddph+cdph*sddph; cdph=tmp;
-    tmp=cddph*cdddph-sddph*sdddph; sddph=sddph*cdddph+cddph*sdddph; cddph=tmp;
-    tmp=cdddph*cddddph-sdddph*sddddph; sdddph=sdddph*cddddph+cdddph*sddddph; cdddph=tmp;
-  }
-  ph=p0+CountEn*(p1+CountEn*(p2+CountEn*(p3+CountEn*p4)));
-}//Sinusoid*/
 
 /**
   function Sinusoids: recursive harmonic multi-sinusoid generator
 
   In: st, en

 }//Sinusoids*/
 
 /**
   function Sinusoid: synthesizes sinusoid piece from trinomial frequency and amplitude coefficients.
 
-  In: CountSt, CountEn
+  In: T
       aa, ab, ac, ad: trinomial coefficients of amplitude.
       fa, fb, fc, fd: trinomial coefficients of frequency.
       ph0, ph2: phase angles at 0 and CountEn.
       add: specifies if the resynthesized sinusoid is to be added to or to replace the content of output buffer
-  Out: data[CountSt:CountEn-1]: output buffer.
-
-  No return value.
-*/
-void Sinusoid(double* data, int CountSt, int CountEn, double aa, double ab, double ac, double ad,
+  Out: data[T]: output buffer.
+
+  No return value.
+*/
+void Sinusoid(int T, double* data, double aa, double ab, double ac, double ad,
   double fa, double fb, double fc, double fd, double ph0, double ph2, bool add)
 {
-  double pend=ph0+2*M_PI*CountEn*(fd+CountEn*(fc/2+CountEn*(fb/3+CountEn*fa/4)));
+  double pend=ph0+2*M_PI*T*(fd+T*(fc/2+T*(fb/3+T*fa/4)));
   int k=floor((pend-ph2)/2/M_PI+0.5);
   double d=ph2-pend+2*M_PI*k;
-  double p=-2*d/CountEn/CountEn/CountEn;
-  double q=3*d/CountEn/CountEn, a, ph;
-  for (int i=CountSt; i<CountEn; i++)
+  double p=-2*d/T/T/T;
+  double q=3*d/T/T, a, ph;
+  for (int i=0; i<T; i++)
   {
     a=ad+i*(ac+i*(ab+i*aa)); if (a<0) a=0;
     ph=ph0+2*M_PI*i*(fd+i*((fc/2)+i*((fb/3)+i*fa/4)))+i*i*(q+i*p);
     if (add) data[i]+=a*cos(ph);
     else data[i]=a*cos(ph);

 
 /**
   function Sinusoid: synthesizes sinusoid piece from trinomial frequency and amplitude coefficients,
   returning sinusoid coefficients instead of waveform.
 
-  In: CountSt, CountEn
+  In: T
       aa, ab, ac, ad: trinomial coefficients of amplitude (or log amplitude if LogA=true)
       fa, fb, fc, fd: trinomial coefficients of frequency.
       ph0, ph2: phase angles at 0 and CountEn.
       LogA: specifies whether log amplitude or amplitude is a trinomial
   Out: f[CountSt:CountEn-1], a[CountSt:CountEn-1], ph[CountSt:CountEn-1]: synthesized sinusoid parameters
        da[CountSt:CountEn-1]: derivative of synthesized amplitude, optional
 
   No return value.
 */
-void Sinusoid(double* f, double* a, double* ph, double* da,	int CountSt, int CountEn, double aa, double ab,
+void Sinusoid(int T, double* f, double* a, double* ph, double* da,	double aa, double ab,
   double ac, double ad,	double fa, double fb, double fc, double fd, double ph0, double ph2, bool LogA)
 {
-  double pend=ph0+2*M_PI*CountEn*(fd+CountEn*(fc/2+CountEn*(fb/3+CountEn*fa/4)));
+  double pend=ph0+2*M_PI*T*(fd+T*(fc/2+T*(fb/3+T*fa/4)));
   int k=floor((pend-ph2)/2/M_PI+0.5);
   double d=ph2-pend+2*M_PI*k;
-  double p=-2*d/CountEn/CountEn/CountEn;
-  double q=3*d/CountEn/CountEn;
-	if (LogA) for (int i=CountSt; i<CountEn; i++)
+  double p=-2*d/T/T/T;
+  double q=3*d/T/T;
+  if (LogA) for (int i=0; i<T; i++)
 	{
 		a[i]=exp(ad+i*(ac+i*(ab+i*aa)));
 		if (da) da[i]=a[i]*(ac+i*(2*ab+i*3*aa));
 		f[i]=fd+i*(fc+i*(fb+i*fa))+i*(2*q+3*i*p)/(2*M_PI);
     ph[i]=ph0+2*M_PI*i*(fd+i*((fc/2)+i*((fb/3)+i*fa/4)))+i*i*(q+i*p);
 	}
-	else for (int i=CountSt; i<CountEn; i++)
+  else for (int i=0; i<T; i++)
 	{
 		a[i]=ad+i*(ac+i*(ab+i*aa));
 		if (da) da[i]=ac+i*(2*ab+i*3*aa);
 		f[i]=fd+i*(fc+i*(fb+i*fa))+i*(2*q+3*i*p)/(2*M_PI);
     ph[i]=ph0+2*M_PI*i*(fd+i*((fc/2)+i*((fb/3)+i*fa/4)))+i*i*(q+i*p);

 }//Sinusoid
 
 /**
   function Sinusoid: generates trinomial frequency and phase with phase correction.
 
-  In: CountSt, CountEn
+  In: T
       fa, fb, fc, fd: trinomial coefficients of frequency.
       ph0, ph2: phase angles at 0 and CountEn.
-  Out: f[CountSt:CountEn-1], ph[CountSt:CountEn-1]: output buffers holding frequency and phase.
-
-  No return value.
-*/
-void Sinusoid(double* f, double* ph, int CountSt, int CountEn, double fa, double fb,
+  Out: f[T], ph[T]: output buffers holding frequency and phase.
+
+  No return value.
+*/
+void Sinusoid(int T, double* f, double* ph, double fa, double fb,
   double fc, double fd, double ph0, double ph2)
 {
-  double pend=ph0+2*M_PI*CountEn*(fd+CountEn*(fc/2+CountEn*(fb/3+CountEn*fa/4)));
+  double pend=ph0+2*M_PI*T*(fd+T*(fc/2+T*(fb/3+T*fa/4)));
   int k=floor((pend-ph2)/2/M_PI+0.5);
   double d=ph2-pend+2*M_PI*k;
-  double p=-2*d/CountEn/CountEn/CountEn;
-  double q=3*d/CountEn/CountEn;
-	for (int i=CountSt; i<CountEn; i++)
+  double p=-2*d/T/T/T;
+  double q=3*d/T/T;
+  for (int i=0; i<T; i++)
 	{
 		f[i]=fd+i*(fc+i*(fb+i*fa))+i*(2*q+3*i*p)/(2*M_PI);
     ph[i]=ph0+2*M_PI*i*(fd+i*((fc/2)+i*((fb/3)+i*fa/4)))+i*i*(q+i*p);
 	}
 }//Sinusoid

 {
   double *f3=new double[Fr*8], *f2=&f3[Fr], *f1=&f3[Fr*2], *f0=&f3[Fr*3],
          *a3=&f3[Fr*4], *a2=&a3[Fr], *a1=&a3[Fr*2], *a0=&a3[Fr*3];
   CubicSpline(Fr-1, f3, f2, f1, f0, xs, fs, 1, 1);
   CubicSpline(Fr-1, a3, a2, a1, a0, xs, as, 1, 1);
-  for (int fr=0; fr<Fr-1; fr++) Sinusoid(&xrecm[(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], phs[fr], phs[fr+1], add);
+  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);
   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);
   //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].
   tmpph=phs[Fr-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]);
   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]);