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 <math.h>
Chris@2: #include <string.h>
Chris@3: #include <stddef.h>
xue@1: #include "procedures.h"
xue@1: #include "matrix.h"
xue@1: #include "opt.h"
Chris@2: #include "sinest.h"
xue@1: 
Chris@5: /** \file procedures.h */
Chris@5: 
xue@1: //---------------------------------------------------------------------------
xue@1: //TGMM methods
xue@1: 
xue@1: //method TGMM::TGMM: default constructor
xue@1: TGMM::TGMM()
xue@1: {
xue@1:   p=0, m=dev=0;
xue@1: }//TGMM
xue@1: 
xue@1: //method GMM:~TGMM: default destructor
xue@1: TGMM::~TGMM()
xue@1: {
xue@1:   ReleaseGMM(p, m, dev)
xue@1: };
xue@1: 
xue@1: //---------------------------------------------------------------------------
xue@1: //TFSpans methods
xue@1: 
xue@1: //method TTFSpans: default constructor
xue@1: TTFSpans::TTFSpans()
xue@1: {
xue@1:   Count=0;
xue@1:   Capacity=100;
xue@1:   Items=new TTFSpan[Capacity];
xue@1: }//TTFSpans
xue@1: 
xue@1: //method ~TTFSpans: default destructor
xue@1: TTFSpans::~TTFSpans()
xue@1: {
xue@1:   delete[] Items;
xue@1: }//~TTFSpans
xue@1: 
Chris@5: /**
xue@1:   method Add: add a new span to the list
xue@1: 
xue@1:   In: ATFSpan: the new span to add
xue@1: */
xue@1: void TTFSpans::Add(TTFSpan& ATFSpan)
xue@1: {
xue@1:   if (Count==Capacity)
xue@1:   {
xue@1:     int OldCapacity=Capacity;
xue@1:     Capacity+=50;
xue@1:     TTFSpan* NewItems=new TTFSpan[Capacity];
xue@1:     memcpy(NewItems, Items, sizeof(TTFSpan)*OldCapacity);
xue@1:     delete[] Items;
xue@1:     Items=NewItems;
xue@1:   }
xue@1:   Items[Count]=ATFSpan;
xue@1:   Count++;
xue@1: }//Add
xue@1: 
Chris@5: /**
xue@1:   method Clear: discard the current content without freeing memory.
xue@1: */
xue@1: void TTFSpans::Clear()
xue@1: {
xue@1:   Count=0;
xue@1: }//Clear
xue@1: 
Chris@5: /**
xue@1:   method Delete: delete a span from current list
xue@1: 
xue@1:   In: Index: index to the span to delete
xue@1: */
xue@1: int TTFSpans::Delete(int Index)
xue@1: {
xue@1:   if (Index<0 || Index>=Count)
xue@1:     return 0;
xue@1:   memmove(&Items[Index], &Items[Index+1], sizeof(TTFSpan)*(Count-1-Index));
xue@1:   Count--;
xue@1:   return 1;
xue@1: }//Delete
xue@1: 
xue@1: //---------------------------------------------------------------------------
xue@1: //SpecTrack methods
xue@1: 
Chris@5: /**
xue@1:   method TSpecTrack::Add: adds a SpecPeak to the track.
xue@1: 
xue@1:   In: APeak: the SpecPeak to add.
xue@1: */
xue@1: int TSpecTrack::Add(TSpecPeak& APeak)
xue@1: {
xue@1:   if (Count>=Capacity)
xue@1:   {
xue@1:     Peaks=(TSpecPeak*)realloc(Peaks, sizeof(TSpecPeak)*(Capacity*2));
xue@1:     Capacity*=2;
xue@1:   }
xue@1:   int ind=LocatePeak(APeak);
xue@1:   if (ind<0)
xue@1:   {
xue@1:     InsertPeak(APeak, -ind-1);
xue@1:     ind=-ind-1;
xue@1:   }
xue@1: 
xue@1:   int t=APeak.t;
xue@1:   double f=APeak.f;
xue@1:   if (Count==1) t1=t2=t, fmin=fmax=f;
xue@1:   else
xue@1:   {
xue@1:     if (t<t1) t1=t;
xue@1:     else if (t>t2) t2=t;
xue@1:     if (f<fmin) fmin=f;
xue@1:     else if (f>fmax) fmax=f;
xue@1:   }
xue@1:   return ind;
xue@1: }//Add
xue@1: 
Chris@5: /**
xue@1:   method TSpecTrack::TSpecTrack: creates a SpecTrack with an inital capacity.
xue@1: 
xue@1:   In: ACapacity: initial capacity, i.e. the number SpecPeak's to allocate storage space for.
xue@1: */
xue@1: TSpecTrack::TSpecTrack(int ACapacity)
xue@1: {
xue@1:   Count=0;
xue@1:   Capacity=ACapacity;
xue@1:   Peaks=new TSpecPeak[Capacity];
xue@1: }//TSpecTrack
xue@1: 
xue@1: //method TSpecTrack::~TSpecTrack: default destructor.
xue@1: TSpecTrack::~TSpecTrack()
xue@1: {
xue@1:   delete[] Peaks;
xue@1: }//TSpecTrack
xue@1: 
Chris@5: /**
xue@1:   method InsertPeak: inserts a new SpecPeak into the track at a given index. Internal use only.
xue@1: 
xue@1:   In: APeak: the SpecPeak to insert.
xue@1:       index: the position in the list to place the new SpecPeak. Original SpecPeak's at and after this
xue@1:         position are shifted by 1 posiiton.
xue@1: */
xue@1: void TSpecTrack::InsertPeak(TSpecPeak& APeak, int index)
xue@1: {
xue@1:   memmove(&Peaks[index+1], &Peaks[index], sizeof(TSpecPeak)*(Count-index));
xue@1:   Peaks[index]=APeak;
xue@1:   Count++;
xue@1: }//InsertPeak
xue@1: 
Chris@5: /**
xue@1:   method TSpecTrack::LocatePeak: looks for a SpecPeak in the track that has the same time (t) as APeak.
xue@1: 
xue@1:   In: APeak: a SpecPeak
xue@1: 
xue@1:   Returns the index in this track of the first SpecPeak that has the same time stamp as APeak. However,
xue@1:   if there is no peak with that time stamp, the method returns -1 if APeaks comes before the first
xue@1:   SpecPeak, -2 if between 1st and 2nd SpecPeak's, -3 if between 2nd and 3rd SpecPeak's, etc.
xue@1: */
xue@1: int TSpecTrack::LocatePeak(TSpecPeak& APeak)
xue@1: {
xue@1:   if (APeak.t<Peaks[0].t) return -1;
xue@1:   if (APeak.t>Peaks[Count-1].t) return -Count-1;
xue@1: 
xue@1:   if (APeak.t==Peaks[0].t) return 0;
xue@1:   else if (APeak.t==Peaks[Count-1].t) return Count-1;
xue@1: 
xue@1:   int a=0, b=Count-1, c=(a+b)/2;
xue@1:   while (a<c)
xue@1:   {
xue@1:     if (APeak.t==Peaks[c].t) return c;
xue@1:     else if (APeak.t<Peaks[c].t) {b=c; c=(a+b)/2;}
xue@1:     else {a=c; c=(a+b)/2;}
xue@1:   }
xue@1:   return -a-2;
xue@1: }//LocatePeak
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function: ACPower: AC power
xue@1: 
xue@1:   In: data[Count]: a signal
xue@1: 
xue@1:   Returns the power of its AC content.
xue@1: */
xue@1: double ACPower(double* data, int Count, void*)
xue@1: {
xue@1:   if (Count<=0) return 0;
xue@1:   double power=0, avg=0, tmp;
xue@1:   for (int i=0; i<Count; i++)
xue@1:   {
xue@1:     tmp=*(data++);
xue@1:     power+=tmp*tmp;
xue@1:     avg+=tmp;
xue@1:   }
xue@8:   power=(power-avg*avg/Count)/Count;
xue@1:   return power;
xue@1: }//ACPower
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function Add: vector addition
xue@1: 
xue@1:   In: dest[Count], source[Count]: two vectors
xue@1:   Out: dest[Count]: their sum
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void Add(double* dest, double* source, int Count)
xue@1: {
xue@1:   for (int i=0; i<Count; i++) *(dest++)+=*(source++);
xue@1: }//Add
xue@1: 
Chris@5: /**
xue@1:   function Add: vector addition
xue@1: 
xue@1:   In: addend[count], adder[count]: two vectors
xue@1:   Out: out[count]: their sum
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void Add(double* out, double* addend, double* adder, int count)
xue@1: {
xue@1:   for (int i=0; i<count; i++) *(out++)=*(addend++)+*(adder++);
xue@1: }//Add
xue@1: 
xue@1: //---------------------------------------------------------------------------
xue@1: 
Chris@5: /**
xue@1:   function ApplyWindow: applies window function to signal buffer.
xue@1: 
xue@1:   In: Buffer[Size]: signal to be windowed
xue@1:       Weight[Size]: the window
xue@1:   Out: OutBuffer[Size]: windowed signal
xue@1: 
xue@1:   No return value;
xue@1: */
xue@1: void ApplyWindow(double* OutBuffer, double* Buffer, double* Weights, int Size)
xue@1: {
xue@1:   for (int i=0; i<Size; i++) *(OutBuffer++)=*(Buffer++)**(Weights++);
xue@1: }//ApplyWindow
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function Avg: average
xue@1: 
xue@1:   In: data[Count]: a data array
xue@1: 
xue@1:   Returns the average of the array.
xue@1: */
xue@1: double Avg(double* data, int Count, void*)
xue@1: {
xue@1:   if (Count<=0) return 0;
xue@1:   double avg=0;
xue@1:   for (int i=0; i<Count; i++) avg+=*(data++);
xue@1:   avg/=Count;
xue@1:   return avg;
xue@1: }//Avg
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function AvgFilter: get slow-varying wave trace by averaging
xue@1: 
xue@1:   In: data[Count]: input signal
xue@1:       HWid: half the size of the averaging window
xue@1:   Out: datout[Count]: the slow-varying part of data[].
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void AvgFilter(double* dataout, double* data, int Count, int HWid)
xue@1: {
xue@1:   double sum=0;
xue@1: 
xue@1:   dataout[0]=data[0];
xue@1: 
xue@1:   for (int i=1; i<=HWid; i++)
xue@1:   {
xue@1:     sum+=data[2*i-1]+data[2*i];
xue@1:     dataout[i]=sum/(2*i+1);
xue@1:   }
xue@1: 
xue@1:   for (int i=HWid+1; i<Count-HWid; i++)
xue@1:   {
xue@1:     sum=sum+data[i+HWid]-data[i-HWid-1];
xue@1:     dataout[i]=sum/(2*HWid+1);
xue@1:   }
xue@1: 
xue@1:   for (int i=Count-HWid; i<Count; i++)
xue@1:   {
xue@1:     sum=sum-data[2*i-Count-1]-data[2*i-Count];
xue@1:     dataout[i]=sum/(2*(Count-i)-1);
xue@1:   }
xue@1: }//AvgFilter
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function CalculateSpectrogram: computes the spectrogram of a signal
xue@1: 
xue@1:   In: data[Count]: the time-domain signal
xue@1:       start, end: start and end points marking the section for which the spectrogram is to be computed
xue@1:       Wid, Offst: frame size and hop size
xue@1:       Window: window function
xue@1:       amp: a pre-amplifier
xue@1:       half: specifies if the spectral values at Wid/2 are to be retried
xue@1:   Out: Spec[][Wid/2] or Spec[][Wid/2+1]: amplitude spectrogram
xue@1:        ph[][][Wid/2] or Ph[][Wid/2+1]: phase spectrogram
xue@1: 
xue@1:   No return value. The caller is repsonse to arrange storage spance of output buffers.
xue@1: */
xue@1: void CalculateSpectrogram(double* data, int Count, int start, int end, int Wid, int Offst, double* Window, double** Spec, double** Ph, double amp, bool half)
xue@1: {
xue@1:   AllocateFFTBuffer(Wid, fft, w, x);
xue@1: 
xue@1:   int Fr=(end-start-Wid)/Offst+1;
xue@1: 
xue@1:   for (int i=0; i<Fr; i++)
xue@1:   {
xue@11:     RFFTCW(&data[i*Offst+start], Window, 0, 0, Log2(Wid), w, x);
xue@1: 
xue@1:     if (Spec)
xue@1:     {
xue@1:       for (int j=0; j<Wid/2; j++)
xue@1:         Spec[i][j]=sqrt(x[j].x*x[j].x+x[j].y*x[j].y)*amp;
xue@1:       if (half)
xue@1:         Spec[i][Wid/2]=sqrt(x[Wid/2].x*x[Wid/2].x+x[Wid/2].y*x[Wid/2].y)*amp;
xue@1:     }
xue@1:     if (Ph)
xue@1:     {
xue@1:       for (int j=0; j<=Wid/2; j++)
xue@1:         Ph[i][j]=Atan2(x[j].y, x[j].x);
xue@1:       if (half)
xue@1:         Ph[i][Wid/2]=Atan2(x[Wid/2].y, x[Wid/2].x);
xue@1:     }
xue@1:   }
xue@1:   FreeFFTBuffer(fft);
xue@1: }//CalculateSpectrogram
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function Conv: simple convolution
xue@1: 
xue@1:   In: in1[N1], in2[N2]: two sequences
xue@1:   Out: out[N1+N2-1]: their convolution
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void Conv(double* out, int N1, double* in1, int N2, double* in2)
xue@1: {
xue@1:   int N=N1+N1-1;
xue@1:   memset(out, 0, sizeof(double)*N);
xue@1:   for (int n1=0; n1<N1; n1++)
xue@1:     for (int n2=0; n2<N2; n2++)
xue@1:       out[n1+n2]+=in1[n1]*in2[n2];
xue@1: }//Conv
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function Correlation: computes correlation coefficient of 2 vectors a & b, equals cos(aOb).
xue@1: 
xue@1:   In: a[Count], b[Count]: two vectors
xue@1: 
xue@1:   Returns their correlation coefficient.
xue@1: */
xue@1: double Correlation(double* a, double* b, int Count)
xue@1: {
xue@1:   double aa=0, bb=0, ab=0;
xue@1:   for (int i=0; i<Count; i++)
xue@1:   {
xue@1:     aa+=*a**a;
xue@1:     bb+=*b**b;
xue@1:     ab+=*(a++)**(b++);
xue@1:   }
xue@1:   return ab/sqrt(aa*bb);
xue@1: }//Correlation
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function DCAmplitude: DC amplitude
xue@1: 
xue@1:   In: data[Count]: a signal
xue@1: 
xue@1:   Returns its DC amplitude (=AC amplitude without DC removing)
xue@1: */
xue@1: double DCAmplitude(double* data, int Count, void*)
xue@1: {
xue@1:   if (Count<=0) return 0;
xue@1:   double power=0, tmp;
xue@1:   for (int i=0; i<Count; i++)
xue@1:   {
xue@1:     tmp=*(data++);
xue@1:     power+=tmp*tmp;
xue@1:   }
xue@1:   power/=Count;
xue@1:   return sqrt(2*power);
xue@1: }//DCAmplitude
xue@1: 
Chris@5: /**
xue@1:   function DCPower: DC power
xue@1: 
xue@1:   In: data[Count]: a signal
xue@1: 
xue@1:   Returns its DC power.
xue@1: */
xue@1: double DCPower(double* data, int Count, void*)
xue@1: {
xue@1:   if (Count<=0) return 0;
xue@1:   double power=0, tmp;
xue@1:   for (int i=0; i<Count; i++)
xue@1:   {
xue@1:     tmp=*(data++);
xue@1:     power+=tmp*tmp;
xue@1:   }
xue@1:   power/=Count;
xue@1:   return power;
xue@1: }//DCPower
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   DCT: discrete cosine transform, direct computation. For fast DCT, see fft.cpp.
xue@1: 
xue@1:   In: input[N]: a signal
xue@1:   Out: output[N]: its DCT
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void DCT( double* output, double* input, int N)
xue@1: {
xue@1:   double Wn;
xue@1: 
xue@1:   for (int n=0; n<N; n++)
xue@1:   {
xue@1:     output[n]=0;
xue@1:     Wn=n*M_PI/2/N;
xue@1:     for (int k=0; k<N; k++)
xue@1:       output[n]+=input[k]*cos((2*k+1)*Wn);
xue@1:     if (n==0) output[n]*=1.4142135623730950488016887242097/N;
xue@1:     else output[n]*=2.0/N;
xue@1:   }
xue@1: }//DCT
xue@1: 
Chris@5: /**
xue@1:   function IDCT: inverse discrete cosine transform, direct computation. For fast IDCT, see fft.cpp.
xue@1: 
xue@1:   In: input[N]: a signal
xue@1:   Out: output[N]: its IDCT
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void IDCT(double* output, double* input, int N)
xue@1: {
xue@1:   for (int k=0; k<N; k++)
xue@1:   {
xue@1:     double Wk=(2*k+1)*M_PI/2/N;
xue@1:     output[k]=input[0]/1.4142135623730950488016887242097;
xue@1:     for (int n=1; n<N; n++)
xue@1:       output[k]+=input[n]*cos(n*Wk);
xue@1:   }
xue@1: }//IDCT
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function DeDC: removes DC component of a signal
xue@1: 
xue@1:   In: data[Count]: the signal
xue@1:       HWid: half of averaging window size
xue@1:   Out: data[Count]: de-DC-ed signal
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void DeDC(double* data, int Count, int HWid)
xue@1: {
xue@1:   double* data2=new double[Count];
xue@1:   AvgFilter(data2, data, Count, HWid);
xue@1:   for (int i=0; i<Count; i++)
xue@1:     *(data++)-=*(data2++);
xue@1:   delete[] data2;
xue@1: }//DeDC
xue@1: 
Chris@5: /**
xue@1:   function DeDC_static: removes DC component statically
xue@1: 
xue@1:   In: data[Count]: the signal
xue@1:   Out: data[Count]: DC-removed signal
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void DeDC_static(double* data, int Count)
xue@1: {
xue@1:   double avg=Avg(data, Count, 0);
xue@1:   for (int i=0; i<Count; i++) *(data++)-=avg;
xue@1: }//DeDC_static
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function DoubleToInt: converts double-precision floating point array to integer array
xue@1: 
xue@1:   In: in[Count]: the double array
xue@1:       BytesPerSample: bytes per sample of destination integers
xue@1:   Out: out[Count]: the integer array
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void DoubleToInt(void* out, int BytesPerSample, double* in, int Count)
xue@1: {
xue@1:   if (BytesPerSample==1){unsigned char* out8=(unsigned char*)out; for (int k=0; k<Count; k++) *(out8++)=*(in++)+128.5;}
xue@1:   else {__int16* out16=(__int16*)out; for (int k=0; k<Count; k++) *(out16++)=floor(*(in++)+0.5);}
xue@1: }//DoubleToInt
xue@1: 
Chris@5: /**
xue@1:   function DoubleToIntAdd: adds double-precision floating point array to integer array
xue@1: 
xue@1:   In: in[Count]: the double array
xue@1:       out[Count]: the integer array
xue@1:       BytesPerSample: bytes per sample of destination integers
xue@1:   Out: out[Count]: the sum of the two arrays
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void DoubleToIntAdd(void* out, int BytesPerSample, double* in, int Count)
xue@1: {
xue@1:   if (BytesPerSample==1)
xue@1:   {
xue@1:     unsigned char* out8=(unsigned char*)out;
xue@1:     for (int k=0; k<Count; k++){*out8=*out8+*in+128.5; out8++; in++;}
xue@1:   }
xue@1:   else
xue@1:   {
xue@1:     __int16* out16=(__int16*)out;
xue@1:     for (int k=0; k<Count; k++){*out16=*out16+floor(*in+0.5); out16++; in++;}
xue@1:   }
xue@1: }//DoubleToIntAdd
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   DPower: in-frame power variation
xue@1: 
xue@1:   In: data[Count]: a signal
xue@1: 
xue@1:   returns the different between AC powers of its first and second halves.
xue@1: */
xue@1: double DPower(double* data, int Count, void*)
xue@1: {
xue@1:   double ene1=ACPower(data, Count/2, 0);
xue@1:   double ene2=ACPower(&data[Count/2], Count/2, 0);
xue@1:   return ene2-ene1;
xue@1: }//DPower
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
Chris@5:   function Energy: energy
xue@1: 
xue@1:   In: data[Count]: a signal
xue@1: 
xue@1:   Returns its total energy
xue@1: */
xue@1: double Energy(double* data, int Count)
xue@1: {
xue@1:   double result=0;
xue@1:   for (int i=0; i<Count; i++) result+=data[i]*data[i];
xue@1:   return result;
xue@1: }//Energy
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function ExpOnsetFilter: onset filter with exponential impulse response h(t)=Aexp(-t/Tr)-Bexp(-t/Ta),
xue@1:   A=1-exp(-1/Tr), B=1-exp(-1/Ta).
xue@1: 
xue@1:   In: data[Count]: signal to filter
xue@1:       Tr, Ta: time constants of h(t)
xue@1:   Out: dataout[Count]: filtered signal, normalized by multiplying a factor.
xue@1: 
xue@1:   Returns the normalization factor. Identical data and dataout is allowed.
xue@1: */
xue@1: double ExpOnsetFilter(double* dataout, double* data, int Count, double Tr, double Ta)
xue@1: {
xue@1:   double FA=0, FB=0;
xue@1:   double EA=exp(-1.0/Tr), EB=exp(-1.0/Ta);
xue@1:   double A=1-EA, B=1-EB;
xue@1:   double NormFactor=1/sqrt((1-EA)*(1-EA)/(1-EA*EA)+(1-EB)*(1-EB)/(1-EB*EB)-2*(1-EA)*(1-EB)/(1-EA*EB));
xue@1:   for (int i=0; i<Count; i++)
xue@1:   {
xue@1:     FA=FA*EA+*data;
xue@1:     FB=FB*EB+*(data++);
xue@1:     *(dataout++)=(A*FA-B*FB)*NormFactor;
xue@1:   }
xue@1:   return NormFactor;
xue@1: }//ExpOnsetFilter
xue@1: 
Chris@5: /**
xue@1:   function ExpOnsetFilter_balanced: exponential onset filter without starting step response. It
xue@1:   extends the input signal at the front end by bal*Ta samples by repeating the its value at 0, then
xue@1:   applies the onset filter on the extended signal instead.
xue@1: 
xue@1:   In: data[Count]: signal to filter
xue@1:       Tr, Ta: time constants to the impulse response of onset filter, see ExpOnsetFilter().
xue@1:       bal: balancing factor
xue@1:   Out: dataout[Count]: filtered signal, normalized by multiplying a factor.
xue@1: 
xue@1:   Returns the normalization factor. Identical data and dataout is allowed.
xue@1: */
xue@1: double ExpOnsetFilter_balanced(double* dataout, double* data, int Count, double Tr, double Ta, int bal)
xue@1: {
xue@1:   double* tmpdata=new double[int(Count+bal*Ta)];
xue@1:   double* ltmpdata=tmpdata;
xue@1:   for (int i=0; i<bal*Ta; i++) *(ltmpdata++)=data[0];
xue@1:   memcpy(ltmpdata, data, sizeof(double)*Count);
xue@1:   double result=ExpOnsetFilter(tmpdata, tmpdata, bal*Ta+Count, Tr, Ta);
xue@1:   memcpy(dataout, ltmpdata, sizeof(double)*Count);
xue@1:   delete[] tmpdata;
xue@1:   return result;
xue@1: }//ExpOnsetFilter_balanced
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function ExtractLinearComponent: Legendre linear component
xue@1: 
xue@1:   In: data[Count+1]: a signal
xue@1:   Out: dataout[Count+1]: its Legendre linear component, optional.
xue@1: 
xue@1:   Returns the coefficient to the linear component.
xue@1: */
xue@1: double ExtractLinearComponent(double* dataout, double* data, int Count)
xue@1: {
xue@1:   double tmp=0;
xue@1:   int N=Count*2;
xue@1:   for (int n=0; n<=Count; n++) tmp+=n**(data++);
xue@1:   tmp=tmp*24/N/(N+1)/(N+2);
xue@1:   if (dataout)
xue@1:     for (int n=0; n<=Count; n++) *(dataout++)=tmp*n;
xue@1:   return tmp;
xue@1: }//ExtractLinearComponent
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function FFTConv: fast convolution of two series by FFT overlap-add. In an overlap-add scheme it is
xue@1:   assumed that one of the convolvends is short compared to the other one, which can be potentially
xue@1:   infinitely long. The long convolvend is devided into short segments, each of which is convolved with
xue@1:   the short convolvend, the results of which are then assembled into the final result. The minimal delay
xue@1:   from input to output is the amount of overlap, which is the size of the short convolvend minus 1.
xue@1: 
xue@1:   In: source1[size1]: convolvend
xue@1:       source2[size2]: second convolvend
xue@1:       zero: position of first point in convoluton result, relative to main output buffer.
xue@1:       pre_buffer[-zero]: buffer hosting values to be overlap-added to the start of the result.
xue@1:   Out: dest[size1]: the middle part of convolution result
xue@1:        pre_buffer[-zero]: now updated by adding beginning part of the convolution result
xue@1:        post_buffer[size2+zero]: end part of the convolution result
xue@1: 
xue@1:   No return value. Identical dest and source1 allowed.
xue@1: 
xue@1:   The convolution result has length size1+size2 (counting one trailing zero). If zero lies in the range
xue@1:   between -size2 and 0, then the first -zero samples are added to pre_buffer[], next size1 samples are
xue@1:   saved to dest[], and the last size2+zero sampled are saved to post_buffer[]; if not, the middle size1
xue@1:   samples are saved to dest[], while pre_buffer[] and post_buffer[] are not used.
xue@1: */
xue@1: void FFTConv(double* dest, double* source1, int size1, double* source2, int size2, int zero, double* pre_buffer, double* post_buffer)
xue@1: {
xue@11:   int order=Log2(size2-1)+1+1;
xue@1:   int Wid=1<<order;
xue@1:   int HWid=Wid/2;
xue@1:   int Fr=size1/HWid;
xue@1:   int res=size1-HWid*Fr;
xue@1:   bool trunc=false;
xue@1:   if (zero<-size2+1 || zero>0)  zero=-size2/2, trunc=true;
xue@1:   if (pre_buffer==NULL || (post_buffer==NULL && size2+zero!=0)) trunc=true;
xue@1: 
xue@1:   AllocateFFTBuffer(Wid, fft, w, x1);
xue@1:   int* hbitinv=CreateBitInvTable(order-1);
xue@1:   cdouble* x2=new cdouble[Wid];
xue@1:   double* tmp=new double[HWid];
xue@1:   memset(tmp, 0, sizeof(double)*HWid);
xue@1: 
xue@1:   memcpy(fft, source2, sizeof(double)*size2);
xue@1:   memset(&fft[size2], 0, sizeof(double)*(Wid-size2));
xue@1:   RFFTC(fft, 0, 0, order, w, x2, hbitinv);
xue@1: 
xue@1:   double r1, r2, i1, i2;
xue@1:   int ind, ind_;
xue@1:   for (int i=0; i<Fr; i++)
xue@1:   {
xue@1:     memcpy(fft, &source1[i*HWid], sizeof(double)*HWid);
xue@1:     memset(&fft[HWid], 0, sizeof(double)*HWid);
xue@1: 
xue@1:     RFFTC(fft, 0, 0, order, w, x1, hbitinv);
xue@1: 
xue@1:     for (int j=0; j<Wid; j++)
xue@1:     {
xue@1:       r1=x1[j].x, r2=x2[j].x, i1=x1[j].y, i2=x2[j].y;
xue@1:       x1[j].x=r1*r2-i1*i2;
xue@1:       x1[j].y=r1*i2+r2*i1;
xue@1:     }
xue@1:     CIFFTR(x1, order, w, fft, hbitinv);
xue@1:     for (int j=0; j<HWid; j++) tmp[j]+=fft[j];
xue@1: 
xue@1:     ind=i*HWid+zero; //(i+1)*HWid<=size1
xue@1:     ind_=ind+HWid; //ind_=(i+1)*HWid+zero<=size1
xue@1:     if (ind<0)
xue@1:     {
xue@1:       if (!trunc)
xue@1:         memdoubleadd(pre_buffer, tmp, -ind);
xue@1:       memcpy(dest, &tmp[-ind], sizeof(double)*(HWid+ind));
xue@1:     }
xue@1:     else
xue@1:       memcpy(&dest[ind], tmp, sizeof(double)*HWid);
xue@1:     memcpy(tmp, &fft[HWid], sizeof(double)*HWid);
xue@1:   }
xue@1: 
xue@1:   if (res>0)
xue@1:   {
xue@1:     memcpy(fft, &source1[Fr*HWid], sizeof(double)*res);
xue@1:     memset(&fft[res], 0, sizeof(double)*(Wid-res));
xue@1: 
xue@1:     RFFTC(fft, 0, 0, order, w, x1, hbitinv);
xue@1: 
xue@1:     for (int j=0; j<Wid; j++)
xue@1:     {
xue@1:       r1=x1[j].x, r2=x2[j].x, i1=x1[j].y, i2=x2[j].y;
xue@1:       x1[j].x=r1*r2-i1*i2;
xue@1:       x1[j].y=r1*i2+r2*i1;
xue@1:     }
xue@1:     CIFFTR(x1, order, w, fft, hbitinv);
xue@1:     for (int j=0; j<HWid; j++)
xue@1:       tmp[j]+=fft[j];
xue@1: 
xue@1:     ind=Fr*HWid+zero; //Fr*HWid=size1-res, ind=size1-res+zero<size1
xue@1:     ind_=ind+HWid;  //ind_=size1 -res+zero+HWid
xue@1:     if (ind<0)
xue@1:     {
xue@1:       if (!trunc)
xue@1:         memdoubleadd(pre_buffer, tmp, -ind);
xue@1:       memcpy(dest, &tmp[-ind], sizeof(double)*(HWid+ind));
xue@1:     }
xue@1:     else if (ind_>size1)
xue@1:     {
xue@1:       memcpy(&dest[ind], tmp, sizeof(double)*(size1-ind));
xue@1:       if (!trunc && post_buffer)
xue@1:       {
xue@1:         if (ind_>size1+size2+zero)
xue@1:           memcpy(post_buffer, &tmp[size1-ind], sizeof(double)*(size2+zero));
xue@1:         else
xue@1:           memcpy(post_buffer, &tmp[size1-ind], sizeof(double)*(ind_-size1));
xue@1:       }
xue@1:     }
xue@1:     else
xue@1:       memcpy(&dest[ind], tmp, sizeof(double)*HWid);
xue@1:     memcpy(tmp, &fft[HWid], sizeof(double)*HWid);
xue@1:     Fr++;
xue@1:   }
xue@1: 
xue@1:   ind=Fr*HWid+zero;
xue@1:   ind_=ind+HWid;
xue@1: 
xue@1:   if (ind<size1)
xue@1:   {
xue@1:     if (ind_>size1)
xue@1:     {
xue@1:       memcpy(&dest[ind], tmp, sizeof(double)*(size1-ind));
xue@1:       if (!trunc && post_buffer)
xue@1:       {
xue@1:         if (ind_>size1+size2+zero)
xue@1:           memcpy(post_buffer, &tmp[size1-ind], sizeof(double)*(size2+zero));
xue@1:         else
xue@1:           memcpy(post_buffer, &tmp[size1-ind], sizeof(double)*(ind_-size1));
xue@1:       }
xue@1:     }
xue@1:     else
xue@1:       memcpy(&dest[ind], tmp, sizeof(double)*HWid);
xue@1:   }
xue@1:   else //ind>=size1 => ind_>=size1+size2+zero
xue@1:   {
xue@1:     if (!trunc && post_buffer)
xue@1:       memcpy(&post_buffer[ind-size1], tmp, sizeof(double)*(size1+size2+zero-ind));
xue@1:   }
xue@1: 
xue@1:   FreeFFTBuffer(fft);
xue@1:   delete[] x2;
xue@1:   delete[] tmp;
xue@1:   delete[] hbitinv;
xue@1: }//FFTConv
xue@1: 
Chris@5: /**
xue@1:   function FFTConv: the simplified version using two output buffers instead of three. This is almost
xue@1:   equivalent to setting zero=-size2 in the three-output-buffer version (so that post_buffer no longer
xue@1:   exists), except that this version requires size2 (renamed HWid) be a power of 2, and pre_buffer point
xue@1:   to the END of the storage (so that pre_buffer=dest automatically connects the two buffers in a
xue@1:   continuous memory block).
xue@1: 
xue@1:   In: source1[size1]: convolvend
xue@1:       source2[HWid]: second convolved, HWid be a power of 2
xue@1:       pre_buffer[-HWid:-1]: buffer hosting values to be overlap-added to the start of the result.
xue@1:   Out: dest[size1]: main output buffer, now hosting end part of the result (after HWid samples).
xue@1:        pre_buffer[-HWid:-1]: now updated by added the start of the result
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void FFTConv(double* dest, double* source1, int size1, double* source2, int HWid, double* pre_buffer)
xue@1: {
xue@1:   int Wid=HWid*2;
xue@11:   int order=Log2(Wid);
xue@1:   int Fr=size1/HWid;
xue@1:   int res=size1-HWid*Fr;
xue@1: 
xue@1:   AllocateFFTBuffer(Wid, fft, w, x1);
xue@1:   cdouble *x2=new cdouble[Wid];
xue@1:   double *tmp=new double[HWid];
xue@1:   int* hbitinv=CreateBitInvTable(order-1);
xue@1: 
xue@1:   memcpy(fft, source2, sizeof(double)*HWid);
xue@1:   memset(&fft[HWid], 0, sizeof(double)*HWid);
xue@1:   RFFTC(fft, 0, 0, order, w, x2, hbitinv);
xue@1: 
xue@1:   double r1, r2, i1, i2;
xue@1:   for (int i=0; i<Fr; i++)
xue@1:   {
xue@1:     memcpy(fft, &source1[i*HWid], sizeof(double)*HWid);
xue@1:     memset(&fft[HWid], 0, sizeof(double)*HWid);
xue@1: 
xue@1:     RFFTC(fft, 0, 0, order, w, x1, hbitinv);
xue@1: 
xue@1:     for (int j=0; j<Wid; j++)
xue@1:     {
xue@1:       r1=x1[j].x, r2=x2[j].x, i1=x1[j].y, i2=x2[j].y;
xue@1:       x1[j].x=r1*r2-i1*i2;
xue@1:       x1[j].y=r1*i2+r2*i1;
xue@1:     }
xue@1:     CIFFTR(x1, order, w, fft, hbitinv);
xue@1: 
xue@1:     if (i==0)
xue@1:     {
xue@1:       if (pre_buffer!=NULL)
xue@1:       {
xue@1:         double* destl=&pre_buffer[-HWid+1];
xue@1:         for (int j=0; j<HWid-1; j++) destl[j]+=fft[j];
xue@1:       }
xue@1:     }
xue@1:     else
xue@1:     {
xue@1:       for (int j=0; j<HWid-1; j++) tmp[j+1]+=fft[j];
xue@1:       memcpy(&dest[(i-1)*HWid], tmp, sizeof(double)*HWid);
xue@1:     }
xue@1:     memcpy(tmp, &fft[HWid-1], sizeof(double)*HWid);
xue@1:   }
xue@1: 
xue@1:   if (res>0)
xue@1:   {
xue@1:     if (Fr==0) memset(tmp, 0, sizeof(double)*HWid);
xue@1: 
xue@1:     memcpy(fft, &source1[Fr*HWid], sizeof(double)*res);
xue@1:     memset(&fft[res], 0, sizeof(double)*(Wid-res));
xue@1: 
xue@1:     RFFTC(fft, 0, 0, order, w, x1, hbitinv);
xue@1:     for (int j=0; j<Wid; j++)
xue@1:     {
xue@1:       r1=x1[j].x, r2=x2[j].x, i1=x1[j].y, i2=x2[j].y;
xue@1:       x1[j].x=r1*r2-i1*i2;
xue@1:       x1[j].y=r1*i2+r2*i1;
xue@1:     }
xue@1:     CIFFTR(x1, order, w, fft, hbitinv);
xue@1: 
xue@1:     if (Fr==0)
xue@1:     {
xue@1:       if (pre_buffer!=NULL)
xue@1:       {
xue@1:         double* destl=&pre_buffer[-HWid+1];
xue@1:         for (int j=0; j<HWid-1; j++) destl[j]+=fft[j];
xue@1:       }
xue@1:     }
xue@1:     else
xue@1:     {
xue@1:       for (int j=0; j<HWid-1; j++) tmp[j+1]+=fft[j];
xue@1:       memcpy(&dest[(Fr-1)*HWid], tmp, sizeof(double)*HWid);
xue@1:     }
xue@1: 
xue@1:     memcpy(&dest[Fr*HWid], &fft[HWid-1], sizeof(double)*res);
xue@1:   }
xue@1:   else
xue@1:     memcpy(&dest[(Fr-1)*HWid], tmp, sizeof(double)*HWid);
xue@1: 
xue@1:   FreeFFTBuffer(fft);
xue@1:   delete[] x2; delete[] tmp; delete[] hbitinv;
xue@1: }//FFTConv
xue@1: 
Chris@5: /**
xue@1:   function FFTConv: fast convolution of two series by FFT overlap-add. Same as the three-output-buffer
xue@1:   version above but using integer output buffers as well as integer source1.
xue@1: 
xue@1:   In: source1[size1]: convolvend
xue@1:       bps: bytes per sample of integer units in source1[].
xue@1:       source2[size2]: second convolvend
xue@1:       zero: position of first point in convoluton result, relative to main output buffer.
xue@1:       pre_buffer[-zero]: buffer hosting values to be overlap-added to the start of the result.
xue@1:   Out: dest[size1]: the middle part of convolution result
xue@1:        pre_buffer[-zero]: now updated by adding beginning part of the convolution result
xue@1:        post_buffer[size2+zero]: end part of the convolution result
xue@1: 
xue@1:   No return value. Identical dest and source1 allowed.
xue@1: */
xue@1: void FFTConv(unsigned char* dest, unsigned char* source1, int bps, int size1, double* source2, int size2, int zero, unsigned char* pre_buffer, unsigned char* post_buffer)
xue@1: {
xue@11:   int order=Log2(size2-1)+1+1;
xue@1:   int Wid=1<<order;
xue@1:   int HWid=Wid/2;
xue@1:   int Fr=size1/HWid;
xue@1:   int res=size1-HWid*Fr;
xue@1:   bool trunc=false;
xue@1:   if (zero<-size2+1 || zero>0)  zero=-size2/2, trunc=true;
xue@1:   if (pre_buffer==NULL || (post_buffer==NULL && size2+zero!=0)) trunc=true;
xue@1: 
xue@1:   AllocateFFTBuffer(Wid, fft, w, x1);
xue@1:   cdouble* x2=new cdouble[Wid];
xue@1:   double* tmp=new double[HWid];
xue@1:   memset(tmp, 0, sizeof(double)*HWid);
xue@1:   int* hbitinv=CreateBitInvTable(order-1);
xue@1: 
xue@1:   memcpy(fft, source2, sizeof(double)*size2);
xue@1:   memset(&fft[size2], 0, sizeof(double)*(Wid-size2));
xue@1:   RFFTC(fft, 0, 0, order, w, x2, hbitinv);
xue@1: 
xue@1:   double r1, r2, i1, i2;
xue@1:   int ind, ind_;
xue@1:   for (int i=0; i<Fr; i++)
xue@1:   {
xue@1:     IntToDouble(fft, &source1[i*HWid*bps], bps, HWid);
xue@1:     memset(&fft[HWid], 0, sizeof(double)*HWid);
xue@1: 
xue@1:     RFFTC(fft, 0, 0, order, w, x1, hbitinv);
xue@1: 
xue@1:     for (int j=0; j<Wid; j++)
xue@1:     {
xue@1:       r1=x1[j].x, r2=x2[j].x, i1=x1[j].y, i2=x2[j].y;
xue@1:       x1[j].x=r1*r2-i1*i2;
xue@1:       x1[j].y=r1*i2+r2*i1;
xue@1:     }
xue@1:     CIFFTR(x1, order, w, fft, hbitinv);
xue@1:     for (int j=0; j<HWid; j++) tmp[j]+=fft[j];
xue@1: 
xue@1:     ind=i*HWid+zero; //(i+1)*HWid<=size1
xue@1:     ind_=ind+HWid; //ind_=(i+1)*HWid+zero<=size1
xue@1:     if (ind<0)
xue@1:     {
xue@1:       if (!trunc)
xue@1:         DoubleToIntAdd(pre_buffer, bps, tmp, -ind);
xue@1:       DoubleToInt(dest, bps, &tmp[-ind], HWid+ind);
xue@1:     }
xue@1:     else
xue@1:       DoubleToInt(&dest[ind*bps], bps, tmp, HWid);
xue@1:     memcpy(tmp, &fft[HWid], sizeof(double)*HWid);
xue@1:   }
xue@1: 
xue@1:   if (res>0)
xue@1:   {
xue@1:     IntToDouble(fft, &source1[Fr*HWid*bps], bps, res);
xue@1:     memset(&fft[res], 0, sizeof(double)*(Wid-res));
xue@1: 
xue@1:     RFFTC(fft, 0, 0, order, w, x1, hbitinv);
xue@1: 
xue@1:     for (int j=0; j<Wid; j++)
xue@1:     {
xue@1:       r1=x1[j].x, r2=x2[j].x, i1=x1[j].y, i2=x2[j].y;
xue@1:       x1[j].x=r1*r2-i1*i2;
xue@1:       x1[j].y=r1*i2+r2*i1;
xue@1:     }
xue@1:     CIFFTR(x1, order, w, fft, hbitinv);
xue@1:     for (int j=0; j<HWid; j++)
xue@1:       tmp[j]+=fft[j];
xue@1: 
xue@1:     ind=Fr*HWid+zero; //Fr*HWid=size1-res, ind=size1-res+zero<size1
xue@1:     ind_=ind+HWid;  //ind_=size1 -res+zero+HWid
xue@1:     if (ind<0)
xue@1:     {
xue@1:       if (!trunc)
xue@1:         DoubleToIntAdd(pre_buffer, bps, tmp, -ind);
xue@1:       DoubleToInt(dest, bps, &tmp[-ind], HWid+ind);
xue@1:     }
xue@1:     else if (ind_>size1)
xue@1:     {
xue@1:       DoubleToInt(&dest[ind*bps], bps, tmp, size1-ind);
xue@1:       if (!trunc && post_buffer)
xue@1:       {
xue@1:         if (ind_>size1+size2+zero)
xue@1:           DoubleToInt(post_buffer, bps, &tmp[size1-ind], size2+zero);
xue@1:         else
xue@1:           DoubleToInt(post_buffer, bps, &tmp[size1-ind], ind_-size1);
xue@1:       }
xue@1:     }
xue@1:     else
xue@1:       DoubleToInt(&dest[ind*bps], bps, tmp, HWid);
xue@1:     memcpy(tmp, &fft[HWid], sizeof(double)*HWid);
xue@1:     Fr++;
xue@1:   }
xue@1: 
xue@1:   ind=Fr*HWid+zero;
xue@1:   ind_=ind+HWid;
xue@1: 
xue@1:   if (ind<size1)
xue@1:   {
xue@1:     if (ind_>size1)
xue@1:     {
xue@1:       DoubleToInt(&dest[ind*bps], bps, tmp, size1-ind);
xue@1:       if (!trunc && post_buffer)
xue@1:       {
xue@1:         if (ind_>size1+size2+zero)
xue@1:           DoubleToInt(post_buffer, bps, &tmp[size1-ind], size2+zero);
xue@1:         else
xue@1:           DoubleToInt(post_buffer, bps, &tmp[size1-ind], ind_-size1);
xue@1:       }
xue@1:     }
xue@1:     else
xue@1:       DoubleToInt(&dest[ind*bps], bps, tmp, HWid);
xue@1:   }
xue@1:   else //ind>=size1 => ind_>=size1+size2+zero
xue@1:   {
xue@1:     if (!trunc && post_buffer)
xue@1:       DoubleToInt(&post_buffer[(ind-size1)*bps], bps, tmp, size1+size2+zero-ind);
xue@1:   }
xue@1: 
xue@1:   FreeFFTBuffer(fft);
xue@1:   delete[] x2;
xue@1:   delete[] tmp;
xue@1:   delete[] hbitinv;
xue@1: }//FFTConv
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function FFTFilter: FFT with cosine-window overlap-add: This FFT filter is not an actural LTI system,
xue@1:   but an block processing with overlap-add. In this function the blocks are overlapped by 50% and summed
xue@1:   up with Hann windowing.
xue@1: 
xue@1:   In: data[Count]: input data
xue@1:       Wid: DFT size
xue@1:       On, Off: cut-off frequencies of FFT filter. On<Off: band-pass; On>Off: band-stop.
xue@1:   Out: dataout[Count]: filtered data
xue@1: 
xue@1:   No return value. Identical data and dataout allowed
xue@1: */
xue@1: void FFTFilter(double* dataout, double* data, int Count, int Wid, int On, int Off)
xue@1: {
xue@11:   int Order=Log2(Wid);
xue@1:   int HWid=Wid/2;
xue@1:   int Fr=(Count-Wid)/HWid+1;
xue@1:   AllocateFFTBuffer(Wid, ldata, w, x);
xue@1: 
xue@1:   double* win=new double[Wid];
xue@1:   for (int i=0; i<Wid; i++) win[i]=sqrt((1-cos(2*M_PI*i/Wid))/2);
xue@1:   double* tmpdata=new double[HWid];
xue@1:   memset(tmpdata, 0, HWid*sizeof(double));
xue@1: 
xue@1:   for (int i=0; i<Fr; i++)
xue@1:   {
xue@1:     memcpy(ldata, &data[i*HWid], Wid*sizeof(double));
xue@1:     if (i>0)
xue@1:       for (int k=0; k<HWid; k++)
xue@1:         ldata[k]=ldata[k]*win[k];
xue@1:     for (int k=HWid; k<Wid; k++)
xue@1:       ldata[k]=ldata[k]*win[k];
xue@1: 
xue@1:     RFFTC(ldata, NULL, NULL, Order, w, x, 0);
xue@1: 
xue@1:     if (On<Off) //band pass: keep [On, Off) and set other bins to zero
xue@1:     {
xue@1:       memset(x, 0, On*sizeof(cdouble));
xue@1:       if (On>=1)
xue@1:         memset(&x[Wid-On+1], 0, (On-1)*sizeof(cdouble));
xue@1:       if (Off*2<=Wid)
xue@1:         memset(&x[Off], 0, (Wid-Off*2+1)*sizeof(cdouble));
xue@1:     }
xue@1:     else //band stop: set [Off, On) to zero.
xue@1:     {
xue@1:       memset(&x[Off], 0, sizeof(cdouble)*(On-Off));
xue@1:       memset(&x[Wid-On+1], 0, sizeof(double)*(On-Off));
xue@1:     }
xue@1: 
xue@1:     CIFFTR(x, Order, w, ldata);
xue@1: 
xue@1:     if (i>0) for (int k=0; k<HWid; k++) ldata[k]=ldata[k]*win[k];
xue@1:     for (int k=HWid; k<Wid; k++) ldata[k]=ldata[k]*win[k];
xue@1: 
xue@1:     memcpy(&dataout[i*HWid], tmpdata, HWid*sizeof(double));
xue@1:     for (int k=0; k<HWid; k++) dataout[i*HWid+k]+=ldata[k];
xue@1:     memcpy(tmpdata, &ldata[HWid], HWid*sizeof(double));
xue@1:   }
xue@1: 
xue@1:   memcpy(&dataout[Fr*HWid], tmpdata, HWid*sizeof(double));
xue@1:   memset(&dataout[Fr*HWid+HWid], 0, (Count-Fr*HWid-HWid)*sizeof(double));
xue@1: 
xue@1:   delete[] win;
xue@1:   delete[] tmpdata;
xue@1:   FreeFFTBuffer(ldata);
xue@1: }//FFTFilter
xue@1: 
Chris@5: /**
Chris@5:   function FFTFilterOLA: FFTFilter with overlap-add support. This is a true LTI filter whose impulse
xue@1:   response is constructed using IFFT. The filtering is implemented by fast convolution.
xue@1: 
xue@1:   In: data[Count]: input data
xue@1:       Wid: FFT size
xue@1:       On, Off: cut-off frequencies, in bins, of the filter
xue@1:       pre_buffer[Wid]: buffer hosting sampled to be added with the start of output
xue@1:   Out: dataout[Count]: main output buffer, hosting the last $Count samples of output.
xue@1:        pre_buffer[Wid]: now updated by adding the first Wid samples of output
xue@1: 
xue@1:   No return value. The complete output contains Count+Wid effective samples (including final 0); firt
xue@1:   $Wid are added to pre_buffer[], next Count samples saved to dataout[].
xue@1: */
xue@1: void FFTFilterOLA(double* dataout, double* data, int Count, int Wid, int On, int Off, double* pre_buffer)
xue@1: {
xue@1:   AllocateFFTBuffer(Wid, spec, w, x);
xue@1:   memset(x, 0, sizeof(cdouble)*Wid);
xue@1:   for (int i=On+1; i<Off; i++) x[i].x=x[Wid-i].x=1-2*(i%2);
xue@11:   CIFFTR(x, Log2(Wid), w, spec);
xue@1:   FFTConv(dataout, data, Count, spec, Wid, -Wid, pre_buffer, NULL);
xue@1:   FreeFFTBuffer(spec);
xue@1: }//FFTFilterOLA
xue@1: //version for integer input and output, where BytesPerSample specifies the integer format.
xue@1: void FFTFilterOLA(unsigned char* dataout, unsigned char* data, int BytesPerSample, int Count, int Wid, int On, int Off, unsigned char* pre_buffer)
xue@1: {
xue@1:   AllocateFFTBuffer(Wid, spec, w, x);
xue@1:   memset(x, 0, sizeof(cdouble)*Wid);
xue@1:   for (int i=On+1; i<Off; i++) x[i].x=x[Wid-i].x=1-2*(i%2);
xue@11:   CIFFTR(x, Log2(Wid), w, spec);
xue@1:   FFTConv(dataout, data, BytesPerSample, Count, spec, Wid, -Wid, pre_buffer, NULL);
xue@1:   FreeFFTBuffer(spec);
xue@1: }//FFTFilterOLA
xue@1: 
Chris@5: /**
xue@1:   function FFTFilterOLA: FFT filter with overlap-add support.
xue@1: 
xue@1:   In: data[Count]: input data
xue@1:       amp[0:HWid]: amplitude response
xue@1:       ph[0:HWid]: phase response, where ph[0]=ph[HWid]=0;
xue@1:       pre_buffer[Wid]: buffer hosting sampled to be added to the beginning of the output
xue@1:   Out: dataout[Count]: main output buffer, hosting the middle $Count samples of output.
xue@1:        pre_buffer[Wid]: now updated by adding the first Wid/2 samples of output
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void FFTFilterOLA(double* dataout, double* data, int Count, double* amp, double* ph, int Wid, double* pre_buffer)
xue@1: {
xue@1:   int HWid=Wid/2;
xue@1:   AllocateFFTBuffer(Wid, spec, w, x);
xue@1:   x[0].x=amp[0], x[0].y=0;
xue@1:   for (int i=1; i<HWid; i++)
xue@1:   {
xue@1:     x[i].x=x[Wid-i].x=amp[i]*cos(ph[i]);
xue@1:     x[i].y=amp[i]*sin(ph[i]);
xue@1:     x[Wid-i].y=-x[i].y;
xue@1:   }
xue@1:   x[HWid].x=amp[HWid], x[HWid].y=0;
xue@11:   CIFFTR(x, Log2(Wid), w, spec);
xue@1:   FFTConv(dataout, data, Count, spec, Wid, -Wid, pre_buffer, NULL);
xue@1:   FreeFFTBuffer(spec);
xue@1: }//FFTFilterOLA
xue@1: 
Chris@5: /**
xue@1:   function FFTMask: masks a band of a signal with noise
xue@1: 
xue@1:   In: data[Count]: input signal
xue@1:       DigiOn, DigiOff: cut-off frequences of the band to mask
xue@1:       maskcoef: masking noise amplifier. If set to 1 than the mask level is set to the highest signal
xue@1:         level in the masking band.
xue@1:   Out: dataout[Count]: output data
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: double FFTMask(double* dataout, double* data, int Count, int Wid, double DigiOn, double DigiOff, double maskcoef)
xue@1: {
xue@11:   int Order=Log2(Wid);
xue@1:   int HWid=Wid/2;
xue@1:   int Fr=(Count-Wid)/HWid+1;
xue@1:   int On=Wid*DigiOn, Off=Wid*DigiOff;
xue@1:   AllocateFFTBuffer(Wid, ldata, w, x);
xue@1: 
xue@1:   double* winhann=new double[Wid];
xue@1:   double* winhamm=new double[Wid];
xue@1:   for (int i=0; i<Wid; i++)
xue@1:     {winhamm[i]=0.54-0.46*cos(2*M_PI*i/Wid); winhann[i]=(1-cos(2*M_PI*i/Wid))/2/winhamm[i];}
xue@1:   double* tmpdata=new double[HWid];
xue@1:   memset(tmpdata, 0, HWid*sizeof(double));
xue@1:   double max, randfi;
xue@1: 
xue@1:   max=0;
xue@1:   for (int i=0; i<Fr; i++)
xue@1:   {
xue@1:     memcpy(ldata, &data[i*HWid], Wid*sizeof(double));
xue@1:     if (i>0)
xue@1:       for (int k=0; k<HWid; k++)
xue@1:         ldata[k]=ldata[k]*winhamm[k];
xue@1:     for (int k=HWid; k<Wid; k++)
xue@1:       ldata[k]=ldata[k]*winhamm[k];
xue@1: 
xue@1:     RFFTC(ldata, ldata, NULL, Order, w, x, 0);
xue@1: 
xue@1:     for (int k=On; k<Off; k++)
xue@1:     {
xue@1:       x[k].x=x[Wid-k].x=x[k].y=x[Wid-k].y=0;
xue@1:       if (max<ldata[k]) max=ldata[k];
xue@1:     }
xue@1: 
xue@1:     CIFFTR(x, Order, w, ldata);
xue@1: 
xue@1:     if (i>0)
xue@1:       for (int k=0; k<HWid; k++) ldata[k]=ldata[k]*winhann[k];
xue@1:     for (int k=HWid; k<Wid; k++) ldata[k]=ldata[k]*winhann[k];
xue@1: 
xue@1:     for (int k=0; k<HWid; k++) tmpdata[k]+=ldata[k];
xue@1:     memcpy(&dataout[i*HWid], tmpdata, HWid*sizeof(double));
xue@1:     memcpy(tmpdata, &ldata[HWid], HWid*sizeof(double));
xue@1:   }
xue@1:   memcpy(&dataout[Fr*HWid], tmpdata, HWid*sizeof(double));
xue@1: 
xue@1:   max*=maskcoef;
xue@1: 
xue@1:   for (int i=0; i<Wid; i++)
xue@1:     winhann[i]=winhann[i]*winhamm[i];
xue@1: 
xue@1:   for (int i=0; i<Fr; i++)
xue@1:   {
xue@1:     memset(x, 0, sizeof(cdouble)*Wid);
xue@1:     for (int k=On; k<Off; k++)
xue@1:     {
xue@1:       randfi=rand()*M_PI*2/RAND_MAX;
xue@1:       x[k].x=x[Wid-k].x=max*cos(randfi);
xue@1:       x[k].y=max*sin(randfi);
xue@1:       x[Wid-k].y=-x[k].y;
xue@1:     }
xue@1: 
xue@1:     CIFFTR(x, Order, w, ldata);
xue@1: 
xue@1:     if (i>0)
xue@1:       for (int k=0; k<HWid; k++)
xue@1:         ldata[k]=ldata[k]*winhann[k];
xue@1:     for (int k=HWid; k<Wid; k++)
xue@1:       ldata[k]=ldata[k]*winhann[k];
xue@1: 
xue@1:     for (int k=0; k<Wid; k++) dataout[i*HWid+k]+=ldata[k];
xue@1:   }
xue@1: 
xue@1:   memset(&dataout[Fr*HWid+HWid], 0, (Count-Fr*HWid-HWid)*sizeof(double));
xue@1: 
xue@1:   delete[] winhann;
xue@1:   delete[] winhamm;
xue@1:   delete[] tmpdata;
xue@1:   FreeFFTBuffer(ldata);
xue@1: 
xue@1:   return max;
xue@1: }//FFTMask
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function FindInc: find the element in ordered list data that is closest to value.
xue@1: 
xue@1:   In: data[Count]: a ordered list
xue@1:       value: the value to locate in the list
xue@1: 
xue@1:   Returns the index of the element in the sorted list which is closest to $value.
xue@1: */
xue@1: int FindInc(double value, double* data, int Count)
xue@1: {
xue@1:   if (value>=data[Count-1]) return Count-1;
xue@1:   if (value<data[0]) return 0;
xue@1:   int end=InsertInc(value, data, Count, false);
xue@1:   if (fabs(value-data[end-1])<fabs(value-data[end])) return end-1;
xue@1:   else return end;
xue@1: }//FindInc
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function Gaussian: Gaussian function
xue@1: 
xue@1:   In: Vector[Dim]: a vector
xue@1:       Mean[Dim]: mean of Gaussian function
xue@1:       Dev[Fim]: diagonal autocorrelation matrix of Gaussian function
xue@1: 
xue@1:   Returns the value of Gaussian function at Vector[].
xue@1: */
xue@1: double Gaussian(int Dim, double* Vector, double* Mean, double* Dev)
xue@1: {
xue@1:   double bmt=0, tmp;
xue@1:   for (int dim=0; dim<Dim; dim++)
xue@1:   {
xue@1:     tmp=Vector[dim]-Mean[dim];
xue@1:     bmt+=tmp*tmp/Dev[dim];
xue@1:   }
xue@1:   bmt=-bmt/2;
xue@1:   tmp=log(Dev[0]);
xue@1:   for (int dim=1; dim<Dim; dim++) tmp+=log(Dev[dim]);
xue@1:   bmt-=tmp/2;
xue@1:   bmt-=Dim*log(M_PI*2)/2;
xue@1:   bmt=exp(bmt);
xue@1:   return bmt;
xue@1: }//Gaussian
xue@1: 
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function Hamming: calculates the amplitude spectrum of Hamming window at a given frequency
xue@1: 
xue@1:   In: f: frequency
xue@1:       T: size of Hamming window
xue@1: 
xue@1:   Returns the amplitude spectrum at specified frequency.
xue@1: */
xue@1: double Hamming(double f, double T)
xue@1: {
xue@1:   double omg0=2*M_PI/T;
xue@1:   double omg=f*2*M_PI;
xue@1:   cdouble c1, c2, c3;
xue@1:   cdouble nj(0, -1);
xue@1:   cdouble pj(0, 1);
xue@1:   double a=0.54, b=0.46;
xue@1: 
xue@1:   cdouble c=1.0-exp(nj*T*omg);
xue@1:   double half=0.5;
xue@1: 
xue@1:   if (fabs(omg)<1e-100)
xue@1:     c1=a*T;
xue@1:   else
xue@1:     c1=a*c/(pj*omg);
xue@1: 
xue@1:   if (fabs(omg+omg0)<1e-100)
xue@1:     c2=b*0.5*T;
xue@1:   else
xue@1:     c2=c*b*half/(nj*cdouble(omg+omg0));
xue@1: 
xue@1:   if (fabs(omg-omg0)<1e-100)
xue@1:     c3=b*0.5*T;
xue@1:   else
xue@1:     c3=b*c*half/(nj*cdouble(omg-omg0));
xue@1: 
xue@1:   c=c1+c2+c3;
xue@1:   return abs(c);
xue@1: }//Hamming*/
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function HannSq: computes the square norm of Hann window spectrum (window-size-normalized)
xue@1: 
xue@1:   In: x: frequency, in bins
xue@1:       N: size of Hann window
xue@1: 
xue@1:   Return the square norm.
xue@1: */
xue@1: double HannSq(double x, double N)
xue@1: {
xue@1:   double re, im;
xue@1:   double pim=M_PI*x;
xue@1:   double pimf=pim/N;
xue@1:   double pif=M_PI/N;
xue@1: 
xue@1:   double sinpim=sin(pim);
xue@1:   double sinpimf=sin(pimf);
xue@1:   double sinpimplus1f=sin(pimf+pif);
xue@1:   double sinpimminus1f=sin(pimf-pif);
xue@1: 
xue@1:   double spmdivbyspmf, spmdivbyspmpf, spmdivbyspmmf;
xue@1: 
xue@1:   if (sinpimf==0)
xue@1:     spmdivbyspmf=N, spmdivbyspmpf=spmdivbyspmmf=0;
xue@1:   else if (sinpimplus1f==0)
xue@1:     spmdivbyspmpf=-N, spmdivbyspmf=spmdivbyspmmf=0;
xue@1:   else if (sinpimminus1f==0)
xue@1:     spmdivbyspmmf=-N, spmdivbyspmf=spmdivbyspmpf=0;
xue@1:   else
xue@1:     spmdivbyspmf=sinpim/sinpimf, spmdivbyspmpf=sinpim/sinpimplus1f, spmdivbyspmmf=sinpim/sinpimminus1f;
xue@1: 
xue@1:   re=0.5*spmdivbyspmf-0.25*cos(pif)*(spmdivbyspmpf+spmdivbyspmmf);
xue@1:   im=0.25*sin(pif)*(-spmdivbyspmpf+spmdivbyspmmf);
xue@1: 
xue@1:   return (re*re+im*im)/(N*N);
xue@1: }//HannSq
xue@1: 
Chris@5: /**
xue@1:   function Hann: computes the Hann window amplitude spectrum (window-size-normalized).
xue@1: 
xue@1:   In: x: frequency, in bins
xue@1:       N: size of Hann window
xue@1: 
xue@1:   Return the amplitude spectrum evaluated at x. Maximum 0.5 is reached at x=0. Time 2 to normalize
xue@1:   maximum to 1.
xue@1: */
xue@1: double Hann(double x, double N)
xue@1: {
xue@1:   double pim=M_PI*x;
xue@1:   double pif=M_PI/N;
xue@1:   double pimf=pif*x;
xue@1: 
xue@1:   double sinpim=sin(pim);
xue@1:   double tanpimf=tan(pimf);
xue@1:   double tanpimplus1f=tan(pimf+pif);
xue@1:   double tanpimminus1f=tan(pimf-pif);
xue@1: 
xue@1:   double spmdivbyspmf, spmdivbyspmpf, spmdivbyspmmf;
xue@1: 
xue@1:   if (fabs(tanpimf)<1e-10)
xue@1:     spmdivbyspmf=N, spmdivbyspmpf=spmdivbyspmmf=0;
xue@1:   else if (fabs(tanpimplus1f)<1e-10)
xue@1:     spmdivbyspmpf=-N, spmdivbyspmf=spmdivbyspmmf=0;
xue@1:   else if (fabs(tanpimminus1f)<1e-10)
xue@1:     spmdivbyspmmf=-N, spmdivbyspmf=spmdivbyspmpf=0;
xue@1:   else
xue@1:     spmdivbyspmf=sinpim/tanpimf, spmdivbyspmpf=sinpim/tanpimplus1f, spmdivbyspmmf=sinpim/tanpimminus1f;
xue@1: 
xue@1:   double result=0.5*spmdivbyspmf-0.25*(spmdivbyspmpf+spmdivbyspmmf);
xue@1: 
xue@1:   return result/N;
xue@1: }//HannC
xue@1: 
Chris@5: /**
xue@1:   function HxPeak2: fine spectral peak detection. This does detection and high-precision LSE estimation
xue@1:   in one go. However, since in practise most peaks are spurious, LSE estimation is not necessary on
xue@1:   them. Accordingly, HxPeak2 is deprecated in favour of faster but coarser peak picking methods, such as
xue@1:   QIFFT, which leaves fine estimation to a later stage of processing.
xue@1: 
xue@1:   In: F, dF, ddF: pointers to functions that compute LSE peak energy for, plus its 1st and 2nd
xue@1:         derivatives against, a given frequency.
xue@1:       params: pointer to a data structure (l_hx) hosting input data fed to F, dF, and ddF
xue@1:       (st, en): frequency range, in bins, to search for peaks in
xue@1:       epf: convergence detection threshold
xue@1:   Out: hps[return value]: peak frequencies
xue@1:        vps[return value]; peak amplitudes
xue@1: 
xue@1:   Returns the number of peaks detected.
xue@1: */
xue@1: int HxPeak2(double*& hps, double*& vhps, double (*F)(double, void*), double (*dF)(double, void*), double(*ddF)(double, void*), void* params, double st, double en, double epf)
xue@1: {
xue@1:   struct l_hx {int N; union {double B; struct {int k1; int k2;};}; cdouble* x; double dhxpeak; double hxpeak;} *p=(l_hx *)params;
Chris@3:   int dfshift=offsetof(l_hx, dhxpeak);
Chris@3:   int fshift=offsetof(l_hx, hxpeak);
xue@1:   double B=p->B;
xue@1:   int count=0;
xue@1: 
xue@1:   int den=ceil(en), dst=floor(st);
xue@1:   if (den-dst<3) den++, dst--;
xue@1:   if (den-dst<3) den++, dst--;
xue@1:   if (dst<1) dst=1;
xue@1: 
xue@1:   double step=0.5;
xue@1:   int num=(den-dst)/step+1;
xue@1:   bool allochps=false, allocvhps=false;
xue@1:   if (hps==NULL) allochps=true, hps=new double[num];
xue@1:   if (vhps==NULL) allocvhps=true, vhps=new double[num];
xue@1: 
xue@1:   {
xue@1:     double* inp=new double[num];
xue@1:     for (int i=0; i<num; i++)
xue@1:     {
xue@1:       double lf=dst+step*i;
xue@1:       p->k1=ceil(lf-B); if (p->k1<0) p->k1=0;
xue@1:       p->k2=floor(lf+B); if (p->k2>=p->N/2) p->k2=p->N/2-1;
xue@1:       inp[i]=F(lf, params);
xue@1:     }
xue@1: 
xue@1:     for (int i=1; i<num-1; i++)
xue@1:     {
xue@1:       if (inp[i]>=inp[i-1] && inp[i]>=inp[i+1])  //inp[i]=F(dst+step*i)
xue@1:       {
xue@1:         if (inp[i]==inp[i-1] && inp[i]==inp[i+1]) continue;
xue@1:         double fa=dst+step*(i-1), fb=dst+step*(i+1);
xue@1:         double ff=dst+step*i;
xue@1:         p->k1=ceil(ff-B); if (p->k1<0) p->k1=0;
xue@1:         p->k2=floor(ff+B); if (p->k2>=p->N/2) p->k2=p->N/2-1;
xue@1: 
xue@1:         double tmp=Newton1dmax(ff, fa, fb, ddF, params, dfshift, fshift, dF, dfshift, epf);
xue@1: 
xue@1:         //although we have selected inp[i] to be a local maximum, different truncation
xue@1:         //  of local spectrum implies it may not hold as the truncation of inp[i] is
xue@1:         //  used for recalculating inp[i-1] and inp[i+1] in init_Newton method. In this
xue@1:         //  case we retry the sub-maximal frequency to see if it becomes a local maximum
xue@1:         //  when the spectrum is truncated to centre on it.
xue@1: 
xue@1:         if (tmp==-0.5 || tmp==-0.7) //y(fa)<=y(ff)<y(fb) or y(ff)<y(fa)<y(fb)
xue@1:         {
xue@1:           tmp=Newton1dmax(fb, ff, 2*fb-ff, ddF, params, dfshift, fshift, dF, dfshift, epf);
xue@1:           if (tmp==-0.5 || tmp==-0.7) continue;
xue@1:           /*
xue@1:           double ff2=(ff+fb)/2;
xue@1:           p->k1=ceil(ff2-B); if (p->k1<0) p->k1=0;
xue@1:           p->k2=floor(ff2+B); if (p->k2>=p->N/2) p->k2=p->N/2-1;
xue@1:           tmp=Newton1dmax(ff2, ff, fb, ddF, params, dfshift, fshift, dF, dfshift, epf);
xue@1:           p->k1=ceil(ff-B); if (p->k1<0) p->k1=0;
xue@1:           p->k2=floor(ff+B); if (p->k2>=p->N/2) p->k2=p->N/2-1;    */
xue@1:         }
xue@1:         else if (tmp==-0.6 || tmp==-0.8) //y(fb)<=y(ff)<y(fa)
xue@1:         {
xue@1:           tmp=Newton1dmax(fa, 2*fa-ff, ff, ddF, params, dfshift, fshift, dF, dfshift, epf);
xue@1:           if (tmp==-0.6 || tmp==-0.8) continue;
xue@1:         }
xue@1:         if (tmp<0 /*tmp==-0.5 || tmp==-0.6 || tmp==-1 || tmp==-2 || tmp==-3*/)
xue@1:         {
xue@1:           Search1Dmax(ff, params, F, dst+step*(i-1), dst+step*(i+1), &vhps[count], epf);
xue@1:         }
xue@1:         else
xue@1:         {
xue@1:           vhps[count]=p->hxpeak;
xue@1:         }
xue@1:         if (ff>=st && ff<=en && ff>dst+step*(i-0.99) && ff<dst+step*(i+0.99))
xue@1:         {
xue@1: //          if (count==0 || fabs(tmp-hps[count-1])>0.1)
xue@1: //          {
xue@1:             hps[count]=ff;
xue@1:             count++;
xue@1: //          }
xue@1:         }
xue@1:       }
xue@1:     }
xue@1:     delete[] inp;
xue@1:   }
xue@1: 
xue@1:   if (allochps) hps=(double*)realloc(hps, sizeof(double)*count);
xue@1:   if (allocvhps) vhps=(double*)realloc(vhps, sizeof(double)*count);
xue@1:   return count;
xue@1: }//HxPeak2
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function InsertDec: inserts value into sorted decreasing list
xue@1: 
xue@1:   In: data[Count]: a sorted decreasing list.
xue@1:       value: the value to be added
xue@1:   Out: data[Count]: the list with $value inserted if the latter is larger than its last entry, in which
xue@1:         case the original last entry is discarded.
xue@1: 
xue@1:   Returns the index where $value is located in data[], or -1 if $value is smaller than or equal to
xue@1:   data[Count-1].
xue@1: */
xue@1: int InsertDec(int value, int* data, int Count)
xue@1: {
xue@1:   if (Count<=0) return -1;
xue@1:   if (value<=data[Count-1]) return -1;
xue@1:   if (value>data[0])
xue@1:   {
xue@1:     memmove(&data[1], &data[0], sizeof(int)*(Count-1));
xue@1:     data[0]=value;
xue@1:     return 0;
xue@1:   }
xue@1: 
xue@1:   //now Count>=2
xue@1:   int head=0, end=Count-1, mid;
xue@1: 
xue@1:   //D(head)>=value>D(end)
xue@1:   while (end-head>1)
xue@1:   {
xue@1:     mid=(head+end)/2;
xue@1:     if (value<=data[mid]) head=mid;
xue@1:     else end=mid;
xue@1:   }
xue@1: 
xue@1:   //D(head=end-1)>=value>D(end)
xue@1:   memmove(&data[end+1], &data[end], sizeof(int)*(Count-end-1));
xue@1:   data[end]=value;
xue@1:   return end;
xue@1: }//InsertDec
xue@1: //the double version
xue@1: int InsertDec(double value, double* data, int Count)
xue@1: {
xue@1:   if (Count<=0) return -1;
xue@1:   if (value<=data[Count-1]) return -1;
xue@1:   if (value>data[0])
xue@1:   {
xue@1:     memmove(&data[1], &data[0], sizeof(double)*(Count-1));
xue@1:     data[0]=value;
xue@1:     return 0;
xue@1:   }
xue@1: 
xue@1:   //now Count>=2
xue@1:   int head=0, end=Count-1, mid;
xue@1: 
xue@1:   //D(head)>=value>D(end)
xue@1:   while (end-head>1)
xue@1:   {
xue@1:     mid=(head+end)/2;
xue@1:     if (value<=data[mid]) head=mid;
xue@1:     else end=mid;
xue@1:   }
xue@1: 
xue@1:   //D(head=end-1)>=value>D(end)
xue@1:   memmove(&data[end+1], &data[end], sizeof(double)*(Count-end-1));
xue@1:   data[end]=value;
xue@1:   return end;
xue@1: }//InsertDec
xue@1: 
Chris@5: /**
xue@1:   function InsertDec: inserts value and attached integer into sorted decreasing list
xue@1: 
xue@1:   In: data[Count]: a sorted decreasing list
xue@1:       indices[Count]: a list of integers attached to entries of data[]
xue@1:       value, index: the value to be added and its attached integer
xue@1:   Out: data[Count], indices[Count]: the lists with $value and $index inserted if $value is larger than
xue@1:         the last entry of data[], in which case the original last entries are discarded.
xue@1: 
xue@1:   Returns the index where $value is located in data[], or -1 if $value is smaller than or equal to
xue@1:   data[Count-1].
xue@1: */
xue@1: int InsertDec(double value, int index, double* data, int* indices, int Count)
xue@1: {
xue@1:   if (Count<=0) return -1;
xue@1:   if (value<=data[Count-1]) return -1;
xue@1:   if (value>data[0])
xue@1:   {
xue@1:     memmove(&data[1], data, sizeof(double)*(Count-1));
xue@1:     memmove(&indices[1], indices, sizeof(int)*(Count-1));
xue@1:     data[0]=value, indices[0]=index;
xue@1:     return 0;
xue@1:   }
xue@1: 
xue@1:   //now Count>=2
xue@1:   int head=0, end=Count-1, mid;
xue@1: 
xue@1:   //D(head)>=value>D(end)
xue@1:   while (end-head>1)
xue@1:   {
xue@1:     mid=(head+end)/2;
xue@1:     if (value<=data[mid]) head=mid;
xue@1:     else end=mid;
xue@1:   }
xue@1: 
xue@1:   //D(head=end-1)>=value>D(end)
xue@1:   memmove(&data[end+1], &data[end], sizeof(double)*(Count-end-1));
xue@1:   memmove(&indices[end+1], &indices[end], sizeof(int)*(Count-end-1));
xue@1:   data[end]=value, indices[end]=index;
xue@1:   return end;
xue@1: }//InsertDec
xue@1: 
Chris@5: /**
xue@1:   InsertInc: inserts value into sorted increasing list.
xue@1: 
xue@1:   In: data[Count]: a sorted increasing list.
xue@1:       Capacity: maximal size of data[]
xue@1:       value: the value to be added
xue@1:       Compare: pointer to function that compare two values
xue@1:   Out: data[Count]: the list with $value inserted. If the original list is full (Count=Capacity) then
xue@1:         either $value, or the last entry of data[], whichever is larger, is discarded.
xue@1: 
xue@1:   Returns the index where $value is located in data[], or -1 if it is not inserted, which happens if
xue@1:   Count=Capacity and $value is larger than or equal to the last entry in data[Capacity].
xue@1: */
xue@1: int InsertInc(void* value, void** data, int Count, int Capacity, int (*Compare)(void*, void*))
xue@1: {
xue@1:   if (Capacity<=0) return -1;
xue@1:   if (Count>Capacity) Count=Capacity;
xue@1: 
xue@1: //Compare(A,B)<0 if A<B, =0 if A=B, >0 if A>B
xue@1:   int PosToInsert;
xue@1:   if (Count==0) PosToInsert=0;
xue@1:   else if (Compare(value, data[Count-1])>0) PosToInsert=Count;
xue@1:   else if (Compare(value, data[0])<0) PosToInsert=0;
xue@1:   else
xue@1:   {
xue@1:     //now Count>=2
xue@1:     int head=0, end=Count-1, mid;
xue@1: 
xue@1:     //D(head)<=value<D(end)
xue@1:     while (end-head>1)
xue@1:     {
xue@1:       mid=(head+end)/2;
xue@1:       if (Compare(value, data[mid])>=0) head=mid;
xue@1:       else end=mid;
xue@1:     }
xue@1:     //D(head=end-1)<=value<D(end)
xue@1:     PosToInsert=end;
xue@1:   }
xue@1: 
xue@1:   if (Count<Capacity)
xue@1:   {
xue@1:     memmove(&data[PosToInsert+1], &data[PosToInsert], sizeof(void*)*(Count-PosToInsert));
xue@1:     data[PosToInsert]=value;
xue@1:   }
xue@1:   else //Count==Capacity
xue@1:   {
xue@1:     if (PosToInsert>=Capacity) return -1;
xue@1:     memmove(&data[PosToInsert+1], &data[PosToInsert], sizeof(void*)*(Count-PosToInsert-1));
xue@1:     data[PosToInsert]=value;
xue@1:   }
xue@1:   return PosToInsert;
xue@1: }//InsertInc
xue@1: 
Chris@5: /**
xue@1:   function InsertInc: inserts value into sorted increasing list
xue@1: 
xue@1:   In: data[Count]: a sorted increasing list.
xue@1:       value: the value to be added
xue@1:       doinsert: specifies whether the actually insertion is to be performed
xue@1:   Out: data[Count]: the list with $value inserted if the latter is smaller than its last entry, in which
xue@1:         case the original last entry of data[] is discarded.
xue@1: 
xue@1:   Returns the index where $value is located in data[], or -1 if value is larger than or equal to
xue@1:   data[Count-1].
xue@1: */
xue@1: int InsertInc(double value, double* data, int Count, bool doinsert)
xue@1: {
xue@1:   if (Count<=0) return -1;
xue@1:   if (value>=data[Count-1]) return -1;
xue@1:   if (value<data[0])
xue@1:   {
xue@1:     memmove(&data[1], &data[0], sizeof(double)*(Count-1));
xue@1:     if (doinsert) data[0]=value;
xue@1:     return 0;
xue@1:   }
xue@1: 
xue@1:   //now Count>=2
xue@1:   int head=0, end=Count-1, mid;
xue@1: 
xue@1:   //D(head)<=value<D(end)
xue@1:   while (end-head>1)
xue@1:   {
xue@1:     mid=(head+end)/2;
xue@1:     if (value>=data[mid]) head=mid;
xue@1:     else end=mid;
xue@1:   }
xue@1: 
xue@1:   //D(head=end-1)<=value<D(end)
xue@1:   if (doinsert)
xue@1:   {
xue@1:     memmove(&data[end+1], &data[end], sizeof(double)*(Count-end-1));
xue@1:     data[end]=value;
xue@1:   }
xue@1:   return end;
xue@1: }//InsertInc
xue@1: //version where data[] is int.
xue@1: int InsertInc(double value, int* data, int Count, bool doinsert)
xue@1: {
xue@1:   if (Count<=0) return -1;
xue@1:   if (value>=data[Count-1]) return -1;
xue@1:   if (value<data[0])
xue@1:   {
xue@1:     memmove(&data[1], &data[0], sizeof(int)*(Count-1));
xue@1:     if (doinsert) data[0]=value;
xue@1:     return 0;
xue@1:   }
xue@1: 
xue@1:   //now Count>=2
xue@1:   int head=0, end=Count-1, mid;
xue@1: 
xue@1:   //D(head)<=value<D(end)
xue@1:   while (end-head>1)
xue@1:   {
xue@1:     mid=(head+end)/2;
xue@1:     if (value>=data[mid]) head=mid;
xue@1:     else end=mid;
xue@1:   }
xue@1: 
xue@1:   //D(head=end-1)<=value<D(end)
xue@1:   if (doinsert)
xue@1:   {
xue@1:     memmove(&data[end+1], &data[end], sizeof(int)*(Count-end-1));
xue@1:     data[end]=value;
xue@1:   }
xue@1:   return end;
xue@1: }//InsertInc
xue@1: 
Chris@5: /**
xue@1:   function InsertInc: inserts value and attached integer into sorted increasing list
xue@1: 
xue@1:   In: data[Count]: a sorted increasing list
xue@1:       indices[Count]: a list of integers attached to entries of data[]
xue@1:       value, index: the value to be added and its attached integer
xue@1:   Out: data[Count], indices[Count]: the lists with $value and $index inserted if $value is smaller than
xue@1:         the last entry of data[], in which case the original last entries are discarded.
xue@1: 
xue@1:   Returns the index where $value is located in data[], or -1 if $value is larger than or equal to
xue@1:   data[Count-1].
xue@1: */
xue@1: int InsertInc(double value, int index, double* data, int* indices, int Count)
xue@1: {
xue@1:   if (Count<=0) return -1;
xue@1:   if (value>=data[Count-1]) return -1;
xue@1:   if (value<data[0])
xue@1:   {
xue@1:     memmove(&data[1], data, sizeof(double)*(Count-1));
xue@1:     memmove(&indices[1], indices, sizeof(int)*(Count-1));
xue@1:     data[0]=value, indices[0]=index;
xue@1:     return 0;
xue@1:   }
xue@1: 
xue@1:   //now Count>=2
xue@1:   int head=0, end=Count-1, mid;
xue@1: 
xue@1:   //D(head)>=value>D(end)
xue@1:   while (end-head>1)
xue@1:   {
xue@1:     mid=(head+end)/2;
xue@1:     if (value>=data[mid]) head=mid;
xue@1:     else end=mid;
xue@1:   }
xue@1: 
xue@1:   //D(head=end-1)>=value>D(end)
xue@1:   memmove(&data[end+1], &data[end], sizeof(double)*(Count-end-1));
xue@1:   memmove(&indices[end+1], &indices[end], sizeof(int)*(Count-end-1));
xue@1:   data[end]=value, indices[end]=index;
xue@1:   return end;
xue@1: }//InsertInc
xue@1: //version where indices[] is double-precision floating point.
xue@1: int InsertInc(double value, double index, double* data, double* indices, int Count)
xue@1: {
xue@1:   if (Count<=0) return -1;
xue@1:   if (value>=data[Count-1]) return -1;
xue@1:   if (value<data[0])
xue@1:   {
xue@1:     memmove(&data[1], data, sizeof(double)*(Count-1));
xue@1:     memmove(&indices[1], indices, sizeof(double)*(Count-1));
xue@1:     data[0]=value, indices[0]=index;
xue@1:     return 0;
xue@1:   }
xue@1: 
xue@1:   //now Count>=2
xue@1:   int head=0, end=Count-1, mid;
xue@1: 
xue@1:   //D(head)>=value>D(end)
xue@1:   while (end-head>1)
xue@1:   {
xue@1:     mid=(head+end)/2;
xue@1:     if (value>=data[mid]) head=mid;
xue@1:     else end=mid;
xue@1:   }
xue@1: 
xue@1:   //D(head=end-1)>=value>D(end)
xue@1:   memmove(&data[end+1], &data[end], sizeof(double)*(Count-end-1));
xue@1:   memmove(&indices[end+1], &indices[end], sizeof(double)*(Count-end-1));
xue@1:   data[end]=value, indices[end]=index;
xue@1:   return end;
xue@1: }//InsertInc
xue@1: 
Chris@5: /**
xue@1:   function InsertIncApp: inserts value into flexible-length sorted increasing list
xue@1: 
xue@1:   In: data[Count]: a sorted increasing list.
xue@1:       value: the value to be added
xue@1:   Out: data[Count+1]: the list with $value inserted.
xue@1: 
xue@1:   Returns the index where $value is located in data[], or -1 if Count<0. data[] must have Count+1
xue@1:   storage units on calling.
xue@1: */
xue@1: int InsertIncApp(double value, double* data, int Count)
xue@1: {
xue@1:   if (Count<0) return -1;
xue@1:   if (Count==0){data[0]=value; return 0;}
xue@1:   if (value>=data[Count-1]){data[Count]=value; return Count;}
xue@1:   if (value<data[0])
xue@1:   {
xue@1:     memmove(&data[1], &data[0], sizeof(double)*Count);
xue@1:     data[0]=value;
xue@1:     return 0;
xue@1:   }
xue@1: 
xue@1:   //now Count>=2
xue@1:   int head=0, end=Count-1, mid;
xue@1: 
xue@1:   //D(head)<=value<D(end)
xue@1:   while (end-head>1)
xue@1:   {
xue@1:     mid=(head+end)/2;
xue@1:     if (value>=data[mid]) head=mid;
xue@1:     else end=mid;
xue@1:   }
xue@1: 
xue@1:   //D(head=end-1)<=value<D(end)
xue@1:   memmove(&data[end+1], &data[end], sizeof(double)*(Count-end));
xue@1:   data[end]=value;
xue@1: 
xue@1:   return end;
xue@1: }//InsertIncApp
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function InstantFreq; calculates instantaneous frequency from spectrum, evaluated at bin k
xue@1:   
xue@1:   In: x[hwid]: spectrum with scale 2hwid
xue@1:       k: reference frequency, in bins
xue@1:       mode: must be 1.
xue@1:   
xue@1:   Returns an instantaneous frequency near bin k.
xue@1: */
xue@1: double InstantFreq(int k, int hwid, cdouble* x, int mode)
xue@1: {
xue@1:   double result;
xue@1:   switch(mode)
xue@1:   {
xue@1:     //mode 1: the phase vocoder method, based on J. Brown, where the spectrogram
xue@1:     //  MUST be calculated using rectangular window
xue@1:     case 1:
xue@1:     {
xue@1:       if (k<1) k=1;
xue@1:       if (k>hwid-2) k=hwid-2;
xue@1:       double hr=0.5*(x[k].x-0.5*(x[k+1].x+x[k-1].x)), hi=0.5*(x[k].y-0.5*(x[k+1].y+x[k-1].y));
xue@1:       double ph0=Atan2(hi, hr);
xue@1:       double c=cos(M_PI/hwid), s=sin(M_PI/hwid);
xue@1:       hr=0.5*(x[k].x-0.5*(x[k+1].x*c-x[k+1].y*s+x[k-1].x*c+x[k-1].y*s));
xue@1:       hi=0.5*(x[k].y-0.5*(x[k+1].y*c+x[k+1].x*s+x[k-1].y*c-x[k-1].x*s));
xue@1:       double ph1=Atan2(hi, hr);
xue@1:       result=(ph1-ph0)/(2*M_PI);
xue@1:       if (result<-0.5) result+=1;
xue@1:       if (result>0.5) result-=1;
xue@1:       result+=k*0.5/hwid;
xue@1:       break;
xue@1:     }
xue@1:     case 2:
xue@1:       break;
xue@1:   }
xue@1:   return result;
xue@1: }//InstantFreq
xue@1: 
Chris@5: /**
xue@1:   function InstantFreq; calculates "frequency spectrum", a sequence of frequencies evaluated at DFT bins
xue@1:   
xue@1:   In: x[hwid]: spectrum with scale 2hwid
xue@1:       mode: must be 1.
xue@1:   Out: freqspec[hwid]: the frequency spectrum
xue@1:     
xue@1:   No return value.
xue@1: */
xue@1: void InstantFreq(double* freqspec, int hwid, cdouble* x, int mode)
xue@1: {
xue@1:   for (int i=0; i<hwid; i++)
xue@1:     freqspec[i]=InstantFreq(i, hwid, x, mode);
xue@1: }//InstantFreq
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function IntToDouble: copy content of integer array to double array
xue@1: 
xue@1:   In: in: pointer to integer array
xue@1:       BytesPerSample: number of bytes each integer takes
xue@1:       Count: size of integer array, in integers
xue@1:   Out: vector out[Count].
xue@1: 
xue@1:   No return value.
xue@1: 
xue@1:   This version is currently commented out in favour of the version implemented in QuickSpec.cpp which
xue@1:   supports 24-bit integers.
xue@1: *//*
xue@1: void IntToDouble(double* out, void* in, int BytesPerSample, int Count)
xue@1: {
xue@1:   if (BytesPerSample==1){unsigned char* in8=(unsigned char*)in; for (int k=0; k<Count; k++) *(out++)=*(in8++)-128.0;}
xue@1:   else {__int16* in16=(__int16*)in; for (int k=0; k<Count; k++) *(out++)=*(in16++);}
xue@1: }//IntToDouble*/
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function IPHannC: inner product with Hann window spectrum
xue@1: 
xue@1:   In: x[N]: spectrum
xue@1:       f: reference frequency
xue@1:       K1, K2: spectral truncation bounds
xue@1: 
xue@1:   Returns the absolute value of the inner product of x[K1:K2] with the corresponding band of the
xue@1:   spectrum of a sinusoid at frequency f.
xue@1: */
xue@1: double IPHannC(double f, cdouble* x, int N, int K1, int K2)
xue@1: {
xue@1:   int M; double c[4], iH2;
xue@1:   windowspec(wtHann, N, &M, c, &iH2);
xue@1:   return abs(IPWindowC(f, x, N, M, c, iH2, K1, K2));
xue@1: }//IPHannC
xue@1: 
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function lse: linear regression y=ax+b
xue@1: 
xue@1:   In: x[Count], y[Count]: input points
xue@1:   Out: a, b: LSE estimation of coefficients in y=ax+b
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void lse(double* x, double* y, int Count, double& a, double& b)
xue@1: {
xue@1:   double sx=0, sy=0, sxx=0, sxy=0;
xue@1:   for (int i=0; i<Count; i++)
xue@1:   {
xue@1:     sx+=x[i];
xue@1:     sy+=y[i];
xue@1:     sxx+=x[i]*x[i];
xue@1:     sxy+=x[i]*y[i];
xue@1:   }
xue@1:   b=(sxx*sy-sx*sxy)/(Count*sxx-sx*sx);
xue@1:   a=(sy-Count*b)/sx;
xue@1: }//lse
xue@1: 
xue@1: //--------------------------------------------------------------------------
Chris@5: /**
xue@1:   memdoubleadd: vector addition
xue@1: 
xue@1:   In: dest[count], source[count]: addends
xue@1:   Out: dest[count]: sum
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void memdoubleadd(double* dest, double* source, int count)
xue@1: {
xue@1:   for (int i=0; i<count; i++){*dest=*dest+*source; dest++; source++;}
xue@1: }//memdoubleadd
xue@1: 
xue@1: //--------------------------------------------------------------------------
Chris@5: /**
xue@1:   function Mel: converts frequency in Hz to frequency in mel.
xue@1: 
xue@1:   In: f: frequency, in Hz
xue@1: 
xue@1:   Returns the frequency measured on mel scale.
xue@1: */
xue@1: double Mel(double f)
xue@1: {
xue@1:   return 1127.01048*log(1+f/700);
xue@1: }//Mel
xue@1: 
Chris@5: /**
xue@1:   function InvMel: converts frequency in mel to frequency in Hz.
xue@1: 
xue@1:   In: f: frequency, in mel.
xue@1: 
xue@1:   Returns the frequency in Hz.
xue@1: */
xue@1: double InvMel(double mel)
xue@1: {
xue@1:   return 700*(exp(mel/1127.01048)-1);
xue@1: }//InvMel
xue@1: 
Chris@5: /**
xue@1:   function MFCC: calculates MFCC.
xue@1: 
xue@1:   In: Data[FrameWidth]: data
xue@1:       NumBands: number of frequency bands on mel scale
xue@1:       Bands[3*NumBands]: mel frequency bands, comes as $NumBands triples, each containing the lower,
xue@1:         middle and high frequencies, in bins, of one band, from which a weighting window is created to
xue@1:         weight individual bins.
xue@1:       Ceps_Order: number of MFC coefficients (i.e. DCT coefficients)
xue@1:       W, X: FFT buffers
xue@1:   Out: C[Ceps_Order]: MFCC
xue@1:        Amps[NumBands]: log spectrum on MF bands
xue@1: 
xue@1:   No return value. Use MFCCPrepareBands() to retrieve Bands[].
xue@1: */
xue@1: void MFCC(int FrameWidth, int NumBands, int Ceps_Order, double* Data, double* Bands, double* C, double* Amps, cdouble* W, cdouble* X)
xue@1: {
xue@1:   double tmp, b2s, b2c, b2e;
xue@1: 
xue@11:   RFFTC(Data, 0, 0, Log2(FrameWidth), W, X, 0);
xue@1:   for (int i=0; i<=FrameWidth/2; i++) Amps[i]=X[i].x*X[i].x+X[i].y*X[i].y;
xue@1:   
xue@1:   for (int i=0; i<NumBands; i++)
xue@1:   {
xue@1:     tmp=0;
xue@1:     b2s=Bands[3*i], b2c=Bands[3*i+1], b2e=Bands[3*i+2];
xue@1: 
xue@1:     for (int j=ceil(b2s); j<ceil(b2c); j++)
xue@1:       tmp+=Amps[j]*(j-b2s)/(b2c-b2s);
xue@1:     for (int j=ceil(b2c); j<b2e; j++)
xue@1:       tmp+=Amps[j]*(b2e-j)/(b2e-b2c);
xue@1: 
xue@1:     if (tmp<3.7200759760208359629596958038631e-44)
xue@1:       Amps[i]=-100;
xue@1:     else
xue@1:       Amps[i]=log(tmp);
xue@1:   }
xue@1: 
xue@1:   for (int i=0; i<Ceps_Order; i++)
xue@1:   {
xue@1:     tmp=Amps[0]*cos(M_PI*(i+1)/2/NumBands);
xue@1:     for (int j=1; j<NumBands; j++)
xue@1:       tmp+=Amps[j]*cos(M_PI*(i+0.5)*(j+0.5)/NumBands);
xue@1:     C[i]=tmp;
xue@1:   }
xue@1: }//MFCC
xue@1: 
Chris@5: /**
xue@1:   function MFCCPrepareBands: returns a array of OVERLAPPING bands given in triples, whose 1st and 3rd
xue@1:   entries are the start and end of a band, in bins, and the 2nd is a middle frequency.
xue@1: 
xue@1:   In: SamplesPerSec: sampling rate
xue@1:       NumberOfBins: FFT size
xue@1:       NumberOfBands: number of mel-frequency bands
xue@1: 
xue@1:   Returns pointer to the array of triples.
xue@1: */
xue@1: double* MFCCPrepareBands(int NumberOfBands, int SamplesPerSec, int NumberOfBins)
xue@1: {
xue@1:   double* Bands=new double[NumberOfBands*3];
xue@1:   double naqfreq=SamplesPerSec/2.0; //naqvist freq
xue@1:   double binwid=SamplesPerSec*1.0/NumberOfBins;
xue@1:   double naqmel=Mel(naqfreq);
xue@1:   double b=naqmel/(NumberOfBands+1);
xue@1: 
xue@1:   Bands[0]=0;
xue@1:   Bands[1]=InvMel(b)/binwid;
xue@1:   Bands[2]=InvMel(b*2)/binwid;
xue@1:   for (int i=1; i<NumberOfBands; i++)
xue@1:   {
xue@1:     Bands[3*i]=Bands[3*i-2];
xue@1:     Bands[3*i+1]=Bands[3*i-1];
xue@1:     Bands[3*i+2]=InvMel(b*(i+2))/binwid;
xue@1:   }
xue@1:   return Bands;
xue@1: }//MFCCPrepareBands
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function Multi: vector-constant multiplication
xue@1: 
xue@1:   In: data[count]: a vector
xue@1:       mul: a constant
xue@1:   Out: data[count]: their product
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void Multi(double* data, int count, double mul)
xue@1: {
xue@1:   for (int i=0; i<count; i++){*data=*data*mul; data++;}
xue@1: }//Multi
xue@1: 
Chris@5: /**
xue@1:   function Multi: vector-constant multiplication
xue@1: 
xue@1:   In: in[count]: a vector
xue@1:       mul: a constant
xue@1:   Out: out[count]: their product
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void Multi(double* out, double* in, int count, double mul)
xue@1: {
xue@1:   for (int i=0; i<count; i++) *(out++)=*(in++)*mul;
xue@1: }//Multi
xue@1: 
Chris@5: /**
xue@1:   function Multi: vector-constant multiply-addition
xue@1: 
xue@1:   In: in[count], adder[count]: vectors
xue@1:       mul: a constant
xue@1:   Out: out[count]: in[]+adder[]*mul
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void MultiAdd(double* out, double* in, double* adder, int count, double mul)
xue@1: {
xue@1:   for (int i=0; i<count; i++) *(out++)=*(in++)+*(adder++)*mul;
xue@1: }//MultiAdd
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function NearestPeak: finds a peak value in an array that is nearest to a given index
xue@1: 
xue@1:   In: data[count]: an array
xue@1:       anindex: an index
xue@1: 
xue@1:   Returns the index to a peak of data[] that is closest to anindex. In case of two cloest indices,
xue@1:   returns the index to the higher peak of the two.
xue@1: */
xue@1: int NearestPeak(double* data, int count, int anindex)
xue@1: {
xue@1:   int upind=anindex, downind=anindex;
xue@1:   if (anindex<1) anindex=1;
xue@1:   if (anindex>count-2) anindex=count-2;
xue@1: 
xue@1:   if (data[anindex]>data[anindex-1] && data[anindex]>data[anindex+1]) return anindex;
xue@1: 
xue@1:   if (data[anindex]<data[anindex-1])
xue@1:     while (downind>0 && data[downind-1]>data[downind]) downind--;
xue@1:   if (data[anindex]<data[anindex+1])
xue@1:     while (upind<count-1 && data[upind]<data[upind+1]) upind++;
xue@1: 
xue@1:   if (upind==anindex) return downind;
xue@1:   if (downind==anindex) return upind;
xue@1: 
xue@1:   if (anindex-downind<upind-anindex) return downind;
xue@1:   else if (anindex-downind>upind-anindex) return upind;
xue@1:   else if (data[upind]<data[downind]) return downind;
xue@1:   else return upind;
xue@1: }//NearestPeak
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function NegativeExp: fits the curve y=1-x^lmd.
xue@1: 
xue@1:   In: x[Count], y[Count]: sample points to fit, x[0]=0, x[Count-1]=1, y[0]=1, y[Count-1]=0
xue@1:   Out: lmd: coefficient of y=1-x^lmd.
xue@1: 
xue@1:   Returns rms fitting error.
xue@1: */
xue@1: double NegativeExp(double* x, double* y, int Count, double& lmd, int sample, double step, double eps, int maxiter)
xue@1: {
xue@1:   lmd=0;
xue@1:   for (int i=1; i<Count-1; i++)
xue@1:   {
xue@1:     if (y[i]<1)
xue@1:       lmd+=log(1-y[i])/log(x[i]);
xue@1:     else
xue@1:       lmd+=-50/log(x[i]);
xue@1:   }
xue@1:   lmd/=(Count-2);
xue@1: 
xue@1: //lmd has been initialized
xue@1: //coming up will be the recursive calculation of lmd by lgg
xue@1: 
xue@1:   int iter=0;
xue@1:   double laste, lastdel, e=0, del=0, tmp;
xue@1:   do
xue@1:   {
xue@1:     iter++;
xue@1:     laste=e;
xue@1:     lastdel=del;
xue@1:     e=0, del=0;
xue@1:     for (int i=1; i<Count-1; i+=sample)
xue@1:     {
xue@1:       tmp=pow(x[i], lmd);
xue@1:       e=e+(y[i]+tmp-1)*(y[i]+tmp-1);
xue@1:       del=del+(y[i]+tmp-1)*tmp*log(x[i]);
xue@1:     }
xue@1:     if (laste && e>laste) lmd+=step*lastdel, step/=2;
xue@1:     else lmd+=-step*sample*del;
xue@1:   }
xue@1:   while (e && iter<=maxiter && (!laste || fabs(laste-e)/e>eps));
xue@1:   return sqrt(e/Count);
xue@1: }//NegativeExp
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function: NL: noise level, calculated on 5% of total frames with least energy
xue@1: 
xue@1:   In: data[Count]:
xue@1:       Wid: window width for power level estimation
xue@1: 
xue@1:   Returns noise level, in rms.
xue@1: */
xue@1: double NL(double* data, int Count, int Wid)
xue@1: {
xue@1:   int Fr=Count/Wid;
xue@1:   int Num=Fr/20+1;
xue@1:   double* ene=new double[Num], tmp;
xue@1:   for (int i=0; i<Num; i++) ene[i]=1e30;
xue@1:   for (int i=0; i<Fr; i++)
xue@1:   {
xue@1:     tmp=DCPower(&data[i*Wid], Wid, 0);
xue@1:     InsertInc(tmp, ene, Num);
xue@1:   }
xue@1:   double result=Avg(ene, Num, 0);
xue@1:   delete[] ene;
xue@1:   result=sqrt(result);
xue@1:   return result;
xue@1: }//NL
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function Normalize: normalizes data to [-Maxi, Maxi], without zero shift
xue@1: 
xue@1:   In: data[Count]: data to be normalized
xue@1:       Maxi: destination maximal absolute value
xue@1:   Out: data[Count]: normalized data
xue@1: 
xue@1:   Returns the original maximal absolute value.
xue@1: */
xue@1: double Normalize(double* data, int Count, double Maxi)
xue@1: {
xue@1:   double max=0;
xue@1:   double* ldata=data;
xue@1:   for (int i=0; i<Count; i++)
xue@1:   {
xue@1:     if (*ldata>max) max=*ldata;
xue@1:     else if (-*ldata>max) max=-*ldata;
xue@1:     ldata++;
xue@1:   }
xue@1:   if (max>0)
xue@1:   {
xue@1:     Maxi=Maxi/max;
xue@1:     for (int i=0; i<Count; i++) *(data++)*=Maxi;
xue@1:   }
xue@1:   return max;
xue@1: }//Normalize
xue@1: 
Chris@5: /**
xue@1:   function Normalize2: normalizes data to a specified Euclidian norm
xue@1: 
xue@1:   In: data[Count]: data to normalize
xue@1:       Norm: destination Euclidian norm
xue@1:   Out: data[Count]: normalized data.
xue@1: 
xue@1:   Returns the original Euclidian norm.
xue@1: */
xue@1: double Normalize2(double* data, int Count, double Norm)
xue@1: {
xue@1:   double norm=0;
xue@1:   for (int i=0; i<Count; i++) norm+=data[i]*data[i];
xue@1:   norm=sqrt(norm);
xue@1:   double mul=norm/Norm;
xue@1:   if (mul!=0) for (int i=0; i<Count; i++) data[i]/=mul;
xue@1:   return norm;
xue@1: }//Normalize2
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function PhaseSpan: computes the unwrapped phase variation across the Nyquist range
xue@1: 
xue@1:   In: data[Count]: time-domain data
xue@1:       aparams: a fftparams structure
xue@1: 
xue@1:   Returns the difference between unwrapped phase angles at 0 and Nyquist frequency.
xue@1: */
xue@1: double PhaseSpan(double* data, int Count, void* aparams)
xue@1: {
xue@1:   int Pad=1;
xue@1:   fftparams* params=(fftparams*)aparams;
xue@1:   double* Arg=new double[Count*Pad];
xue@1:   AllocateFFTBuffer(Count*Pad, Amp, w, x);
xue@1:   memset(Amp, 0, sizeof(double)*Count*Pad);
xue@1:   memcpy(&Amp[Count*(Pad-1)/2], data, sizeof(double)*Count);
xue@1:   ApplyWindow(Amp, Amp, params->Amp, Count);
xue@11:   RFFTC(Amp, Amp, Arg, Log2(Count*Pad), w, x, 0);
xue@1: 
xue@1:   SmoothPhase(Arg, Count*Pad/2+1);
xue@1:   double result=Arg[Count*Pad/2]-Arg[0];
xue@1:   delete[] Arg;
xue@1:   FreeFFTBuffer(Amp);
xue@1:   return result;
xue@1: }//PhaseSpan
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function PolyFit: least square polynomial fitting y=sum(i){a[i]*x^i}
xue@1: 
xue@1:   In: x[N], y[N]: sample points
xue@1:       P: order of polynomial, P<N
xue@1:   Out: a[P+1]: coefficients of polynomial
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void PolyFit(int P, double* a, int N, double* x, double* y)
xue@1: {
xue@1:   Alloc2(P+1, P+1, aa);
xue@1:   double ai0, bi, *bb=new double[P+1], *s=new double[N], *r=new double[N];
xue@1:   aa[0][0]=N; bi=0; for (int i=0; i<N; i++) s[i]=1, r[i]=y[i], bi+=y[i]; bb[0]=bi;
xue@1: 
xue@1:   for (int i=1; i<=P; i++)
xue@1:   {
xue@1:     ai0=bi=0; for (int j=0; j<N; j++) {s[j]*=x[j], r[j]*=x[j]; ai0+=s[j], bi+=r[j];}
xue@1:     for (int j=0; j<=i; j++) aa[i-j][j]=ai0; bb[i]=bi;
xue@1:   }
xue@1:   for (int i=P+1; i<=2*P; i++)
xue@1:   {
xue@1:     ai0=0; for (int j=0; j<N; j++) {s[j]*=x[j]; ai0+=s[j];}
xue@1:     for (int j=i-P; j<=P; j++) aa[i-j][j]=ai0;
xue@1:   }
xue@1:   GESCP(P+1, a, aa, bb);
xue@1:   DeAlloc2(aa); delete[] bb; delete[] s; delete[] r;
xue@1: }//PolyFit
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function Pow: vector power function
xue@1: 
xue@1:   In: data[Count]: a vector
xue@1:       ex: expontnet
xue@1:   Out: data[Count]: point-wise $ex-th power of data[]
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void Pow(double* data, int Count, double ex)
xue@1: {
xue@1:   for (int i=0; i<Count; i++)
xue@1:     data[i]=pow(data[i], ex);
xue@1: }//Power
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   Rectify: semi-wave rectification
xue@1: 
xue@1:   In: data[Count]: data to rectify
xue@1:       th: rectification threshold: values below th are rectified to th
xue@1:   Out: data[Count]: recitified data
xue@1: 
xue@1:   Returns number of preserved (i.e. not rectified) samples.
xue@1: */
xue@1: int Rectify(double* data, int Count, double th)
xue@1: {
xue@1:   int Result=0;
xue@1:   for (int i=0; i<Count; i++)
xue@1:   {
xue@1:     if (data[i]<=th) data[i]=th;
xue@1:     else Result++;
xue@1:   }
xue@1:   return Result;
xue@1: }//Rectify
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function Res: minimum absolute residue.
xue@1: 
xue@1:   In: in: a number
xue@1:       mod: modulus
xue@1: 
xue@1:   Returns the minimal absolute residue of $in devided by $mod.
xue@1: */
xue@1: double Res(double in, double mod)
xue@1: {
xue@1:   int i=in/mod;
xue@1:   in=in-i*mod;
xue@1:   if (in>mod/2) in-=mod;
xue@1:   if (in<-mod/2) in+=mod;
xue@1:   return in;
xue@1: }//Res
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function Romberg: Romberg algorithm for numerical integration
xue@1: 
xue@1:   In: f: function to integrate
xue@1:       params: extra argument for calling f
xue@1:       a, b: integration boundaries
xue@1:       n: depth of sampling
xue@1: 
xue@1:   Returns the integral of f(*, params) over [a, b].
xue@1: */
xue@1: double Romberg(int n, double(*f)(double, void*), double a, double b, void* params)
xue@1: {
xue@1:   int np=1;
Chris@3:   double* r1=new double[n+1];
xue@1:   double* r2=new double[n+1];
xue@1:   double h=b-a, *swp;
xue@1:   r1[1]=h*(f(a, params)+f(b, params))/2;
xue@1:   for (int i=2; i<=n; i++)
xue@1:   {
xue@1:     double akh=a+0.5*h; r2[1]=f(akh, params);
xue@1:     for (int k=2; k<=np; k++) {akh+=h; r2[1]+=f(akh, params);} //akh=a+(k-0.5)h
xue@1:     r2[1]=0.5*(r1[1]+h*r2[1]);
xue@1:     double fj=4;
xue@1:     for (int j=2; j<=i; j++) {r2[j]=(fj*r2[j-1]-r1[j-1])/(fj-1); fj*=4;}  //fj=4^(j-1)
xue@1:     h/=2; np*=2;
xue@1:     swp=r1; r1=r2; r2=swp;
xue@1:   }
xue@1:   h=r1[n];
xue@1:   delete[] r1;
xue@1:   delete[] r2;
xue@1:   return h;
xue@1: }//Romberg
xue@1: 
Chris@5: /**
xue@1:   function Romberg: Romberg algorithm for numerical integration, may return before specified depth on
xue@1:   convergence.
xue@1: 
xue@1:   In: f: function to integrate
xue@1:       params: extra argument for calling f
xue@1:       a, b: integration boundaries
xue@1:       n: depth of sampling
xue@1:       ep: convergence test threshold
xue@1: 
xue@1:   Returns the integral of f(*, params) over [a, b].
xue@1: */
xue@1: double Romberg(double(*f)(double, void*), double a, double b, void* params, int n, double ep)
xue@1: {
xue@1:   int i, np=1;
xue@1:   double* r1=new double[n+1];
xue@1:   double* r2=new double[n+1];
xue@1:   double h=b-a, *swp;
xue@1:   r1[1]=h*(f(a, params)+f(b, params))/2;
xue@1:   bool mep=false;
xue@1:   for (i=2; i<=n; i++)
xue@1:   {
xue@1:     double akh=a+0.5*h; r2[1]=f(akh, params);
xue@1:     for (int k=2; k<=np; k++) {akh+=h; r2[1]+=f(akh, params);} //akh=a+(k-0.5)h
xue@1:     r2[1]=0.5*(r1[1]+h*r2[1]);
xue@1:     double fj=4;
xue@1:     for (int j=2; j<=i; j++) {r2[j]=(fj*r2[j-1]-r1[j-1])/(fj-1); fj*=4;}  //fj=4^(j-1)
xue@1: 
xue@1:     if (fabs(r2[i]-r1[i-1])<ep)
xue@1:     {
xue@1:       if (mep) break;
xue@1:       else mep=true;
xue@1:     }
xue@1:     else mep=false;
xue@1: 
xue@1:     h/=2; np*=2;
xue@1:     swp=r1; r1=r2; r2=swp;
xue@1:   }
xue@1:   if (i<=n) h=r2[i];
xue@1:   else h=r1[n];
xue@1:   delete[] r1;
xue@1:   delete[] r2;
xue@1:   return h;
xue@1: }//Romberg
xue@1: 
xue@1: //---------------------------------------------------------------------------
xue@1: //analog and digital sinc functions
xue@1: 
xue@1: //sinca(0)=1, sincd(0)=N, sinca(1)=sincd(1)=0.
Chris@5: /**
xue@1:   function sinca: analog sinc function.
xue@1: 
xue@1:   In: x: frequency
xue@1: 
xue@1:   Returns sinc(x)=sin(pi*x)/(pi*x), sinca(0)=1, sinca(1)=0
xue@1: */
xue@1: double sinca(double x)
xue@1: {
xue@1:   if (x==0) return 1;
xue@1:   return sin(M_PI*x)/(M_PI*x);
xue@1: }//sinca
xue@1: 
Chris@5: /**
xue@1:   function sincd_unn: unnormalized discrete sinc function
xue@1: 
xue@1:   In: x: frequency
xue@1:       N: scale (window size, DFT size)
xue@1: 
xue@1:   Returns sinc(x)=sin(pi*x)/sin(pi*x/N), sincd(0)=N, sincd(1)=0.
xue@1: */
xue@1: double sincd_unn(double x, int N)
xue@1: {
xue@1:   if (x==0) return N;
xue@1:   return sin(M_PI*x)/sin(M_PI*x/N);
xue@1: }//sincd
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   SmoothPhase: phase unwrapping on module mpi*PI, 2PI by default
xue@1: 
xue@1:   In: Arg[Count]: phase angles to unwrap
xue@1:       mpi: unwrapping modulus, in pi's
xue@1:   Out: Arg[Count]: unwrapped phase
xue@1: 
xue@1:   Returns the amount of unwrap, in pi's, of the last phase angle
xue@1: */
xue@1: double SmoothPhase(double* Arg, int Count, int mpi)
xue@1: {
xue@1:   double m2pi=mpi*M_PI;
xue@1:   for (int i=1; i<Count-1; i++)
xue@1:     Arg[i]=Arg[i-1]+Res(Arg[i]-Arg[i-1], m2pi);
xue@1:   double tmp=Res(Arg[Count-1]-Arg[Count-2], m2pi);
xue@1:   double result=(Arg[Count-1]-Arg[Count-2]-tmp)/m2pi;
xue@1:   Arg[Count-1]=Arg[Count-2]+tmp;
xue@1: 
xue@1:   return result;
xue@1: }//SmoothPhase
xue@1: 
xue@1: //---------------------------------------------------------------------------
xue@1: //the stiff string partial frequency model f[m]=mf[1]*sqrt(1+B(m*m-1)).
xue@1: 
Chris@5: /**
xue@1:   StiffB: computes stiffness coefficient from fundamental and another partial frequency based on the
xue@1:   stiff string partial frequency model f[m]=mf[1]*sqrt(1+B(m*m-1)).
xue@1: 
xue@1:   In: f0: fundamental frequency
xue@1:       m: 1-based partial index
xue@1:       fm: frequency of partial m
xue@1: 
xue@1:   Returns stiffness coefficient B.
xue@1: */
xue@1: double StiffB(double f0, double fm, int m)
xue@1: {
xue@1:   double f2=fm/m/f0;
xue@1:   return (f2*f2-1)/(m*m-1);
xue@1: }//StiffB
xue@1: 
xue@1: //StiffF: partial frequency of a stiff string
Chris@5: /**
xue@1:   StiffFm: computes a partial frequency from fundamental frequency and partial index based on the stiff
xue@1:   string partial frequency model f[m]=mf[1]*sqrt(1+B(m*m-1)).
xue@1: 
xue@1:   In: f0: fundamental frequency
xue@1:       m: 1-based partial index
xue@1:       B: stiffness coefficient
xue@1: 
xue@1:   Returns frequency of the m-th partial.
xue@1: */
xue@1: double StiffFm(double f0, int m, double B)
xue@1: {
xue@1:   return m*f0*sqrt(1+B*(m*m-1));
xue@1: }//StiffFm
xue@1: 
Chris@5: /**
xue@1:   StiffF0: computes fundamental frequency from another partial frequency and stiffness coefficient based
xue@1:   on the stiff string partial frequency model f[m]=mf[1]*sqrt(1+B(m*m-1)).
xue@1: 
xue@1:   In: m: 1-based partial index
xue@1:       fm: frequency of partial m
xue@1:       B: stiffness coefficient
xue@1: 
xue@1:   Returns the fundamental frequency.
xue@1: */
xue@1: double StiffF0(double fm, int m, double B)
xue@1: {
xue@1:   return fm/m/sqrt(1+B*(m*m-1));
xue@1: }//StiffF0
xue@1: 
Chris@5: /**
xue@1:   StiffM: computes 1-based partial index from partial frequency, fundamental frequency and stiffness
xue@1:   coefficient based on the stiff string partial frequency model f[m]=mf[1]*sqrt(1+B(m*m-1)).
xue@1: 
xue@1:   In: f0: fundamental freqency
xue@1:       fm: a partial frequency
xue@1:       B: stiffness coefficient
xue@1: 
xue@1:   Returns the 1-based partial index which will generate the specified partial frequency from the model
xue@1:   which, however, does not have to be an integer in this function.
xue@1: */
xue@1: double StiffM(double f0, double fm, double B)
xue@1: {
xue@1:   if (B<1e-14) return fm/f0;
xue@1:   double b1=B-1, ff=fm/f0;
xue@1:   double delta=b1*b1+4*B*ff*ff;
xue@1:   if (delta<0)
xue@1:     return sqrt(b1/2/B);
xue@1:   else
xue@1:     return sqrt((b1+sqrt(delta))/2/B);
xue@1: }//StiffMd
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   TFFilter: time-frequency filtering with Hann-windowed overlap-add.
xue@1: 
xue@1:   In: data[Count]: input data
xue@1:       Spans: time-frequency spans
xue@1:       wt, windp: type and extra parameter of FFT window
xue@1:       Sps: sampling rate
xue@1:       TOffst: optional offset applied to all time values in Spans, set to Spans timing of of data[0].
xue@1:       Pass: set to pass T-F content covered by Spans, clear to stop T-F content covered by Spans
xue@1:   Out: dataout[Count]: filtered data
xue@1: 
xue@1:   No return value. Identical data and dataout allowed.
xue@1: */
xue@1: void TFFilter(double* data, double* dataout, int Count, int Wid, TTFSpans* Spans, bool Pass, WindowType wt, double windp, int Sps, int TOffst)
xue@1: {
xue@1:   int HWid=Wid/2;
xue@1:   int Fr=Count/HWid-1;
xue@11:   int Order=Log2(Wid);
xue@1: 
xue@1:   int** lspan=new int*[Fr];
xue@1:   double* lxspan=new double[Fr];
xue@1: 
xue@1:   lspan[0]=new int[Fr*Wid];
xue@1:   for (int i=1; i<Fr; i++)
xue@1:     lspan[i]=&lspan[0][i*Wid];
xue@1: 
xue@1:   //fill local filter span table
xue@1:   if (Pass)
xue@1:     memset(lspan[0], 0, sizeof(int)*Fr*Wid);
xue@1:   else
xue@1:     for (int i=0; i<Fr; i++)
xue@1:       for (int j=0; j<Wid; j++)
xue@1:         lspan[i][j]=1;
xue@1: 
xue@1:   for (int i=0; i<Spans->Count; i++)
xue@1:   {
xue@1:     int lx1, lx2, ly1, ly2;
xue@1:     lx1=(Spans->Items[i].T1-TOffst)/HWid-1;
xue@1:     lx2=(Spans->Items[i].T2-1-TOffst)/HWid;
xue@1:     ly1=Spans->Items[i].F1*2/Sps*HWid+0.001;
xue@1:     ly2=Spans->Items[i].F2*2/Sps*HWid+1;
xue@1:     if (lx1<0) lx1=0;
xue@1:     if (lx2>=Fr) lx2=Fr-1;
xue@1:     if (ly1<0) ly1=0;
xue@1:     if (ly1>HWid) ly1=HWid;
xue@1:     if (Pass)
xue@1:       for (int x=lx1; x<=lx2; x++)
xue@1:         for (int y=ly1; y<=ly2; y++)
xue@1:           lspan[x][y]=1;
xue@1:     else
xue@1:       for (int x=lx1; x<=lx2; x++)
xue@1:         for (int y=ly1; y<=ly2; y++)
xue@1:           lspan[x][y]=0;
xue@1:   }
xue@1:   for (int i=0; i<Fr; i++)
xue@1:   {
xue@1:     lxspan[i]=0;
xue@1:     for (int j=0; j<=HWid; j++)
xue@1:     {
xue@1:       if (lspan[i][j])
xue@1:         lxspan[i]=lxspan[i]+1;
xue@1:     }
xue@1:     lxspan[i]/=(HWid+1);
xue@1:   }
xue@1:   double* winf=NewWindow(wt, Wid, 0, &windp);
xue@1:   double* wini=NewWindow(wtHann, Wid, NULL, NULL);
xue@1:   for (int i=0; i<Wid; i++)
xue@1:     if (winf[i]!=0) wini[i]=wini[i]/winf[i];
xue@1:   AllocateFFTBuffer(Wid, ldata, w, x);
xue@1:   double* tmpdata=new double[HWid];
xue@1:   memset(tmpdata, 0, HWid*sizeof(double));
xue@1: 
xue@1:   for (int i=0; i<Fr; i++)
xue@1:   {
xue@1:     if (lxspan[i]==0)
xue@1:     {
xue@1:       memcpy(&dataout[i*HWid], tmpdata, sizeof(double)*HWid);
xue@1:       memset(tmpdata, 0, sizeof(double)*HWid);
xue@1:       continue;
xue@1:     }
xue@1:     if (lxspan[i]==1)
xue@1:     {
xue@1:       memcpy(ldata, &data[i*HWid], Wid*sizeof(double));
xue@1:       if (i>0)
xue@1:         for (int k=0; k<HWid; k++)
xue@1:           ldata[k]=ldata[k]*winf[k]*wini[k];
xue@1:       for (int k=HWid; k<Wid; k++)
xue@1:         ldata[k]=ldata[k]*winf[k]*wini[k];
xue@1: 
xue@1:       memcpy(&dataout[i*HWid], tmpdata, HWid*sizeof(double));
xue@1:       for (int k=0; k<HWid; k++)
xue@1:         dataout[i*HWid+k]+=ldata[k];
xue@1:       memcpy(tmpdata, &ldata[HWid], HWid*sizeof(double));
xue@1:       continue;
xue@1:     }
xue@1:     memcpy(ldata, &data[i*HWid], Wid*sizeof(double));
xue@1:     if (i>0)
xue@1:       for (int k=0; k<HWid; k++)
xue@1:         ldata[k]=ldata[k]*winf[k];
xue@1:     for (int k=HWid; k<Wid; k++)
xue@1:       ldata[k]=ldata[k]*winf[k];
xue@1: 
xue@1:     RFFTC(ldata, NULL, NULL, Order, w, x, 0);
xue@1: 
xue@1:     if (!lspan[i][0]) x[0].x=x[0].y=0;
xue@1:     for (int j=1; j<=HWid; j++)
xue@1:       if (!lspan[i][j]) x[j].x=x[Wid-j].x=x[j].y=x[Wid-j].y=0;
xue@1: 
xue@1:     CIFFTR(x, Order, w, ldata);
xue@1: 
xue@1:     if (i>0)
xue@1:       for (int k=0; k<HWid; k++)
xue@1:         ldata[k]=ldata[k]*wini[k];
xue@1:     for (int k=HWid; k<Wid; k++)
xue@1:       ldata[k]=ldata[k]*wini[k];
xue@1: 
xue@1:     memcpy(&dataout[i*HWid], tmpdata, HWid*sizeof(double));
xue@1:     for (int k=0; k<HWid; k++)
xue@1:       dataout[i*HWid+k]+=ldata[k];
xue@1:     memcpy(tmpdata, &ldata[HWid], HWid*sizeof(double));
xue@1:   }
xue@1:   memcpy(&dataout[Fr*HWid], tmpdata, sizeof(double)*HWid);
xue@1:   memset(&dataout[Fr*HWid+HWid], 0, sizeof(double)*(Count-Fr*HWid-HWid));
xue@1: 
xue@1:   FreeFFTBuffer(ldata);
xue@1:   delete[] lspan[0];
xue@1:   delete[] lspan;
xue@1:   delete[] lxspan;
xue@1:   delete[] tmpdata;
xue@1:   delete[] winf;
xue@1:   delete[] wini;
xue@1: }//TFFilter
xue@1: //version on integer data, where BytesPerSample specified the integer format.
xue@1: void TFFilter(void* data, void* dataout, int BytesPerSample, int Count, int Wid, TTFSpans* Spans, bool Pass, WindowType wt, double windp, int Sps, int TOffst)
xue@1: {
xue@1:   int HWid=Wid/2;
xue@1:   int Fr=Count/HWid-1;
xue@11:   int Order=Log2(Wid);
xue@1: 
xue@1:   int** lspan=new int*[Fr];
xue@1:   double* lxspan=new double[Fr];
xue@1: 
xue@1:   lspan[0]=new int[Fr*Wid];
xue@1:   for (int i=1; i<Fr; i++)
xue@1:     lspan[i]=&lspan[0][i*Wid];
xue@1: 
xue@1:   //fill local filter span table
xue@1:   if (Pass)
xue@1:     memset(lspan[0], 0, sizeof(int)*Fr*Wid);
xue@1:   else
xue@1:     for (int i=0; i<Fr; i++)
xue@1:       for (int j=0; j<Wid; j++)
xue@1:         lspan[i][j]=1;
xue@1: 
xue@1:   for (int i=0; i<Spans->Count; i++)
xue@1:   {
xue@1:     int lx1, lx2, ly1, ly2;
xue@1:     lx1=(Spans->Items[i].T1-TOffst)/HWid-1;
xue@1:     lx2=(Spans->Items[i].T2-1-TOffst)/HWid;
xue@1:     ly1=Spans->Items[i].F1*2/Sps*HWid+0.001;
xue@1:     ly2=Spans->Items[i].F2*2/Sps*HWid+1;
xue@1:     if (lx1<0) lx1=0;
xue@1:     if (lx2>=Fr) lx2=Fr-1;
xue@1:     if (ly1<0) ly1=0;
xue@1:     if (ly1>HWid) ly1=HWid;
xue@1:     if (Pass)
xue@1:       for (int x=lx1; x<=lx2; x++)
xue@1:         for (int y=ly1; y<=ly2; y++)
xue@1:           lspan[x][y]=1;
xue@1:     else
xue@1:       for (int x=lx1; x<=lx2; x++)
xue@1:         for (int y=ly1; y<=ly2; y++)
xue@1:           lspan[x][y]=0;
xue@1:   }
xue@1:   for (int i=0; i<Fr; i++)
xue@1:   {
xue@1:     lxspan[i]=0;
xue@1:     for (int j=0; j<=HWid; j++)
xue@1:     {
xue@1:       if (lspan[i][j])
xue@1:         lxspan[i]=lxspan[i]+1;
xue@1:     }
xue@1:     lxspan[i]/=(HWid+1);
xue@1:   }
xue@1:   double* winf=NewWindow(wt, Wid, 0, &windp);
xue@1:   double* wini=NewWindow(wtHann, Wid, NULL, NULL);
xue@1:   for (int i=0; i<Wid; i++)
xue@1:     if (winf[i]!=0) wini[i]=wini[i]/winf[i];
xue@1:   AllocateFFTBuffer(Wid, ldata, w, x);
xue@1:   double* tmpdata=new double[HWid];
xue@1:   memset(tmpdata, 0, HWid*sizeof(double));
xue@1: 
xue@1:   for (int i=0; i<Fr; i++)
xue@1:   {
xue@1:     if (lxspan[i]==0)
xue@1:     {
xue@1:       DoubleToInt(&((char*)dataout)[i*HWid*BytesPerSample], BytesPerSample, tmpdata, HWid);
xue@1:       memset(tmpdata, 0, sizeof(double)*HWid);
xue@1:       continue;
xue@1:     }
xue@1:     if (lxspan[i]==1)
xue@1:     {
xue@1:       IntToDouble(ldata, &((char*)data)[i*HWid*BytesPerSample], BytesPerSample, Wid);
xue@1:       if (i>0)
xue@1:         for (int k=0; k<HWid; k++)
xue@1:           ldata[k]=ldata[k]*winf[k]*wini[k];
xue@1:       for (int k=HWid; k<Wid; k++)
xue@1:         ldata[k]=ldata[k]*winf[k]*wini[k];
xue@1: 
xue@1:       for (int k=0; k<HWid; k++) tmpdata[k]+=ldata[k];
xue@1:       DoubleToInt(&((char*)dataout)[i*HWid*BytesPerSample], BytesPerSample, tmpdata, HWid);
xue@1:       memcpy(tmpdata, &ldata[HWid], HWid*sizeof(double));
xue@1:       continue;
xue@1:     }
xue@1:     IntToDouble(ldata, &((char*)data)[i*HWid*BytesPerSample], BytesPerSample, Wid);
xue@1:     if (i>0)
xue@1:       for (int k=0; k<HWid; k++)
xue@1:         ldata[k]=ldata[k]*winf[k];
xue@1:     for (int k=HWid; k<Wid; k++)
xue@1:       ldata[k]=ldata[k]*winf[k];
xue@1: 
xue@1:     RFFTC(ldata, NULL, NULL, Order, w, x, 0);
xue@1: 
xue@1:     if (!lspan[i][0]) x[0].x=x[0].y=0;
xue@1:     for (int j=1; j<=HWid; j++)
xue@1:       if (!lspan[i][j]) x[j].x=x[Wid-j].x=x[j].y=x[Wid-j].y=0;
xue@1: 
xue@1:     CIFFTR(x, Order, w, ldata);
xue@1: 
xue@1:     if (i>0)
xue@1:       for (int k=0; k<HWid; k++)
xue@1:         ldata[k]=ldata[k]*wini[k];
xue@1:     for (int k=HWid; k<Wid; k++)
xue@1:       ldata[k]=ldata[k]*wini[k];
xue@1: 
xue@1: 
xue@1:     for (int k=0; k<HWid; k++)
xue@1:       tmpdata[k]+=ldata[k];
xue@1:     DoubleToInt(&((char*)dataout)[i*HWid*BytesPerSample], BytesPerSample, tmpdata, HWid);
xue@1:     memcpy(tmpdata, &ldata[HWid], HWid*sizeof(double));
xue@1:   }
xue@1:   DoubleToInt(&((char*)dataout)[Fr*HWid*BytesPerSample], BytesPerSample, tmpdata, HWid);
xue@1:   memset(&((char*)dataout)[(Fr*HWid+HWid)*BytesPerSample], 0, BytesPerSample*(Count-Fr*HWid-HWid));
xue@1: 
xue@1:   FreeFFTBuffer(ldata);
xue@1: 
xue@1:   delete[] lspan[0];
xue@1:   delete[] lspan;
xue@1:   delete[] lxspan;
xue@1:   delete[] tmpdata;
xue@1:   delete[] winf;
xue@1:   delete[] wini;
xue@1: }//TFFilter
xue@1: 
xue@1: 
xue@1: