comparison procedures.cpp @ 5:5f3c32dc6e17

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

 #include "procedures.h"
 #include "matrix.h"
 #include "opt.h"
 #include "sinest.h"
 
+/** \file procedures.h */
+
 //---------------------------------------------------------------------------
 //TGMM methods
 
 //method TGMM::TGMM: default constructor
 TGMM::TGMM()

 TTFSpans::~TTFSpans()
 {
   delete[] Items;
 }//~TTFSpans
 
-/*
+/**
   method Add: add a new span to the list
 
   In: ATFSpan: the new span to add
 */
 void TTFSpans::Add(TTFSpan& ATFSpan)

   }
   Items[Count]=ATFSpan;
   Count++;
 }//Add
 
-/*
+/**
   method Clear: discard the current content without freeing memory.
 */
 void TTFSpans::Clear()
 {
   Count=0;
 }//Clear
 
-/*
+/**
   method Delete: delete a span from current list
 
   In: Index: index to the span to delete
 */
 int TTFSpans::Delete(int Index)

 }//Delete
 
 //---------------------------------------------------------------------------
 //SpecTrack methods
 
-/*
+/**
   method TSpecTrack::Add: adds a SpecPeak to the track.
 
   In: APeak: the SpecPeak to add.
 */
 int TSpecTrack::Add(TSpecPeak& APeak)

     else if (f>fmax) fmax=f;
   }
   return ind;
 }//Add
 
-/*
+/**
   method TSpecTrack::TSpecTrack: creates a SpecTrack with an inital capacity.
 
   In: ACapacity: initial capacity, i.e. the number SpecPeak's to allocate storage space for.
 */
 TSpecTrack::TSpecTrack(int ACapacity)

 TSpecTrack::~TSpecTrack()
 {
   delete[] Peaks;
 }//TSpecTrack
 
