comparison sinest.cpp @ 5:5f3c32dc6e17

* Adjust comment syntax to permit Doxygen to generate HTML documentation; add Doxyfile

 #include "opt.h"
 #include "sinsyn.h"
 #include "splines.h"
 #include "windowfunctions.h"
 
+/** \file sinest.h */
+
 //---------------------------------------------------------------------------
-/*
+/**
   function dsincd_unn: derivative of unnormalized discrete sinc function
 
   In: x, scale N
 
   Returns the derivative of sincd_unn(x, N)

       r=0;
   }
   return r;
 }//dsincd_unn
 
-/*
+/**
   function ddsincd_unn: 2nd-order derivative of unnormalized discrete sinc function
 
   In: x, scale (equivalently, window size) N
 
   Returns the 2nd-order derivative of sincd_unn(x, N)

   }
   return r;
 }//ddsincd_unn
 
 //---------------------------------------------------------------------------
-/*
+/**
   function Window: calculates the cosine-family-windowed spectrum of a complex sinusoid on [0:N-1] at
   frequency f bins with zero central phase.
 
   In: f: frequency, in bins
       N: window size

     }
   }
   return x;
 }//Window
 
-/*
+/**
   function dWindow: calculates the cosine-family-windowed spectrum and its derivative of a complex
   sinusoid on [0:N-1] at frequency f bins with zero central phase.
 
   In: f: frequency, in bins
       N: window size

       }
     }
   }
 }//dWindow
 
-/*
+/**
   function ddWindow: calculates the cosine-family-windowed spectrum and its 1st and 2nd derivatives of
   a complex sinusoid on [0:N-1] at frequency f bins with zero central phase.
 
   In: f: frequency, in bins
       N: window size

     }
   }
 }//ddWindow
 
 //---------------------------------------------------------------------------
-/*
+/**
   function IPWindow: computes the truncated inner product of a windowed spectrum with that of a sinusoid
   at reference frequency f.
 
   In: x[0:N-1]: input spectrum
       f: reference frequency, in bins

 {
   struct l_ip {int N; int k1; int k2; int M; double* c; double iH2; cdouble* x; double dipwindow; double ipwindow;} *p=(l_ip *)params;
   return IPWindow(f, p->x, p->N, p->M, p->c, p->iH2, p->k1, p->k2, true);
 }//IPWindow
 
-/*
+/**
   function ddIPWindow: computes the norm of the truncated inner product of a windowed spectrum with
   that of a sinusoid at reference frequency f, as well as its 1st and 2nd derivatives.
 
   In: x[0:N-1]: input spectrum
       f: reference frequency, in bins

 	struct l_ip {int N; int k1; int k2; int M; double* c; double iH2; cdouble* x; double dipwindow; double ipwindow;} *p=(l_ip *)params;
 	return ddIPWindow(f, p->x, p->N, p->M, p->c, p->iH2, p->k1, p->k2, p->dipwindow, p->ipwindow);
 }//ddIPWindow
 
 //---------------------------------------------------------------------------
-/*
+/**
   function IPWindowC: computes the truncated inner product of a windowed spectrum with that of a
   sinusoid at reference frequency f.
 
   In: x[0:N-1]: input spectrum
       f: reference frequency, in bins

 	result*=iH2;
 	return result;
 }//IPWindowC
 
 //---------------------------------------------------------------------------
-/*
+/**
   function sIPWindow: computes the total energy of truncated inner products between multiple windowed
   spectra and that of a sinusoid at a reference frequency f. This does not consider phase alignment
   between the spectra, supposedly measured at a sequence of known instants.
 
   In: x[L][N]: input spectra

 {
 	struct l_ip {int N; int k1; int k2; int M; double* c; double iH2; int Fr; cdouble** x; double dipwindow; double ipwindow; cdouble* lmd;} *p=(l_ip *)params;
 	return sIPWindow(f, p->Fr, p->x, p->N, p->M, p->c, p->iH2, p->k1, p->k2, p->lmd);
 }//sIPWindow
 
-/*
+/**
   function dsIPWindow: computes the total energy of truncated inner products between multiple windowed
   spectra and that of a sinusoid at a reference frequency f, as well as its derivative. This does not
   consider phase synchronization between the spectra, supposedly measured at a sequence of known
   instants.
 

 {
 	struct l_ip1 {int N; int k1; int k2; int M; double* c; double iH2; int Fr; cdouble** x; double sip;} *p=(l_ip1 *)params;
 	return dsIPWindow(f, p->Fr, p->x, p->N, p->M, p->c, p->iH2, p->k1, p->k2, p->sip);
 }//dsIPWindow
 
-/*
+/**
   function dsdIPWindow_unn: computes the energy of unnormalized truncated inner products between a given
   windowed spectrum and that of a sinusoid at a reference frequency f, as well as its 1st and 2nd
   derivatives. "Unnormalized" indicates that the inner product cannot be taken as the actual amplitude-
   phase factor of a sinusoid, but deviate from that by an unspecified factor.
 

   dsipwindow=dR2;
   if (w_unn) w_unn->x=rr, w_unn->y=ri;
 	return ddR2;
 }//ddsIPWindow_unn
 
-/*
+/**
   function ddsIPWindow: computes the total energy of truncated inner products between multiple windowed
   spectra and that of a sinusoid at a reference frequency f, as well as its 1st and 2nd derivatives.
   This does not consider phase synchronization between the spectra, supposedly measured at a sequence
   of known instants.
 

   struct l_ip1 {int N; int k1; int k2; int M; double* c; double iH2; int Fr; cdouble** x; double dsip; double sip;} *p=(l_ip1 *)params;
 	return ddsIPWindow(f, p->Fr, p->x, p->N, p->M, p->c, p->iH2, p->k1, p->k2, p->dsip, p->sip);
 }//ddsIPWindow
 
 //---------------------------------------------------------------------------
-/*
+/**
   function sIPWindowC: computes the total energy of truncated inner products between multiple frames of
   a spectrogram and multiple frames of a spectrogram of a sinusoid at a reference frequency f.
 
   In: x[L][N]: the spectrogram
       offst_rel: frame offset, relative to frame size

 {
   struct l_ip {int N; int k1; int k2; int M; double* c; double iH2; int L; double offst_rel; cdouble** x; double dipwindow; double ipwindow;} *p=(l_ip *)params;
 	return sIPWindowC(f, p->L, p->offst_rel, p->x, p->N, p->M, p->c, p->iH2, p->k1, p->k2);
 }//sIPWindowC
 
-/*
+/**
   function dsIPWindowC: computes the total energy of truncated inner products between multiple frames of
   a spectrogram and multiple frames of a spectrogram of a sinusoid at a reference frequency f, together
   with its derivative.
 
   In: x[L][N]: the spectrogram

 {
   struct l_ip {int N; int k1; int k2; int M; double* c; double iH2; int L; double offst_rel; cdouble** x; double sip;} *p=(l_ip *)params;
 	return dsIPWindowC(f, p->L, p->offst_rel, p->x, p->N, p->M, p->c, p->iH2, p->k1, p->k2, p->sip);
 }//dsIPWindowC
 
-/*
+/**
   function ddsIPWindowC: computes the total energy of truncated inner products between multiple frames
   of a spectrogram and multiple frames of a spectrogram of a sinusoid at a reference frequency f,
   together with its 1st and 2nd derivatives.
 
   In: x[L][N]: the spectrogram

   f & epf are given/returned in bins.
 
   Further reading: "Least-square-error estimation of sinusoids.pdf"
 */
 