-/*
+/**
   method InsertPeak: inserts a new SpecPeak into the track at a given index. Internal use only.
 
   In: APeak: the SpecPeak to insert.
       index: the position in the list to place the new SpecPeak. Original SpecPeak's at and after this
         position are shifted by 1 posiiton.

   memmove(&Peaks[index+1], &Peaks[index], sizeof(TSpecPeak)*(Count-index));
   Peaks[index]=APeak;
   Count++;
 }//InsertPeak
 
-/*
+/**
   method TSpecTrack::LocatePeak: looks for a SpecPeak in the track that has the same time (t) as APeak.
 
   In: APeak: a SpecPeak
 
   Returns the index in this track of the first SpecPeak that has the same time stamp as APeak. However,

   }
   return -a-2;
 }//LocatePeak
 
 //---------------------------------------------------------------------------
-/*
+/**
   function: ACPower: AC power
 
   In: data[Count]: a signal
 
   Returns the power of its AC content.

   power=(power-avg*avg)/Count;
   return power;
 }//ACPower
 
 //---------------------------------------------------------------------------
-/*
+/**
   function Add: vector addition
 
   In: dest[Count], source[Count]: two vectors
   Out: dest[Count]: their sum
 

 void Add(double* dest, double* source, int Count)
 {
   for (int i=0; i<Count; i++) *(dest++)+=*(source++);
 }//Add
 
-/*
+/**
   function Add: vector addition
 
   In: addend[count], adder[count]: two vectors
   Out: out[count]: their sum
 

   for (int i=0; i<count; i++) *(out++)=*(addend++)+*(adder++);
 }//Add
 
 //---------------------------------------------------------------------------
 
-/*
+/**
   function ApplyWindow: applies window function to signal buffer.
 
   In: Buffer[Size]: signal to be windowed
       Weight[Size]: the window
   Out: OutBuffer[Size]: windowed signal

 {
   for (int i=0; i<Size; i++) *(OutBuffer++)=*(Buffer++)**(Weights++);
 }//ApplyWindow
 
 //---------------------------------------------------------------------------
-/*
+/**
   function Avg: average
 
   In: data[Count]: a data array
 
   Returns the average of the array.

   avg/=Count;
   return avg;
 }//Avg
 
 //---------------------------------------------------------------------------
-/*
+/**
   function AvgFilter: get slow-varying wave trace by averaging
 
   In: data[Count]: input signal
       HWid: half the size of the averaging window
   Out: datout[Count]: the slow-varying part of data[].

     dataout[i]=sum/(2*(Count-i)-1);
   }
 }//AvgFilter
 
 //---------------------------------------------------------------------------
-/*
+/**
   function CalculateSpectrogram: computes the spectrogram of a signal
 
   In: data[Count]: the time-domain signal
       start, end: start and end points marking the section for which the spectrogram is to be computed
       Wid, Offst: frame size and hop size

   }
   FreeFFTBuffer(fft);
 }//CalculateSpectrogram
 
 //---------------------------------------------------------------------------
-/*
+/**
   function Conv: simple convolution
 
   In: in1[N1], in2[N2]: two sequences
   Out: out[N1+N2-1]: their convolution
 

     for (int n2=0; n2<N2; n2++)
       out[n1+n2]+=in1[n1]*in2[n2];
 }//Conv
 
 //---------------------------------------------------------------------------
-/*
+/**
   function Correlation: computes correlation coefficient of 2 vectors a & b, equals cos(aOb).
 
   In: a[Count], b[Count]: two vectors
 
   Returns their correlation coefficient.

   }
   return ab/sqrt(aa*bb);
 }//Correlation
 
 //---------------------------------------------------------------------------
-/*
+/**
   function DCAmplitude: DC amplitude
 
   In: data[Count]: a signal
 
   Returns its DC amplitude (=AC amplitude without DC removing)

   }
   power/=Count;
   return sqrt(2*power);
 }//DCAmplitude
 
-/*
+/**
   function DCPower: DC power
 
   In: data[Count]: a signal
 
   Returns its DC power.

   power/=Count;
   return power;
 }//DCPower
 
 //---------------------------------------------------------------------------
-/*
+/**
   DCT: discrete cosine transform, direct computation. For fast DCT, see fft.cpp.
 
   In: input[N]: a signal
   Out: output[N]: its DCT
 

     if (n==0) output[n]*=1.4142135623730950488016887242097/N;
     else output[n]*=2.0/N;
   }
 }//DCT
 
-/*
+/**
   function IDCT: inverse discrete cosine transform, direct computation. For fast IDCT, see fft.cpp.
 
   In: input[N]: a signal
   Out: output[N]: its IDCT
 

       output[k]+=input[n]*cos(n*Wk);
   }
 }//IDCT
 
 //---------------------------------------------------------------------------
-/*
+/**
   function DeDC: removes DC component of a signal
 
   In: data[Count]: the signal
       HWid: half of averaging window size
   Out: data[Count]: de-DC-ed signal

   for (int i=0; i<Count; i++)
     *(data++)-=*(data2++);
   delete[] data2;
 }//DeDC
 
-/*
+/**
   function DeDC_static: removes DC component statically
 
   In: data[Count]: the signal
   Out: data[Count]: DC-removed signal
 

   double avg=Avg(data, Count, 0);
   for (int i=0; i<Count; i++) *(data++)-=avg;
 }//DeDC_static
 
 //---------------------------------------------------------------------------
-/*
+/**
   function DoubleToInt: converts double-precision floating point array to integer array
 
   In: in[Count]: the double array
       BytesPerSample: bytes per sample of destination integers
   Out: out[Count]: the integer array

 {
   if (BytesPerSample==1){unsigned char* out8=(unsigned char*)out; for (int k=0; k<Count; k++) *(out8++)=*(in++)+128.5;}
   else {__int16* out16=(__int16*)out; for (int k=0; k<Count; k++) *(out16++)=floor(*(in++)+0.5);}
 }//DoubleToInt
 
-/*
+/**
   function DoubleToIntAdd: adds double-precision floating point array to integer array
 
   In: in[Count]: the double array
       out[Count]: the integer array
       BytesPerSample: bytes per sample of destination integers

     for (int k=0; k<Count; k++){*out16=*out16+floor(*in+0.5); out16++; in++;}
   }
 }//DoubleToIntAdd
 
 //---------------------------------------------------------------------------
-/*
+/**
   DPower: in-frame power variation
 
   In: data[Count]: a signal
 
   returns the different between AC powers of its first and second halves.

   double ene2=ACPower(&data[Count/2], Count/2, 0);
   return ene2-ene1;
 }//DPower
 
 //---------------------------------------------------------------------------
-/*
-  funciton Energy: energy
+/**
+  function Energy: energy
 
   In: data[Count]: a signal
 
   Returns its total energy
 */

   for (int i=0; i<Count; i++) result+=data[i]*data[i];
   return result;
 }//Energy
 
 //---------------------------------------------------------------------------
-/*
+/**
   function ExpOnsetFilter: onset filter with exponential impulse response h(t)=Aexp(-t/Tr)-Bexp(-t/Ta),
   A=1-exp(-1/Tr), B=1-exp(-1/Ta).
 
   In: data[Count]: signal to filter
       Tr, Ta: time constants of h(t)

     *(dataout++)=(A*FA-B*FB)*NormFactor;
   }
   return NormFactor;
 }//ExpOnsetFilter
 
-/*
+/**
   function ExpOnsetFilter_balanced: exponential onset filter without starting step response. It
   extends the input signal at the front end by bal*Ta samples by repeating the its value at 0, then
   applies the onset filter on the extended signal instead.
 
   In: data[Count]: signal to filter

   delete[] tmpdata;
   return result;
 }//ExpOnsetFilter_balanced
 
 //---------------------------------------------------------------------------
-/*
+/**
   function ExtractLinearComponent: Legendre linear component
 
   In: data[Count+1]: a signal
   Out: dataout[Count+1]: its Legendre linear component, optional.
 

     for (int n=0; n<=Count; n++) *(dataout++)=tmp*n;
   return tmp;
 }//ExtractLinearComponent
 
 //---------------------------------------------------------------------------
-/*
+/**
   function FFTConv: fast convolution of two series by FFT overlap-add. In an overlap-add scheme it is
   assumed that one of the convolvends is short compared to the other one, which can be potentially
   infinitely long. The long convolvend is devided into short segments, each of which is convolved with
   the short convolvend, the results of which are then assembled into the final result. The minimal delay
   from input to output is the amount of overlap, which is the size of the short convolvend minus 1.

   delete[] x2;
   delete[] tmp;
   delete[] hbitinv;
 }//FFTConv
 
-/*
+/**
   function FFTConv: the simplified version using two output buffers instead of three. This is almost
   equivalent to setting zero=-size2 in the three-output-buffer version (so that post_buffer no longer
   exists), except that this version requires size2 (renamed HWid) be a power of 2, and pre_buffer point
   to the END of the storage (so that pre_buffer=dest automatically connects the two buffers in a
   continuous memory block).

 
   FreeFFTBuffer(fft);
   delete[] x2; delete[] tmp; delete[] hbitinv;
 }//FFTConv
 
-/*
+/**
   function FFTConv: fast convolution of two series by FFT overlap-add. Same as the three-output-buffer
   version above but using integer output buffers as well as integer source1.
 
   In: source1[size1]: convolvend
       bps: bytes per sample of integer units in source1[].

   delete[] tmp;
   delete[] hbitinv;
 }//FFTConv
 
 //---------------------------------------------------------------------------
-/*
+/**
   function FFTFilter: FFT with cosine-window overlap-add: This FFT filter is not an actural LTI system,
   but an block processing with overlap-add. In this function the blocks are overlapped by 50% and summed
   up with Hann windowing.
 
   In: data[Count]: input data

   delete[] win;
   delete[] tmpdata;
   FreeFFTBuffer(ldata);
 }//FFTFilter
 
-/*
-  funtion FFTFilterOLA: FFTFilter with overlap-add support. This is a true LTI filter whose impulse
+/**
+  function FFTFilterOLA: FFTFilter with overlap-add support. This is a true LTI filter whose impulse
   response is constructed using IFFT. The filtering is implemented by fast convolution.
 
   In: data[Count]: input data
       Wid: FFT size
       On, Off: cut-off frequencies, in bins, of the filter

   CIFFTR(x, log2(Wid), w, spec);
   FFTConv(dataout, data, BytesPerSample, Count, spec, Wid, -Wid, pre_buffer, NULL);
   FreeFFTBuffer(spec);
 }//FFTFilterOLA
 
-/*
+/**
   function FFTFilterOLA: FFT filter with overlap-add support.
 
   In: data[Count]: input data
       amp[0:HWid]: amplitude response
       ph[0:HWid]: phase response, where ph[0]=ph[HWid]=0;

   CIFFTR(x, log2(Wid), w, spec);
   FFTConv(dataout, data, Count, spec, Wid, -Wid, pre_buffer, NULL);
   FreeFFTBuffer(spec);
 }//FFTFilterOLA
 
-/*
+/**
   function FFTMask: masks a band of a signal with noise
 
   In: data[Count]: input signal
       DigiOn, DigiOff: cut-off frequences of the band to mask
       maskcoef: masking noise amplifier. If set to 1 than the mask level is set to the highest signal

 
   return max;
 }//FFTMask
 
 //---------------------------------------------------------------------------
-/*
+/**
   function FindInc: find the element in ordered list data that is closest to value.
 
   In: data[Count]: a ordered list
       value: the value to locate in the list
 

   if (fabs(value-data[end-1])<fabs(value-data[end])) return end-1;
   else return end;
 }//FindInc
 
 //---------------------------------------------------------------------------
-/*
+/**
   function Gaussian: Gaussian function
 
   In: Vector[Dim]: a vector
       Mean[Dim]: mean of Gaussian function
       Dev[Fim]: diagonal autocorrelation matrix of Gaussian function

   return bmt;
 }//Gaussian
 
 
 //---------------------------------------------------------------------------
-/*
+/**
   function Hamming: calculates the amplitude spectrum of Hamming window at a given frequency
 
   In: f: frequency
       T: size of Hamming window
 

   c=c1+c2+c3;
   return abs(c);
 }//Hamming*/
 
 //---------------------------------------------------------------------------