-/*
+/**
   function LSESinusoid: LSE estimation of the predominant stationary sinusoid.
 
   In: x[N]: windowed spectrum
       B: spectral truncation half width, in bins.
       M, c[], iH2: cosine-family window specification parameters

   else
     a=p.hxpeak;
   pp=IPWindow(f, x, N, M, c, iH2, p.k1, p.k2, false);
 }//LSESinusoid
 
-/*
+/**
   function LSESinusoid: LSE estimation of stationary sinusoid predominant within [f1, f2].
 
   In: x[N]: windowed spectrum
       [f1, f2]: frequency range
       B: spectral truncation half width, in bins.

     a=p.hxpeak;
   pp=IPWindow(f, x, N, M, c, iH2, p.k1, p.k2, false);
   return f;
 }//LSESinusoid
 
-/*
+/**
   function LSESinusoid: LSE estimation of stationary sinusoid near a given initial frequency within [f1,
   f2].
 
   In: x[N]: windowed spectrum
       f: initial frequency, in bins

     }
   }
   return result;
 }//LSESinusoid
 
-/*
+/**
   function LSESinusoidMP: LSE estimation of a stationary sinusoid from multi-frames spectrogram without
   considering phase-frequency consistency across frames.
 
   In: x[Fr][N]: spectrogram
       f: initial frequency, in bins

     }
   }
   return errf;
 }//LSESinusoidMP
 
-/*
+/**
   function LSESinusoidMP: LSE estimation of a stationary sinusoid from multi-frames spectrogram without
   considering phase-frequency consistency across frames.
 
   In: x[Fr][N]: spectrogram
       f: initial frequency, in bins

   }
   return errf;
 }//LSESinusoidMPC
 
 //---------------------------------------------------------------------------
-/*
+/**
   function IPMulti: least square estimation of multiple sinusoids, given their frequencies and an energy
   suppression index of eps, i.e. the least square error is minimized with an additional eps*||lmd||^2
   term.
 
   In: x[Wid]: spectrum

   for (int i=0; i<I; i++) whw[i][i]+=eps;
   GECP(I, lmd, whw, whx);
   delete List;
 }//IPMulti
 
-/*
+/**
   function IPMulti: least square estimation of multiple sinusoids, given their frequencies and an energy
   suppression index of eps, and optionally returns residue and sensitivity indicators for each sinusoid.
 
   In: x[Wid]: spectrum
       f[I]: frequencies

     }
   }
   delete List;
 }//IPMulti
 
-/*
+/**
   function IPMultiSens: computes the sensitivity of the least square estimation of multiple sinusoids given
   their frequencies .
 
   In: f[I]: frequencies
       M, c[]: cosine-family window specification parameters

     sens[i]=0; for (int k=0; k<K; k++) sens[i]+=~u[i][k]; sens[i]=sqrt(sens[i]);
   }
   delete List;
 }//IPMultiSens
 
-/*
+/**
   function IPMulti: least square estimation of multi-sinusoids with GIVEN frequencies. This version
   operates in groups at least B bins from each other, rather than LSE all frequencies together.
 
   In: x[Wid]: spectrum
       f[I]: frequencies, must be ordered low to high.

     i++;
   }
   return 0;
 }//IPMulti
 
-/*
+/**
   function IPMulti_Direct: LSE estimation of multiple sinusoids given frequencies AND PHASES (direct
   method)
 
   In: x[Wid]: spectrum
       f[I], ph[I]: frequencies and phase angles.

   double result=Inner(hWid, r, r).x;
   delete List;
   return result;
 }//IPMulti_Direct
 
-/*
+/**
   function IPMulti_GS: LSE estimation of multiple sinusoids given frequencies AND PHASES (Gram-Schmidt method)
 
   In: x[Wid]: spectrum
       f[I], ph[I]: frequencies and phase angles.
       B: spectral truncation, in bins; sinusoids over 3B bins apart are regarded non-interfering

   double result=Inner(hWid, r, r).x;
   delete List;
 	return result;
 }//IPMulti_GS
 
-/*
+/**
   function IPMulti: LSE estimation of I sinusoids given frequency and phase and J sinusoids given
   frequency only
 
   In: x[Wid]: spectrum
       f[I+J], ph[I]: frequencies and phase angles

     Routines for estimation two sinusoids with 1 fixed and 1 flexible frequency
 
     Further reading: "LSE estimation for 2 sinusoids with 1 at a fixed frequency.pdf"
 */
 
-/*
+/**
   function WindowDuo: calcualtes the square norm of the inner product between windowed spectra of two
   sinusoids at frequencies f1 and f2, df=f1-f2.
 
   In: df: frequency difference, in bins
       N: DFT size

 	if (w) w->x=wr, w->y=wi;
 	double result=wr*wr+wi*wi;
 	return result;
 }//WindowDuo
 
-/*
+/**
   function ddWindowDuo: calcualtes the square norm of the inner product between windowed spectra of two
   sinusoids at frequencies f1 and f2, df=f1-f2, with its 1st and 2nd derivatives
 
   In: df: frequency difference, in bins
       N: DFT size

   if (w) w->x=wr, w->y=wi;
 	double ddwindow=2*(wr*ddwr+dwr*dwr+wi*ddwi+dwi*dwi);
   return ddwindow;
 }//ddWindowDuo
 
-/*
+/**
   function sIPWindowDuo: calculates the square norm of the orthogonal projection of a windowed spectrum
   onto the linear span of the windowed spectra of two sinusoids at reference frequencies f1 and f2.
 
   In: x[N]: spectrum
       f1, f2: reference frequencies.

 	struct l_ip {int N; int k1; int k2; double* c; double* d; int M; double iH2; cdouble* x; double f1; double dipwindow; double ipwindow;} *p=(l_ip *)params;
 	cdouble r1, r2;
 	return sIPWindowDuo(p->f1, f2, p->x, p->N, p->c, p->d, p->M, p->iH2, p->k1, p->k2, r1, r2);
 }//sIPWindowDuo
 
-/*
+/**
   function ddsIPWindowDuo: calculates the square norm, and its 1st and 2nd derivatives against f2,, of
   the orthogonal projection of a windowed spectrum onto the linear span of the windowed spectra of two
   sinusoids at reference frequencies f1 and f2.
 
   In: x[N]: spectrum

   ddsIPWindowDuo(ddsip2, p->f1, f2, p->x, p->N, p->c, p->d, p->M, p->iH2, p->k1, p->k2, r1, r2);
   p->dipwindow=ddsip2[1], p->ipwindow=ddsip2[2];
   return ddsip2[0];
 }//ddsIPWindowDuo
 
-/*
+/**
   function LSEDuo: least-square estimation of two sinusoids of which one has a fixed frequency
 
   In: x[N]: the windowed spectrum
       f1: the fixed frequency
       f2: initial value of the flexible frequency

 
     Further reading: A. R?bel, ¡°Estimating partial frequency and frequency slope using reassignment
     operators,¡± in Proc. ICMC¡¯02. G?teborg. 2002.
 */
 