-/*
+/**
   function HannSq: computes the square norm of Hann window spectrum (window-size-normalized)
 
   In: x: frequency, in bins
       N: size of Hann window
 

   im=0.25*sin(pif)*(-spmdivbyspmpf+spmdivbyspmmf);
 
   return (re*re+im*im)/(N*N);
 }//HannSq
 
-/*
+/**
   function Hann: computes the Hann window amplitude spectrum (window-size-normalized).
 
   In: x: frequency, in bins
       N: size of Hann window
 

   double result=0.5*spmdivbyspmf-0.25*(spmdivbyspmpf+spmdivbyspmmf);
 
   return result/N;
 }//HannC
 
-/*
+/**
   function HxPeak2: fine spectral peak detection. This does detection and high-precision LSE estimation
   in one go. However, since in practise most peaks are spurious, LSE estimation is not necessary on
   them. Accordingly, HxPeak2 is deprecated in favour of faster but coarser peak picking methods, such as
   QIFFT, which leaves fine estimation to a later stage of processing.
 

   if (allocvhps) vhps=(double*)realloc(vhps, sizeof(double)*count);
   return count;
 }//HxPeak2
 
 //---------------------------------------------------------------------------
-/*
+/**
   function InsertDec: inserts value into sorted decreasing list
 
   In: data[Count]: a sorted decreasing list.
       value: the value to be added
   Out: data[Count]: the list with $value inserted if the latter is larger than its last entry, in which

   memmove(&data[end+1], &data[end], sizeof(double)*(Count-end-1));
   data[end]=value;
   return end;
 }//InsertDec
 
-/*
+/**
   function InsertDec: inserts value and attached integer into sorted decreasing list
 
   In: data[Count]: a sorted decreasing list
       indices[Count]: a list of integers attached to entries of data[]
       value, index: the value to be added and its attached integer

   memmove(&indices[end+1], &indices[end], sizeof(int)*(Count-end-1));
   data[end]=value, indices[end]=index;
   return end;
 }//InsertDec
 
-/*
+/**
   InsertInc: inserts value into sorted increasing list.
 
   In: data[Count]: a sorted increasing list.
       Capacity: maximal size of data[]
       value: the value to be added

     data[PosToInsert]=value;
   }
   return PosToInsert;
 }//InsertInc
 
-/*
+/**
   function InsertInc: inserts value into sorted increasing list
 
   In: data[Count]: a sorted increasing list.
       value: the value to be added
       doinsert: specifies whether the actually insertion is to be performed

     data[end]=value;
   }
   return end;
 }//InsertInc
 
-/*
+/**
   function InsertInc: inserts value and attached integer into sorted increasing list
 
   In: data[Count]: a sorted increasing list
       indices[Count]: a list of integers attached to entries of data[]
       value, index: the value to be added and its attached integer

   memmove(&indices[end+1], &indices[end], sizeof(double)*(Count-end-1));
   data[end]=value, indices[end]=index;
   return end;
 }//InsertInc
 
-/*
+/**
   function InsertIncApp: inserts value into flexible-length sorted increasing list
 
   In: data[Count]: a sorted increasing list.
       value: the value to be added
   Out: data[Count+1]: the list with $value inserted.

 
   return end;
 }//InsertIncApp
 
 //---------------------------------------------------------------------------
-/*
+/**
   function InstantFreq; calculates instantaneous frequency from spectrum, evaluated at bin k
   
   In: x[hwid]: spectrum with scale 2hwid
       k: reference frequency, in bins
       mode: must be 1.

       break;
   }
   return result;
 }//InstantFreq
 
-/*
+/**
   function InstantFreq; calculates "frequency spectrum", a sequence of frequencies evaluated at DFT bins
   
   In: x[hwid]: spectrum with scale 2hwid
       mode: must be 1.
   Out: freqspec[hwid]: the frequency spectrum

   for (int i=0; i<hwid; i++)
     freqspec[i]=InstantFreq(i, hwid, x, mode);
 }//InstantFreq
 
 //---------------------------------------------------------------------------
-/*
+/**
   function IntToDouble: copy content of integer array to double array
 
   In: in: pointer to integer array
       BytesPerSample: number of bytes each integer takes
       Count: size of integer array, in integers

   if (BytesPerSample==1){unsigned char* in8=(unsigned char*)in; for (int k=0; k<Count; k++) *(out++)=*(in8++)-128.0;}
   else {__int16* in16=(__int16*)in; for (int k=0; k<Count; k++) *(out++)=*(in16++);}
 }//IntToDouble*/
 
 //---------------------------------------------------------------------------