-/*
+/**
   function CDFTW: single-frequency windowed DTFT, centre-aligned
 
   In: data[Wid]: waveform data x
       win[Wid+1]: window function
       k: frequency, in bins, where bin=1/Wid

     tmp.rotate(ph);
     X+=tmp;
   }
 }//CDFTW
 
-/*
+/**
   function CuDFTW: single-frequency windowed DTFT of t*data[t], centre-aligned
 
   In: data[Wid]: waveform data x
       wid[Wid+1]: window function
       k: frequency, in bins

     tmp.rotate(ph);
 		X+=tmp;
   }
 }//CuDFTW
 
-/*
+/**
   function TFReas: time-frequency reassignment
 
   In: data[Wid]: waveform data
       win[Wid+1], dwin[Wid+1], ddwin[Wid+1]: window function and its derivatives
       f, t: initial digital frequency and time

          dwdt=(xtt.y*x.x-xtt.x*x.y)/px-2*(xt.x*x.x+xt.y*x.y)*(xt.y*x.x-xt.x*x.y)/px/px;
   if (dtdt!=0) fslope=dwdt/dtdt/(2*M_PI);
   else fslope=0;
 } //TFReas*/
 
-/*
+/**
   function TFReas: sinusoid estimation using reassignment method
 
   In: data[Wid]: waveform data
       w[Wid+1], dw[Wid+1], ddw[Wid+1]: window function and its derivatives
       win[Wid]: window function used for estimating amplitude and phase by projection onto a chirp

 
     Further reading: Wen X. and M. Sandler, "Additive and multiplicative reestimation schemes
     for the sinusoid modeling of audio," in Proc. EUSIPCO'09, Glasgow, 2009.
 */
 
-/*
+/**
   function AdditiveUpdate: additive reestimation of time-varying sinusoid
 
   In: x[Count]: waveform data
       Wid, Offst: frame size and hop
       fs[Count], as[Count], phs[Count]: initial estimate of sinusoid parameters

     if (fs[i]<0 || fs[i]>0.5){}
 	}
 	delete[] y; delete[] lf; delete[] ref;
 }//AdditiveUpdate
 
-/*
+/**
   function AdditiveAnalyzer: sinusoid analyzer with one additive update
 
   In: x[Count]: waveform data
       Wid, Offst: frame size and hop size
       BasicAnalyzer: pointer to a sinusoid analyzer

 {
 	BasicAnalyzer(fs, as, phs, das, x, Count, Wid, Offst, ref, reserved, LogA);
 	AdditiveUpdate(fs, as, phs, das, x, Count, Wid, Offst, BasicAnalyzer, reserved, LogA);
 }//AdditiveAnalyzer
 
-/*
+/**
   function MultiplicativeUpdate: multiplicative reestimation of time-varying sinusoid
 
   In: x[Count]: waveform data
       Wid, Offst: frame size and hop
       fs[Count], as[Count], phs[Count]: initial estimate of sinusoid parameters

 	}
 
 	delete[] y; delete[] lf; delete[] lref;
 }//MultiplicativeUpdate
 
-/*
+/**
   function MultiplicativeAnalyzer: sinusoid analyzer with one multiplicative update
 
   In: x[Count]: waveform data
       Wid, Offst: frame size and hop size
       BasicAnalyzer: pointer to a sinusoid analyzer

 
     Further reading: Wen X. and M. Sandler, "Evaluating parameters of time-varying
     sinusoids by demodulation," in Proc. DAFx'08, Espoo, 2008.
 */
 
-/*
+/**
   function ReEstFreq: sinusoid reestimation by demodulating frequency.
 
   In: x[Wid+Offst*(FrCount-1)]: waveform data
       FrCount, Wid, Offst: frame count, frame size and hop size
       fbuf[FrCount], ns[FrCount]: initial frequency estiamtes and their timing

         fbuf[fr]=lf/Wid, abuf[fr]=la*2, pbuf[fr]=lp;
     }
   }
 }//ReEstFreq
 
-/*
+/**
   function ReEstFreq_2: sinusoid reestimation by demodulating frequency. This is that same as ReEstFreq(...)
   except that it calls Sinusoid(...) to synthesize the phase track used for demodulation and that it
   does not allow variable window sizes for estimating demodulated sinusoid.
 
   In: x[Wid+Offst*(FrCount-1)]: waveform data

     if (la*2>abuf[fr])
       fbuf[fr]=lf/Wid, abuf[fr]=la*2, pbuf[fr]=lp+centralph;
   }
 }//ReEstFreq_2
 
-/*
+/**
   function ReEstFreqAmp: sinusoid reestimation by demodulating frequency and amplitude.
 
   In: x[Wid+Offst*(FrCount-1)]: waveform data
       FrCount, Wid, Offst: frame count, frame size and hop size
       fbuf[FrCount], abuf[FrCount], ns[FrCount]: initial frequency and amplitude estiamtes and their

       if (la*2>abuf[fr]) fbuf[fr]=lf/Wid, abuf[fr]=la*2, pbuf[fr]=lp;
     }
   }
 }//ReEstFreqAmp
 
-/*
+/**
   function Reestimate2: iterative demodulation method for sinusoid parameter reestimation.
 
   In: x[(FrCount-1)*Offst+Wid]: waveform data
       FrCount, Wid, Offst: frame count, frame size and hop size
       win[Wid]: window function

 
     Further reading: Wen X. and M. Sandler, "Notes on model-based non-stationary sinusoid estimation methods
     using derivatives," in Proc. DAFx'09, Como, 2009.
 */
 
-/*
+/**
   function Derivative: derivative method for estimating amplitude derivative, frequency, and frequency derivative given
   signal and its derivatives.
 
   In: x[Wid], dx[Wid], ddx[Wid]: waveform and its derivatives
       win[Wid]: window function

   *a1=miu0;
 
   FreeFFTBuffer(fft);
 }//Derivative
 
-/*
+/**
   function Xkw: computes windowed spectrum of x and its derivatives up to order K at angular frequency omg,
   from x using window w and its derivatives.
 
   In: x[Wid]: waveform data
       w[K+1][Wid]: window functions and its derivatives up to order K

     cdouble *thisX=&X[k], *lastX=&X[k-1];
 			for (int kk=K-k; kk>=0; kk--) thisX[kk]=-lastX[kk+1]+cdouble(0, omg)*lastX[kk];
   }
 }//Xkw
 
-/*
+/**
   function Xkw: computes windowed spectrum of x and its derivatives up to order K at angular frequency
   omg, from x and its derivatives using window w.
 
   In: x[K+1][Wid]: waveform data and its derivatives up to order K.
       w[Wid]: window function

       X[k]+=tmp.rotate(-ph);
     }
   }
 }//Xkw
 
-/*
+/**
   function Derivative: derivative method for estimating the model log(s)=h[M]'r[M], by discarding extra
   equations
 
   In: s[Wid]: waveform data
       win[][Wid]: window function and its derivatives

 	delete[] ind;
 	DeAlloc2(Smkw);
 	DeAlloc2(A);
 }//Derivative*/
 