-/*
+/**
   function IPHannC: inner product with Hann window spectrum
 
   In: x[N]: spectrum
       f: reference frequency
       K1, K2: spectral truncation bounds

   return abs(IPWindowC(f, x, N, M, c, iH2, K1, K2));
 }//IPHannC
 
 
 //---------------------------------------------------------------------------
-/*
+/**
   function lse: linear regression y=ax+b
 
   In: x[Count], y[Count]: input points
   Out: a, b: LSE estimation of coefficients in y=ax+b
 

   b=(sxx*sy-sx*sxy)/(Count*sxx-sx*sx);
   a=(sy-Count*b)/sx;
 }//lse
 
 //--------------------------------------------------------------------------
-/*
+/**
   memdoubleadd: vector addition
 
   In: dest[count], source[count]: addends
   Out: dest[count]: sum
 

 {
   for (int i=0; i<count; i++){*dest=*dest+*source; dest++; source++;}
 }//memdoubleadd
 
 //--------------------------------------------------------------------------
-/*
+/**
   function Mel: converts frequency in Hz to frequency in mel.
 
   In: f: frequency, in Hz
 
   Returns the frequency measured on mel scale.

 double Mel(double f)
 {
   return 1127.01048*log(1+f/700);
 }//Mel
 
-/*
+/**
   function InvMel: converts frequency in mel to frequency in Hz.
 
   In: f: frequency, in mel.
 
   Returns the frequency in Hz.

 double InvMel(double mel)
 {
   return 700*(exp(mel/1127.01048)-1);
 }//InvMel
 
-/*
+/**
   function MFCC: calculates MFCC.
 
   In: Data[FrameWidth]: data
       NumBands: number of frequency bands on mel scale
       Bands[3*NumBands]: mel frequency bands, comes as $NumBands triples, each containing the lower,

       tmp+=Amps[j]*cos(M_PI*(i+0.5)*(j+0.5)/NumBands);
     C[i]=tmp;
   }
 }//MFCC
 
-/*
+/**
   function MFCCPrepareBands: returns a array of OVERLAPPING bands given in triples, whose 1st and 3rd
   entries are the start and end of a band, in bins, and the 2nd is a middle frequency.
 
   In: SamplesPerSec: sampling rate
       NumberOfBins: FFT size

   }
   return Bands;
 }//MFCCPrepareBands
 
 //---------------------------------------------------------------------------
-/*
+/**
   function Multi: vector-constant multiplication
 
   In: data[count]: a vector
       mul: a constant
   Out: data[count]: their product

 void Multi(double* data, int count, double mul)
 {
   for (int i=0; i<count; i++){*data=*data*mul; data++;}
 }//Multi
 
-/*
+/**
   function Multi: vector-constant multiplication
 
   In: in[count]: a vector
       mul: a constant
   Out: out[count]: their product

 void Multi(double* out, double* in, int count, double mul)
 {
   for (int i=0; i<count; i++) *(out++)=*(in++)*mul;
 }//Multi
 
-/*
+/**
   function Multi: vector-constant multiply-addition
 
   In: in[count], adder[count]: vectors
       mul: a constant
   Out: out[count]: in[]+adder[]*mul

 {
   for (int i=0; i<count; i++) *(out++)=*(in++)+*(adder++)*mul;
 }//MultiAdd
 
 //---------------------------------------------------------------------------
-/*
+/**
   function NearestPeak: finds a peak value in an array that is nearest to a given index
 
   In: data[count]: an array
       anindex: an index
 

   else if (data[upind]<data[downind]) return downind;
   else return upind;
 }//NearestPeak
 
 //---------------------------------------------------------------------------
-/*
+/**
   function NegativeExp: fits the curve y=1-x^lmd.
 
   In: x[Count], y[Count]: sample points to fit, x[0]=0, x[Count-1]=1, y[0]=1, y[Count-1]=0
   Out: lmd: coefficient of y=1-x^lmd.
 

   while (e && iter<=maxiter && (!laste || fabs(laste-e)/e>eps));
   return sqrt(e/Count);
 }//NegativeExp
 
 //---------------------------------------------------------------------------
-/*
+/**
   function: NL: noise level, calculated on 5% of total frames with least energy
 
   In: data[Count]:
       Wid: window width for power level estimation
 

   result=sqrt(result);
   return result;
 }//NL
 
 //---------------------------------------------------------------------------
-/*
+/**
   function Normalize: normalizes data to [-Maxi, Maxi], without zero shift
 
   In: data[Count]: data to be normalized
       Maxi: destination maximal absolute value
   Out: data[Count]: normalized data

     for (int i=0; i<Count; i++) *(data++)*=Maxi;
   }
   return max;
 }//Normalize
 
-/*
+/**
   function Normalize2: normalizes data to a specified Euclidian norm
 
   In: data[Count]: data to normalize
       Norm: destination Euclidian norm
   Out: data[Count]: normalized data.

   if (mul!=0) for (int i=0; i<Count; i++) data[i]/=mul;
   return norm;
 }//Normalize2
 
 //---------------------------------------------------------------------------
-/*
+/**
   function PhaseSpan: computes the unwrapped phase variation across the Nyquist range
 
   In: data[Count]: time-domain data
       aparams: a fftparams structure
 

   FreeFFTBuffer(Amp);
   return result;
 }//PhaseSpan
 
 //---------------------------------------------------------------------------
-/*
+/**
   function PolyFit: least square polynomial fitting y=sum(i){a[i]*x^i}
 
   In: x[N], y[N]: sample points
       P: order of polynomial, P<N
   Out: a[P+1]: coefficients of polynomial

   GESCP(P+1, a, aa, bb);
   DeAlloc2(aa); delete[] bb; delete[] s; delete[] r;
 }//PolyFit
 
 //---------------------------------------------------------------------------
-/*
+/**
   function Pow: vector power function
 
   In: data[Count]: a vector
       ex: expontnet
   Out: data[Count]: point-wise $ex-th power of data[]

   for (int i=0; i<Count; i++)
     data[i]=pow(data[i], ex);
 }//Power
 
 //---------------------------------------------------------------------------
-/*
+/**
   Rectify: semi-wave rectification
 
   In: data[Count]: data to rectify
       th: rectification threshold: values below th are rectified to th
   Out: data[Count]: recitified data

   }
   return Result;
 }//Rectify
 
 //---------------------------------------------------------------------------
-/*
+/**
   function Res: minimum absolute residue.
 
   In: in: a number
       mod: modulus
 

   if (in<-mod/2) in+=mod;
   return in;
 }//Res
 
 //---------------------------------------------------------------------------
-/*
+/**
   function Romberg: Romberg algorithm for numerical integration
 
   In: f: function to integrate
       params: extra argument for calling f
       a, b: integration boundaries

   delete[] r1;
   delete[] r2;
   return h;
 }//Romberg
 
-/*
+/**
   function Romberg: Romberg algorithm for numerical integration, may return before specified depth on
   convergence.
 
   In: f: function to integrate
       params: extra argument for calling f

 
 //---------------------------------------------------------------------------
 //analog and digital sinc functions
 
 //sinca(0)=1, sincd(0)=N, sinca(1)=sincd(1)=0.
-/*
+/**
   function sinca: analog sinc function.
 
   In: x: frequency
 
   Returns sinc(x)=sin(pi*x)/(pi*x), sinca(0)=1, sinca(1)=0

 {
   if (x==0) return 1;
   return sin(M_PI*x)/(M_PI*x);
 }//sinca
 
-/*
+/**
   function sincd_unn: unnormalized discrete sinc function
 
   In: x: frequency
       N: scale (window size, DFT size)
 

   if (x==0) return N;
   return sin(M_PI*x)/sin(M_PI*x/N);
 }//sincd
 
 //---------------------------------------------------------------------------
-/*
+/**
   SmoothPhase: phase unwrapping on module mpi*PI, 2PI by default
 
   In: Arg[Count]: phase angles to unwrap
       mpi: unwrapping modulus, in pi's
   Out: Arg[Count]: unwrapped phase

 }//SmoothPhase
 
 //---------------------------------------------------------------------------
 //the stiff string partial frequency model f[m]=mf[1]*sqrt(1+B(m*m-1)).
 
-/*
+/**
   StiffB: computes stiffness coefficient from fundamental and another partial frequency based on the
   stiff string partial frequency model f[m]=mf[1]*sqrt(1+B(m*m-1)).
 
   In: f0: fundamental frequency
       m: 1-based partial index

   double f2=fm/m/f0;
   return (f2*f2-1)/(m*m-1);
 }//StiffB
 
 //StiffF: partial frequency of a stiff string
-/*
+/**
   StiffFm: computes a partial frequency from fundamental frequency and partial index based on the stiff
   string partial frequency model f[m]=mf[1]*sqrt(1+B(m*m-1)).
 
   In: f0: fundamental frequency
       m: 1-based partial index

 double StiffFm(double f0, int m, double B)
 {
   return m*f0*sqrt(1+B*(m*m-1));
 }//StiffFm
 
-/*
+/**
   StiffF0: computes fundamental frequency from another partial frequency and stiffness coefficient based
   on the stiff string partial frequency model f[m]=mf[1]*sqrt(1+B(m*m-1)).
 
   In: m: 1-based partial index
       fm: frequency of partial m

 double StiffF0(double fm, int m, double B)
 {
   return fm/m/sqrt(1+B*(m*m-1));
 }//StiffF0
 
-/*
+/**
   StiffM: computes 1-based partial index from partial frequency, fundamental frequency and stiffness
   coefficient based on the stiff string partial frequency model f[m]=mf[1]*sqrt(1+B(m*m-1)).
 
   In: f0: fundamental freqency
       fm: a partial frequency

   else
     return sqrt((b1+sqrt(delta))/2/B);
 }//StiffMd
 
 //---------------------------------------------------------------------------
-/*
+/**
   TFFilter: time-frequency filtering with Hann-windowed overlap-add.
 
   In: data[Count]: input data
       Spans: time-frequency spans
       wt, windp: type and extra parameter of FFT window