-/*
+/**
   function DerivativeLS: derivative method for estimating the model log(s)=h[M]'r[M], least-square
   implementation
 
   In: s[Wid]: waveform data
       win[][Wid]: window function and its derivatives

 	delete[] ind;
 	DeAlloc2(Smkw);
 	DeAlloc2(A);
 }//DerivativeLS
 
-/*
+/**
   function DerivativeLS: derivative method for estimating the model log(s)=h[M]'r[M] using Fr
   measurement points a quarter of Wid apart from each other, implemented by least-square.
 
   In: s[Wid+(Fr-1)*Wid/4]: waveform data
       win[][Wid]: window function and its derivatives

 
     Further reading: M. Abe and J. O. Smith III, ¡°AM/FM rate estimation for time-varying sinusoidal
     modeling,¡± in Proc. ICASSP'05, Philadelphia, 2005.
 */
 
-/*
+/**
   function RDFTW: windowed DTFT at frequency k bins
 
   In: data[Wid]: waveform data
       w[Wid]: window function
       k: frequency, in bins

     Xr+=tmp*cos(ph);
     Xi+=tmp*sin(ph); //*/
   }
 }//RDFTW
 
-/*
+/**
   function TFAS05: the Abe-Smith method 2005
 
   In: data[Wid]: waveform data
       w[Wid]: window function
       res: resolution of frequency for QIFFT

 				 yeses=y+zt[6]*alf*alf/p+zt[7]*log(1+beta_p*beta_p),
 				 pheses=ph+zt[8]*alf*alf*beta/p+zt[9]*atan(beta_p); //*/
 	f=feses/Wid, a=exp(yeses), ph=pheses, fslope=2*beta/2/M_PI, aesp=alf;
 }//TFAS05
 
-/*
+/**
   function TFAS05_enh: the Abe-Smith method 2005 enhanced by LSE amplitude and phase estimation
 
   In: data[Wid]: waveform data
       w[Wid]: window function
       res: resolution of frequency for QIFFT

 	double aesp, fslope;
 	TFAS05_enh(f, t, a, ph, aesp, fslope, Wid, data, w, res);
 }//TFAS05_enh
 
 //---------------------------------------------------------------------------
-/*
+/**
   function DerivativeLSv_AmpPh: estimate the constant-term in the local derivative method. This is used
   by the local derivative algorithm, whose implementation is found in the header file as templates.
 
   In: sv0: inner product <s, v0>, where s is the sinusoid being estimated.
       integr_h[M][Wid]: M vectors containing samples of the integral of basis functions h[M].

 
     Further reading: Wen X. and M. Sandler, "Spline exponential approximation of time-varying
     sinusoids," under review.
 */
 
-/*
+/**
   function setv: computes I test functions v[I] by modulation u[I] to frequency f
 
   In: u[I+1][Wid], du[I+1][Wid]: base-band test functions and their derivatives
       f: carrier frequency
   Out: v[I][Wid], dv[I][Wid]: test functions and their derivatives

     if (f>=fbin || I%2==1){v[I-1][c]=u[I-1][c]*rot; dv[I-1][c]=du[I-1][c]*rot+jomg*v[I-1][c];}
     else{v[I-1][c]=u[I][c]*rot; dv[I-1][c]=du[I][c]*rot+jomg*v[I-1][c];}
   }
 }//setv
 
-/*
+/**
   function setvhalf: computes I half-size test functions v[I] by modulation u[I] to frequency f.
 
   In: u[I][hWid*2], du[I][Wid*2]: base-band test functions and their derivatives
       f: carrier frequency
   Out: v[I][hWid], dv[hWid]: half-size test functions and their derivatives

   }
 }//setvhalf
 
 //#define ERROR_CHECK
 
-/*
+/**
   function DerivativePiecewise: Piecewise derivative algorithm. In this implementation of the piecewise
   method the test functions v are constructed from I "basic" (single-frame) test functions, each
   covering the same period of 2T, by shifting these I functions by steps of T. A total number of (L-1)I
   test functions are used.
 

   else LSLinear(L_1*I, N, aita, srv[0], sv[0]);
 
   delete mlist;
 }//DerivativePiecewise
 
-/*
+/**
   function DerivativePiecewise2: Piecewise derivative algorithm in which the real and imaginary parts of
   the exponent are modelled separately. In this implementation of the piecewise method the test
   functions v are constructed from I "basic" (single-frame) test functions, each covering the same
   period of 2T, by shifting these I functions by steps of T. A total number of (L-1)I test functions are
   used.

   }
   DeAlloc2(v); DeAlloc2(dv);
   return ene;
 }//testsv
 
-/*
+/**
   function DerivativePiecewise3: Piecewise derivative algorithm in which the log amplitude and frequeny
   are modeled separately as piecewise functions. In this implementation of the piecewise method the test
   functions v are constructed from I "basic" (single-frame) test functions, each covering the same
   period of 2T, by shifting these I functions by steps of T. A total number of (L-1)I test functions are
   used.

   delete mlist;
 }//DerivativePiecewise3
 
 //initialization routines for the piecewise derivative method
 
-/*
+/**
   function seth: set h[M] to a series of power functions.
 
   In: M, T.
   Out: h[M][T], where h[m] is power function of order m.
 

  {
    hm=h[m]; for (int t=0; t<T; t++) hm[t]=pow(t*1.0, m);
  }
 }//seth
 
-/*
+/**
   function setdh: set dh[M] to the derivative of a series of power functions.
 
   In: M, T.
   Out: dh[M][T], where dh[m] is derivative of the power function of order m.
 

  {
    dhm=dh[m]; for (int t=0; t<T; t++) dhm[t]=m*pow(t*1.0, m-1);
  }
 }//setdh
 
-/*
+/**
   function setdih: set dih[M] to the difference of the integral of a series of power functions.
 
   In: M, I
   Out: dih[M][I], where the accumulation of dih[m] is the integral of the power function of order m.
 

   {
     dihm=dih[m]; for (int t=0; t<T; t++)	dihm[t]=(pow(t+1.0, m+1)-pow(t*1.0, m+1))/(m+1);
   }
 }//setdih
 
-/*
+/**
   function sshLinear: sets M and h[M] for the linear spline model
 
   In: T
   Out: M=2, h[2][T] filled out for linear spline model.
 

 {
   M=2; Allocate2L(double, M, T, h, mlist);
   for (int t=0; t<T; t++) h[0][t]=1, h[1][t]=t;
 }//sshLinear
 
-/*
+/**
   function sdihLinear: sets dih[M] for the linear spline model. For testing only.
 
   In: T
   Out: dih[2][T] filled out for linear spline model.
 

 {
   Allocate2L(double, 2, T, dih, mlist);
   for (int t=0; t<T; t++)	dih[0][t]=1, dih[1][t]=t+0.5;
 }//sdihLinear
 
-/*
+/**
   function sshCubic: sets M and h[M] for cubic spline models.
 
   In: T
   Out: M=4 and h[M] filled out for cubic spline models, including cubic and cubic-Hermite.
 

 {
   M=4; Allocate2L(double, M, T, h, mlist);
   for (int t=0; t<T; t++) h[3][t]=t*t*t, h[2][t]=t*t, h[1][t]=t, h[0][t]=1;
 }//sshCubic
 
-/*
+/**
   function sdihCubic: sets dih[M] for cubic spline models.
 
   In: T
   Out: dih[4] filled out for cubic spline models.
 

   {
     dih[3][t]=t*(t*(t+1.5)+1)+0.25, dih[2][t]=t*(t+1)+1.0/3, dih[1][t]=t+0.5, dih[0][t]=1;
   }
 }//sdihCubic*/
 
-/*
+/**
   function ssALinearSpline: sets N and A[L] for the linear spline model
 
   In: L, M, T
   Out: N=L+1, A[L][M][N] filled out for the linear spline model
 

   Allocate3L(double, L, M, N, A, mlist);
   memset(A[0][0], 0, sizeof(double)*L*M*N);
   double iT=1.0/T; for (int l=0; l<L; l++) A[l][0][l]=1, A[l][1][l]=-iT, A[l][1][l+1]=iT;
 }//ssALinearSpline
 
-/*
+/**
   function ssLinearSpline: sets M, N, h and A for the linear spline model
 
   In: L, M, T
   Out: N and h[][] and A[][][] filled out for the linear spline model
 

 {
   seth(M, T, h, mlist);
   ssALinearSpline(L, T, M, N, A, mlist);
 }//ssLinearSpline
 
-/*
+/**
   function ssACubicHermite: sets N and A[L] for cubic Hermite spline model
 
   In: L, M, T
   Out: N=2(L+1), A[L][M][N] filled out for the cubic Hermite spline
 

     A[l][1][2*l+1]=1;
     A[l][0][2*l]=1;
   }
 }//ssACubicHermite
 
-/*
+/**
   function ssLinearSpline: sets M, N, h and A for the cubic Hermite spline model
 
   In: L, M, T
   Out: N and h[][] and A[][][] filled out for the cubic Hermite spline model
 

 {
   seth(M, T, h, mlist);
   ssACubicHermite(L, T, M, N, A, mlist);
 }//ssCubicHermite
 
-/*
+/**
   function ssACubicSpline: sets N and A[L] for cubic spline model
 
   In: L, M, T
       mode: boundary mode of cubic spline, 0=natural, 1=quadratic run-out, 2=cubic run-out
   Out: N=2(L+1), A[L][M][N] filled out for the cubic spline

     A[l][1][l]-=iT; A[l][1][l+1]+=iT; A[l][0][l]+=1;
   }
   DeAlloc2(ML); DeAlloc2(MR); DeAlloc2(M42);
 }//ssACubicSpline
 
-/*
+/**
   function ssLinearSpline: sets M, N, h and A for the cubic spline model
 
   In: L, M, T
   Out: N and h[][] and A[][][] filled out for the cubic spline model
 

 {
   seth(M, T, h, mlist);
   ssACubicSpline(L, T, M, N, A, mlist, mode);
 }//ssCubicSpline
 
-/*
+/**
   function setu: sets u[I+1] as base-band windowed Fourier atoms, whose frequencies come in the order of
   0, 1, -1, 2, -2, 3, -3, 4, etc, in bins.
 
   In: I, Wid: number and size of atoms to generate.
       WinOrder: order (=vanishing moment) of window function to use (2=Hann, 4=Hann^2, etc.)

     }
   }
   DeAlloc2(wins);
 }//setu
 
-/*
+/**
   function DerivativePiecewiseI: wrapper for DerivativePiecewise(), doing the initialization ,etc.
 
   In: L, T: number and length of pieces
       s[LT]: waveform signal
       ds[LT]: derivative of s[LT], used only when ERROR_CHECK is defined.

 
   DerivativePiecewise(N, aita, L, f, T, s, A, M, h, I, u, du, endmode, ds);
   delete mlist;
 }//DerivativePiecewiseI
 
-/*
+/**
   function DerivativePiecewiseII: wrapper for DerivativePiecewise2(), doing the initialization ,etc.
   This models the derivative of log ampltiude and frequency as separate piecewise polynomials, the first
   specified by SpecifyA, the second by SpecifyB.
 
   In: L, T: number and length of pieces

   DerivativePiecewise2(Np, p, Nq, q, L, f, T, s, A, B, M, h, I, u, du, endmode, ds);
 
   delete mlist;
 }//DerivativePiecewiseII
 
-/*
+/**
   function DerivativePiecewiseIII: wrapper for DerivativePiecewise3(), doing the initialization ,etc.
   Notice that this time the log amplitude, rather than its derivative, is modeled as a piecewise
   polynomial specified by SpecifyA.
 
   In: L, T: number and length of pieces

   DerivativePiecewise3(Np, p, Nq, q, L, f, T, s, A, B, M, h, I, u, du, endmode, ds, dh);
 
   delete mlist;
 }//DerivativePiecewiseIII
 
-/*
+/**
   function AmpPhCorrectionExpA: model-preserving amplitude and phase correction in piecewise derivative
   method.
 
   In: aita[N]: inital independent coefficients
       L, T: number and size of pieces

   }
   double result=s2ph[0];
   delete mlist;
   return result;
 }//AmpPhCorrectionExpA
+