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>
xue@1: #include "vibrato.h"
xue@1: #include "arrayalloc.h"
xue@1: #include "fft.h"
xue@1: #include "hssf.h"
Chris@2: #include "matrix.h"
xue@1: #include "splines.h"
xue@1: #include "xcomplex.h"
Chris@5: 
Chris@5: /** \file vibrato.h */
Chris@5: 
xue@1: //---------------------------------------------------------------------------
xue@1: //method TVo::TVo: default constructor
xue@1: TVo::TVo()
xue@1: {
xue@1: 	memset(this, 0, sizeof(TVo));
xue@1: 	afres=8192;
xue@1: }//TVo
xue@1: 
xue@1: //method TVo::~TVo: default destructor
xue@1: TVo::~TVo()
xue@1: {
xue@1:   free(peakfr);
xue@1: 	delete[] lp;
xue@1:   delete[] F0;
xue@1: 	if (fxr) DeAlloc2(fxr);
xue@1: 	if (fres) DeAlloc2(fres);
xue@1: 	if (LogASp) DeAlloc2(LogASp);
xue@1: }//~TVo
xue@1: 
Chris@5: /**
xue@1:   method TVo::AllocateL: allocates or reallocates storage space whose size depends on L. This uses an
xue@1:   external L value for allocation and updates L itself.
xue@1: */
xue@1: void TVo::AllocateL(int AnL)
xue@1: {
xue@1:   L=AnL;
xue@1:   delete[] F0;
xue@1:   F0=new double[(L+2)*4];
xue@1:   F0C=&F0[L+2], F0D=&F0[(L+2)*2], A0C=&F0[(L+2)*3];
xue@1: }//AllocateL
xue@1: 
Chris@5: /**
xue@1:   method TVo::AllocateP: allocates or reallocates storage space whose size depends on P, using the
xue@1:   interal P.
xue@1: */
xue@1: void TVo::AllocateP()
xue@1: {
xue@1:   peakfr=(int*)realloc(peakfr, sizeof(int)*P);
xue@1:   delete[] lp; lp=new double[(P+2)*3];
xue@1:   Dp=&lp[P+2];
xue@1:   Kp=(int*)&lp[(P+2)*2];
xue@1:   if (LogASp) DeAlloc2(LogASp); Allocate2(double, P, M, LogASp);
xue@1: }//AllocateP
xue@1: 
Chris@5: /**
xue@1:   method TVo::AllocatePK: allocates or reallocates storage space whose size depends on both P and K.
xue@1: */
xue@1: void TVo::AllocatePK()
xue@1: {
xue@1:   if (fxr) DeAlloc2(fxr);
xue@1:   Allocate2(double, P, 2*K+1, fxr);
xue@1:   memset(fxr[0], 0, sizeof(double)*P*(2*K+1));
xue@1:   if (fres) DeAlloc2(fres);
xue@1:   Allocate2(double, P, K, fres);
xue@1: }//AllocatePK
xue@1: 
Chris@5: /**
xue@1:   method TVo::GetOldP: returns the number of "equivalent" cycles of the whole duration. Although P-1
xue@1:   cycles are delineated by lp[P], the parts before lp[0] and after lp[P-1] are also a part of the
xue@1:   harmonic sinusoid being modeled and therefore should be taken into account when necessary.
xue@1: */
xue@1: double TVo::GetOldP()
xue@1: {
xue@1:   return P-1+lp[0]/(lp[1]-lp[0])+(L-lp[P-1])/(lp[P-1]-lp[P-2]);
xue@1: }//GetOldP
xue@1: 
Chris@5: /**
xue@1:   methodTVo::LoadFromFileHandle: reads TVo object from a file stream.
xue@1: */
xue@1: void TVo::LoadFromFileHandle(FILE* file)
xue@1: {
xue@1:   fread(&M, sizeof(int), 1, file);
xue@1:   fread(&L, sizeof(int), 1, file);
xue@1:   fread(&P, sizeof(int), 1, file);
xue@1:   fread(&K, sizeof(int), 1, file);
xue@1:   fread(&h, sizeof(double), 1, file);
xue@1:   fread(&afres, sizeof(int), 1, file);
xue@1:   AllocateL(L); AllocateP(); AllocatePK();
xue@1:   fread(F0C, sizeof(double), L, file);
xue@1:   fread(A0C, sizeof(double), L, file);
xue@1:   fread(LogAF, sizeof(double), 4096, file);
xue@1:   fread(peakfr, sizeof(int), P, file);
xue@1:   fread(lp, sizeof(double), P, file);
xue@1:   fread(Dp, sizeof(double), P, file);
xue@1:   fread(Kp, sizeof(int), P, file);
xue@1:   fread(fxr[0], sizeof(double), P*(2*K+1), file);
xue@1:   fread(LogASp[0], sizeof(double), P*M, file);
xue@1: }//LoadFromFileHandle
xue@1: 
Chris@5: /**
xue@1:   method TVo::ReAllocateL: reallocates storage space whose size depends on L, and transfer the contents.
xue@1:   This uses an external L value for allocation but does not update L itself.
xue@1: */
xue@1: void TVo::ReAllocateL(int newL)
xue@1: {
xue@1:   double* newF0=new double[(newL+2)*4];
xue@1:   F0C=&newF0[newL+2], F0D=&newF0[(newL+2)*2], A0C=&newF0[(newL+2)*3];
xue@1:   memcpy(A0C, &F0[(L+2)*3], sizeof(double)*(L+2));
xue@1:   memcpy(F0D, &F0[(L+2)*2], sizeof(double)*(L+2));
xue@1:   memcpy(F0C, &F0[L+2], sizeof(double)*(L+2));
xue@1:   memcpy(newF0, F0, sizeof(double)*(L+2));
xue@1:   delete[] F0;
xue@1:   F0=newF0;
xue@1: }//ReAllocateL
xue@1: 
Chris@5: /**
xue@1:   method TVo::SaveToFile: saves a TVo object to a file.
xue@1: 
xue@1:   In: filename: path to the destination file.
xue@1: */
xue@1: void TVo::SaveToFile(char* filename)
xue@1: {
xue@1:   FILE* file=fopen(filename, "wb");
xue@1:   if (file)
xue@1:   {
xue@1:     SaveToFileHandle(file);
xue@1:     fclose(file);
xue@1:   }
xue@1:   else
xue@1:   {
xue@1:     throw(strcat((char*)"Cannot open file ", filename));
xue@1:   }
xue@1: }//SaveToFile
xue@1: 
Chris@5: /**
xue@1:   method TVo::SaveToFileHandle: saves a TVo object to an open file stream.
xue@1: */
xue@1: void TVo::SaveToFileHandle(FILE* file)
xue@1: {
xue@1:   fwrite(&M, sizeof(int), 1, file);
xue@1:   fwrite(&L, sizeof(int), 1, file);
xue@1:   fwrite(&P, sizeof(int), 1, file);
xue@1:   fwrite(&K, sizeof(int), 1, file);
xue@1: 	fwrite(&h, sizeof(double), 1, file);
xue@1:   fwrite(&afres, sizeof(int), 1, file);
xue@1: 	fwrite(F0C, sizeof(double), L, file);
xue@1:   fwrite(A0C, sizeof(double), L, file);
xue@1: 	fwrite(LogAF, sizeof(double), 4096, file);
xue@1:   fwrite(peakfr, sizeof(int), P, file);
xue@1:   fwrite(lp, sizeof(double), P, file);
xue@1: 	fwrite(Dp, sizeof(double), P, file);
xue@1:   fwrite(Kp, sizeof(int), P, file);
xue@1:   fwrite(fxr[0], sizeof(double), P*(2*K+1), file);
xue@1: 	fwrite(LogASp[0], sizeof(double), P*M, file);
xue@1: }//SaveToFileHandle
xue@1: 
xue@1: //---------------------------------------------------------------------------
xue@1: //functions
xue@1: 
xue@1: 
Chris@5: /**
xue@1:   function AnalyzeV: calculates all basic and non-basic descriptors of a vibrato from a harmonic
xue@1:   sinusoid
xue@1: 
xue@1:   In: HS: harmonic sinusoid to compute vibrato representation from
xue@1:       sps: sampling rate
xue@1:       h: hop size
xue@1:       SFMode: specifies source-filter estimation criterion, 0=FB (filter bank), 1=SV (slow variation)
xue@1:       SFF: filter response sampling interval
xue@1:       SFScale: set if filter sampling uses mel scale
xue@1:       SFtheta: balancing factor of amplitude and frequency variations, needed for SV approach
xue@1:   Out: V: a TVo object hosting vibrato descriptors (vibrato representation) of HS
xue@1:        cyclefrcount: number of cycles between cycle boundaries, equals V.P-1.
xue@1:        cyclefrs[cyclefrcount], cyclefs[cyclefrcount]: time (in frames) and F0 centroid of cycles. These are
xue@1:          knots from which V.F0C[] is interpolated.
xue@1:        peaky[V.P]: shape score of cycle peaks
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void AnalyzeV(THS& HS, TVo& V, double*& peaky, double*& cyclefrs, double*& cyclefs, double sps, double h, int* cyclefrcount,
xue@1:       int SFMode, double SFF, int SFFScale, double SFtheta)
xue@1: {
xue@1:   int M=HS.M, Fr=HS.Fr;
xue@1:   if (M<=0 || Fr<=0) return;
xue@1:   atom** Partials=HS.Partials;
xue@1: 
xue@1:   //basic descriptor: h
xue@1:   V.h=h;
xue@1: 
xue@1:   //demodulating pitch
xue@1:   //Basic descriptor: L
xue@1:   V.AllocateL(Fr);
xue@1: 
xue@1:   double* AA=new double[Fr];
xue@1:   //find the pitch track
xue@1:   V.F0max=0, V.F0min=1;
xue@1:   for (int fr=0; fr<Fr; fr++)
xue@1:   {
xue@1:     double suma2=0, suma2p=0;
xue@1:     for (int m=0; m<M; m++)
xue@1:     {
xue@1:       double a=Partials[m][fr].a, p=Partials[m][fr].f/(m+1);
xue@1:       if (p<=0) break;
xue@1:       suma2+=a*a, suma2p+=a*a*p;
xue@1:     }
xue@1:     V.F0[fr]=suma2p/suma2;
xue@1:     if (V.F0max<V.F0[fr]) V.F0max=V.F0[fr];
xue@1:     if (V.F0min>V.F0[fr]) V.F0min=V.F0[fr];
xue@1:     AA[fr]=suma2;
xue@1:   }
xue@1: 
xue@1:   //Basic descriptor: M
xue@1:   V.M=0.45/V.F0max;
xue@1:   if (V.M>M) V.M=M;
xue@1: 
xue@1:   //rough estimation of rate
xue@1:   double* f0d=new double[Fr-1]; for (int i=0; i<Fr-1; i++) f0d[i]=V.F0[i+1]-V.F0[i];
xue@1:   RateAndReg(V.rate, V.regularity, f0d, 0, Fr-2, 8, sps, h);
xue@1:   delete[] f0d;
xue@1: 
xue@1:   //find peaks of the pitch track
xue@1:   free(V.peakfr); V.peakfr=(int*)malloc(sizeof(int)*Fr/2);
xue@1:   double periodinframe=sps/(V.rate*h);
xue@1:   FindPeaks(V.peakfr, V.P, V.F0, Fr, periodinframe, peaky);
xue@1:   if (V.P>=1) //de-modulation
xue@1:   {
xue@1:     //Basic descriptor: F0C[]
xue@1:     DeFM(V.F0C, V.F0, AA, V.P, V.peakfr, Fr, V.F0Overall, *cyclefrcount, cyclefrs, cyclefs, true);
xue@1:     //Basic descriptor: A0C[]
xue@1:     DeAM(V.A0C, AA, V.P, V.peakfr, Fr);
xue@1:   }
xue@1: 
xue@1:   V.F0Cmax=0; V.F0Dmax=-1; V.F0Cmin=1; V.F0Dmin=1;
xue@1:   for (int fr=0; fr<Fr; fr++)
xue@1:   {
xue@1:     if (V.F0Cmax<V.F0C[fr]) V.F0Cmax=V.F0C[fr];
xue@1:     if (V.F0Cmin>V.F0C[fr]) V.F0Cmin=V.F0C[fr];
xue@1:     V.F0D[fr]=V.F0[fr]-V.F0C[fr];
xue@1:   }
xue@1: 
xue@1:   V.D=0; int cyclestart=(V.P>=2)?V.peakfr[0]:0, cycleend=(V.P>=2)?V.peakfr[V.P-1]:Fr;
xue@1:   for (int fr=cyclestart; fr<cycleend; fr++)
xue@1:   {
xue@1:     double lf0d=V.F0D[fr];
xue@1:     if (V.F0Dmax<lf0d) V.F0Dmax=lf0d;
xue@1:     if (V.F0Dmin>lf0d) V.F0Dmin=lf0d;
xue@1:     V.D+=lf0d*lf0d;
xue@1:   }
xue@1:   V.D=sqrt(V.D/(cycleend-cyclestart));
xue@1:   delete[] AA;
xue@1: 
xue@1:   //Calculate rate and regularity
xue@1:   //find the section used for calculating regularity and rate
xue@1:   {
xue@1:     int frst=0, fren=Fr-1;
xue@1:     while (frst<fren && V.F0D[frst]*V.F0D[frst+1]>0) frst++;
xue@1:     while (fren>frst && V.F0D[fren]*V.F0D[fren-1]>0) fren--;
xue@1: 
xue@1:     RateAndReg(V.rate, V.regularity, V.F0D, frst, fren, 8, sps, h);
xue@1:   }
xue@1: 
xue@1:   //Basic descriptor: P
xue@1:   //Basic descriptor: peakfr[]
xue@1:   periodinframe=sps/(V.rate*h);
xue@1:   FindPeaks(V.peakfr, V.P, V.F0D, Fr, periodinframe, peaky);
xue@1: 
xue@1:   //Basic descriptor: lp[]
xue@1:   V.AllocateP();
xue@1:   for (int p=0; p<V.P; p++)
xue@1:   {
xue@1:     double startfr; QIE(&V.F0D[V.peakfr[p]], startfr); V.lp[p]=startfr+V.peakfr[p];
xue@1:   }
xue@1: 
xue@1:   //Basic descriptor: LogAF[]
xue@1:   //Source-filter modeling
xue@1:   if (SFMode==0)
xue@1:     S_F(V.afres, V.LogAF, V.LogAS, Fr, M, Partials, V.A0C, V.lp, V.P, V.F0Overall);
xue@1:   else
xue@1:     S_F_SV(V.afres, V.LogAF, V.LogAS, Fr, M, Partials, V.A0C, V.lp, V.P, SFF, SFFScale, SFtheta, sps);
xue@1: 
xue@1: 
xue@1:   //the alternative: find K only
xue@1:   V.K=50;
xue@1:   if (V.K>periodinframe/2) V.K=periodinframe/2;
xue@1: 
xue@1:   //single-cycle analyses
xue@1:   //Basic descriptor: Dp[]
xue@1:   V.AllocatePK();
xue@1:   for (int p=1; p<V.P; p++)
xue@1:   {
xue@1:     //int frst=ceil(V.lp[p-1]), freninc=floor(V.lp[p]);
xue@1:     int frst=V.peakfr[p-1], freninc=V.peakfr[p];
xue@1:     V.Dp[p]=0; for (int fr=frst; fr<=freninc; fr++) V.Dp[p]+=V.F0D[fr]*V.F0D[fr]; V.Dp[p]=sqrt(V.Dp[p]/(freninc-frst+1));
xue@1:     double* f0e=new double[freninc-frst+1];
xue@1:     for (int fr=0; fr<=freninc-frst; fr++) f0e[fr]=V.F0D[frst+fr]/V.Dp[p];
xue@1:     V.Kp[p]=V.K; FS_QR(V.Kp[p], V.fxr[p], &f0e[0], freninc-frst, V.lp[p]-V.lp[p-1], V.lp[p-1]-frst, V.fres[p]);
xue@1:     delete[] f0e;
xue@1: 
xue@1:     memset(V.LogASp[p], 0, sizeof(double)*V.M);
xue@1:     for (int m=0; m<V.M; m++)
xue@1:     {
xue@1:       int ccount=0;
xue@1:       for (int fr=frst; fr<=freninc; fr++)
xue@1:       {
xue@1:         double f=Partials[m][fr].f*V.afres;
xue@1:         if (f<=0) continue;
xue@1:         int intf=floor(f);
xue@1:         V.LogASp[p][m]=V.LogASp[p][m]+log(Partials[m][fr].a/V.A0C[fr])-(V.LogAF[intf]*(intf+1-f)+V.LogAF[intf+1]*(f-intf));
xue@1:         ccount++;
xue@1:       }
xue@1:       if (ccount>0) V.LogASp[p][m]/=ccount;
xue@1:       else V.LogASp[p][m]=0;
xue@1:     }
xue@1:   }
xue@1:   memset(V.FXR, 0, sizeof(double)*2*V.K); memset(V.FRes, 0, sizeof(double)*V.K);
xue@1:   double sumdp=0; for (int p=1; p<V.P; p++) sumdp+=V.Dp[p], V.FXR[0]+=V.fxr[p][0]*V.Dp[p], V.FRes[0]+=V.fres[p][0]*V.Dp[p]; V.FXR[0]/=sumdp, V.FRes[0]/=sumdp;
xue@1:   for (int k=1; k<V.K; k++)
xue@1:   {
xue@1:     double sumdp=0;
xue@1:     for (int p=1; p<V.P; p++)
xue@1:     {
xue@1:       if (k<V.Kp[p])
xue@1:       {
xue@1:         V.FXR[2*k-1]+=V.fxr[p][2*k-1]*V.Dp[p];
xue@1:         V.FXR[2*k]+=V.fxr[p][2*k]*V.Dp[p];
xue@1:         V.FRes[k]+=V.fres[p][k]*V.Dp[p];
xue@1:         sumdp+=V.Dp[p];
xue@1:       }
xue@1:     }
xue@1:     V.FXR[2*k-1]/=sumdp, V.FXR[2*k]/=sumdp, V.FRes[k]/=sumdp;
xue@1:   }
xue@1: }//AnalyzeV
xue@1: 
Chris@5: /**
xue@1:   function copeak: measures the similarity of F0 between frst and fren to a cosine cycle 0~2pi.
xue@1: 
xue@1:   In: F0[peakfr[frst]:peakfr[fren]]
xue@1:       periodinframe: period of the cosine reference, in frames
xue@1: 
xue@1:   Returns the correlation coefficient.
xue@1: */
xue@1: double copeak(double* F0, int* peakfr, int frst, int fren, double periodinframe)
xue@1: {
xue@1:   double ef0=0, rfs=0, rff=0, rss=0; int frcount=0;
xue@1:   for (int fr=peakfr[frst]; fr<=peakfr[fren]; fr++) ef0+=F0[fr], frcount++; ef0/=frcount;
xue@1:   for (int fr=peakfr[frst]; fr<=peakfr[fren]; fr++)
xue@1:   {
xue@1:     double s=cos(2*M_PI*(fr-peakfr[frst])/periodinframe), f=F0[fr];
xue@1:     rfs+=(f-ef0)*s; rff+=(f-ef0)*(f-ef0); rss+=s*s;
xue@1:   }
xue@1:   return rfs/sqrt(rff*rss);
xue@1: }//copeak
xue@1: 
Chris@5: /**
xue@1:   function CorrectFS: time-shifts Fourier series so that maximum is reached at 0. This is done by
xue@1:   finding the actual peak numerically then shifting it to 0.
xue@1: 
xue@1:   In: c[2K-1]: Fourier series coefficients, in the order of 0r, 1i, 1r, 2i, 2r, ..., so that the series
xue@1:         is expanded as c[0]+c[1]sinx+c[2]cosx+c[3]sin2x+c[4]cos2x+...c[2K-2]cos(K-1)x
xue@1:       epx: a small positive value for convergence test for finding the peak
xue@1:       maxiter: maximal number of iterates for finding the peak
xue@1:   Out: c[2K-1]: updated Fourier series coefficients.
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void CorrectFS(int K, double* c, double epx=1e-8, int maxiter=10)
xue@1: {
xue@1:   double y=0, kb=0, kka=0, th2pi=0, dth2pi=0, my, mth2pi;
xue@1:   epx=2*M_PI*epx;
xue@1:   int iter=0;
xue@1:   for (int k=1; k<K; k++) y+=c[2*k], kb+=k*c[2*k-1], kka+=k*k*c[2*k]; my=y;
xue@1:   dth2pi=(kka>0)?(kb/kka):(kb*0.1);
xue@1:   while (iter<maxiter && fabs(dth2pi)>epx)
xue@1:   {
xue@1:     iter++;
xue@1:     th2pi+=dth2pi;
xue@1:     kb=kka=y=0;
xue@1:     for (int k=1; k<K; k++)
xue@1:     {
xue@1:       double cc=cos(k*th2pi), ss=sin(k*th2pi);
xue@1:       y+=(c[2*k]*cc+c[2*k-1]*ss);
xue@1:       kb+=k*(-c[2*k]*ss+c[2*k-1]*cc);
xue@1:       kka+=k*k*(c[2*k]*cc+c[2*k-1]*ss);
xue@1:     }
xue@1:     if (y>my) my=y, mth2pi=th2pi;
xue@1:     dth2pi=(kka>0)?(kb/kka):(kb*0.1);
xue@1:   }
xue@1:   if (iter>=maxiter || th2pi!=mth2pi)
xue@1:   {
xue@1:     th2pi=mth2pi;
xue@1:   }
xue@1:   for (int k=1; k<K; k++)
xue@1:   {
xue@1:     double cc=cos(k*th2pi), ss=sin(k*th2pi);
xue@1:     double tmp=c[2*k]*cc+c[2*k-1]*ss;
xue@1:     c[2*k-1]=-c[2*k]*ss+c[2*k-1]*cc; c[2*k]=tmp;
xue@1:   }
xue@1: }//CorrectFS
xue@1: 
Chris@5: /**
xue@1:   function DeAM: amplitude demodulation based on given segmentation into cycles
xue@1: 
xue@1:   In: A1[Fr]: original amplitude
xue@1:       peakfr[npfr]: cycle boundaries (literally, "frame indices of frequency modulation peaks")
xue@1:   Out: A2[Fr]: demodulated amplitude
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void DeAM(double* A2, double* A1, int npfr, int* peakfr, int Fr)
xue@1: {
xue@1:   if (npfr==1) {memcpy(A2, A1, sizeof(double)*Fr); return;}
xue@1:   double *frs=new double[npfr*2], *a=&frs[npfr];
xue@1:   for (int i=0; i<npfr-1; i++)
xue@1:   {
xue@1:     int frst=peakfr[i], fren=peakfr[i+1];
xue@1:     a[i]=(A1[frst]+A1[fren])*0.5;
xue@1:     frs[i]=(frst*A1[frst]+fren*A1[fren])*0.5;
xue@1:     for (int fr=frst+1; fr<fren; fr++) a[i]+=A1[fr], frs[i]+=fr*A1[fr];
xue@1:     frs[i]/=a[i]; a[i]=log(sqrt(a[i]/(fren-frst)));
xue@1:   }
xue@1:   if (npfr==2)
xue@1:   {
xue@1:     A2[0]=exp(a[0]);
xue@1:     for (int fr=1; fr<Fr; fr++) A2[fr]=A2[0];
xue@1:   }
xue@1:   else if (npfr==3)
xue@1:   {
xue@1:     double alp=(a[1]-a[0])/(frs[1]-frs[0]);
xue@1:     for (int fr=0; fr<Fr; fr++) A2[fr]=exp(a[0]+(fr-frs[0])*alp);
xue@1:   }
xue@1:   else if (npfr==4)
xue@1:   {
xue@1:     double x0=frs[0], x1=frs[1], x2=frs[2], y0=log(a[0]), y1=log(a[1]), y2=log(a[2]);
xue@1:     double ix0_12=y0/((x0-x1)*(x0-x2)), ix1_02=y1/((x1-x0)*(x1-x2)), ix2_01=y2/((x2-x0)*(x2-x1));
xue@1:     double aa=ix0_12+ix1_02+ix2_01,
xue@1:            b=-(x1+x2)*ix0_12-(x0+x2)*ix1_02-(x0+x1)*ix2_01,
xue@1:            c=x1*x2*ix0_12+x0*x2*ix1_02+x0*x1*ix2_01;
xue@1:     for (int fr=0; fr<Fr; fr++) A2[fr]=exp(c+fr*(b+fr*aa));
xue@1:   }
xue@1:   else
xue@1:   {
xue@1:     double *fa=new double[npfr*4], *fb=&fa[npfr], *fc=&fa[npfr*2], *fd=&fa[npfr*3];
xue@1:     CubicSpline(npfr-2, fa, fb, fc, fd, frs, a, 0, 1, A2, 1, 0, Fr-1);
xue@1:     for (int fr=0; fr<Fr; fr++) A2[fr]=exp(A2[fr]);
xue@1:     delete[] fa;
xue@1:   }
xue@1:   delete[] frs;
xue@1: }//DeAM
xue@1: 
Chris@5: /**
xue@1:   function DeAM: wrapper function using floating-point cycle boundaries
xue@1: 
xue@1:   In: A1[Fr]: original amplitude
xue@1:       lp[npfr]: cycle boundaries
xue@1:   Out: A2[Fr]: demodulated amplitude
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void DeAM(double* A2, double* A1, int npfr, double* lp, int Fr)
xue@1: {
xue@11:   int* peakfr=(int*)malloc(sizeof(int)*npfr);
xue@1:   for (int i=0; i<npfr; i++) peakfr[i]=ceil(lp[i]);
xue@1:   DeAM(A2, A1, npfr, peakfr, Fr);
xue@11:   free(peakfr);
xue@1: }//DeAM
xue@1: 
Chris@5: /**
xue@1:   function DeFM: frequency demodulation based on given segmentation into cycles
xue@1: 
xue@1:   In: f1[Fr], AA[Fr]: original frequency and partial-independent amplitude
xue@1:       peakfr[npfr]: cycle boundaries (literally, "frame indices of frequency modulation peaks")
xue@1:       furthursmoothing: specifies if a second round demodulation is to be performed for more smooth
xue@1:         carrier
xue@1:   Out: f2[Fr]: demodulated frequency
xue@1:        f[fcount], frs[fcount]: frequency and time (in frames) of carrier knots, optional
xue@1:        fct: frequency centroid
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void DeFM(double* f2, double* f1, double* AA, int npfr, int* peakfr, int Fr, double& fct, int& fcount, double* &frs, double* &f, bool furthersmoothing)
xue@1: {
xue@1:   double ac=0;
xue@1:   fct=0;
xue@1: 
xue@1:   if (npfr==1)
xue@1:   {
xue@1:     memcpy(f2, f1, sizeof(double)*Fr);
xue@1:     for (int fr=0; fr<Fr; fr++) ac+=AA[fr], fct+=f1[fr]*AA[fr]; fct/=ac;
xue@1:     return;
xue@1:   }
xue@1: 
xue@1:   bool localfrs=(frs==(double*)0xFFFFFFFF), localf=(f==(double*)0xFFFFFFFF);
xue@1:   if (!localfrs) delete[] frs;
xue@1:   if (!localf) delete[] f;
xue@1:   frs=new double[npfr]; f=new double[npfr];
xue@1:   for (int i=0; i<npfr-1; i++)
xue@1:   {
xue@1:     int frst=peakfr[i], fren=peakfr[i+1];
xue@1:     f[i]=(f1[frst]*AA[frst]+f1[fren]*AA[fren])*0.5;
xue@1:     frs[i]=(frst*AA[frst]+fren*AA[fren])*0.5;
xue@1:     double lrec=(AA[frst]+AA[fren])*0.5;
xue@1:     for (int fr=frst+1; fr<fren; fr++) f[i]+=f1[fr]*AA[fr], frs[i]+=fr*AA[fr], lrec+=AA[fr];
xue@1:     ac+=lrec;
xue@1:     fct+=f[i];
xue@1:     f[i]/=lrec, frs[i]/=lrec;
xue@1:   }
xue@1:   fct/=ac;
xue@1: 
xue@1:   if (furthersmoothing)
xue@1:   {
xue@1:     //further smoothing
xue@1:     int newnpfr=1;
xue@1:     for (int pfr=1; pfr<npfr-2; pfr++)
xue@1:     {
xue@1:       if ((f[pfr]>f[pfr-1] && f[pfr]>f[pfr+1]) || (f[pfr]<f[pfr-1] && f[pfr]<f[pfr+1]))
xue@1:       {
xue@1:         //this is a peak or valley
xue@1:         if ((pfr+1<npfr-2 && f[pfr+1]>f[pfr] && f[pfr+1]>f[pfr+2]) || (f[pfr+1]<f[pfr] && f[pfr+1]<f[pfr+2]))
xue@1:         {
xue@1:           frs[newnpfr]=0.5*(frs[pfr]+frs[pfr+1]), f[newnpfr]=0.5*(f[pfr]+f[pfr+1]);
xue@1:           newnpfr++;
xue@1:         }
xue@1:       }
xue@1:       else
xue@1:       {
xue@1:         frs[newnpfr]=frs[pfr], f[newnpfr]=f[pfr];
xue@1:         newnpfr++;
xue@1:       }
xue@1:     }
xue@1:     frs[newnpfr]=frs[npfr-2], f[newnpfr]=f[npfr-2]; newnpfr++; newnpfr++;
xue@1:     npfr=newnpfr;
xue@1:   }
xue@1: 
xue@1:   Smooth_Interpolate(f2, Fr, npfr, f, frs);
xue@1:   fcount=npfr-1;
xue@1: 
xue@1:   if (localfrs) delete[] frs;
xue@1:   if (localf) delete[] f;
xue@1: }//DeFM
xue@1: 
Chris@5: /**
xue@1:   function DeFM: wrapper function using floating-point cycle boundaries
xue@1: 
xue@1:   In: f1[Fr], AA[Fr]: original frequency and partial-independent amplitude
xue@1:       lp[npfr]: cycle boundaries
xue@1:       furthursmoothing: specifies if a second round demodulation is to be performed for more smooth
xue@1:         carrier
xue@1:   Out: f2[Fr]: demodulated frequency
xue@1:        f[fcount], frs[fcount]: frequency and time (in frames) of carrier knots, optional
xue@1:        fct: frequency centroid
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void DeFM(double* f2, double* f1, double* AA, int npfr, double* lp, int Fr, double& fct, int& fcount, double* &frs, double* &f, bool furthersmoothing)
xue@1: {
xue@11:   int* peakfr=(int*)malloc(sizeof(int)*npfr);
xue@1:   for (int i=0; i<npfr; i++) peakfr[i]=ceil(lp[i]);
xue@1:   DeFM(f2, f1, AA, npfr, peakfr, Fr, fct, fcount, frs, f, furthersmoothing);
xue@11:   free(peakfr);
xue@1: }//DeFM
xue@1: 
Chris@5: /**
xue@1:   function FindPeaks: find modulation peaks of signal F0[] roughly periodical at $periodinframe
xue@1: 
xue@1:   In: F0[Fr]: modulated signal
xue@1:       periodinframe: reference modulation period
xue@1:   Out: npfr, peakfr[npfr]: modulation peaks
xue@1:        peaky[npfr]: shape score for the peaks measuring their similarity to cosine peak, optional
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void FindPeaks(int* peakfr, int& npfr, double* F0, int Fr, double periodinframe, double*& peaky)
xue@1: {
xue@1:   npfr=0;
xue@1:   for (int fr=1; fr<Fr-1; fr++)
xue@1:   {
xue@1:     if (F0[fr]>F0[fr-1] && F0[fr]>F0[fr+1])
xue@1:     {
xue@1:       peakfr[npfr++]=fr;
xue@1:     }
xue@1:   }
xue@1: 
xue@1:   bool localpeaky=(peaky==(double*)0xFFFFFFFF); if (!localpeaky) delete[] peaky;
xue@1: 
xue@1:   peaky=new double[npfr];
xue@1:   for (int pfr=0; pfr<npfr; pfr++)
xue@1:   {
xue@1:     double rfs=0, rff=0, rss=0;
xue@1:     if (peakfr[pfr]>periodinframe/2)
xue@1:     {
xue@1:       double ef0=0; int frcount=0; for (int fr=peakfr[pfr]-periodinframe/2; fr<peakfr[pfr]; fr++) ef0+=F0[fr], frcount++; ef0/=frcount;
xue@1:       for (int fr=peakfr[pfr]-periodinframe/2; fr<peakfr[pfr]; fr++)
xue@1:       {
xue@1:         double s=cos(2*M_PI*(fr-peakfr[pfr])/periodinframe), f=F0[fr];
xue@1:         rfs+=(f-ef0)*s; rff+=(f-ef0)*(f-ef0); rss+=s*s;
xue@1:       }
xue@1:     }
xue@1:     if (peakfr[pfr]<Fr-periodinframe/2)
xue@1:     {
xue@1:       double ef0=0; int frcount=0; for (int fr=peakfr[pfr]; fr<peakfr[pfr]+periodinframe/2; fr++) ef0+=F0[fr], frcount++; ef0/=frcount;
xue@1:       for (int fr=peakfr[pfr]; fr<peakfr[pfr]+periodinframe/2; fr++)
xue@1:       {
xue@1:         double s=cos(2*M_PI*(fr-peakfr[pfr])/periodinframe), f=F0[fr];
xue@1:         rfs+=(f-ef0)*s; rff+=(f-ef0)*(f-ef0); rss+=s*s;
xue@1:       }
xue@1:     }
xue@1:     peaky[pfr]=rfs/sqrt(rff*rss);
xue@1:   }
xue@1:   int pst=0; for (int pfr=1; pfr<npfr; pfr++) if (peaky[pst]<peaky[pfr]) pst=pfr;
xue@1: //    return;
xue@1:   //*
xue@1:   int pen=pst+1, pfr=pst+1;
xue@1:   while (pfr<npfr)
xue@1:   {
xue@1:     int impfr=-1;
xue@1:     while (pfr<npfr && peakfr[pfr]<peakfr[pen-1]+periodinframe*0.5) pfr++;
xue@1:     if (pfr<npfr && peakfr[pfr]<peakfr[pen-1]+periodinframe*1.5)
xue@1:     {
xue@1:       double mpfr=peaky[pfr]+copeak(F0, peakfr, pen-1, pfr, periodinframe);
xue@1:       impfr=pfr; pfr++;
xue@1:       while (pfr<npfr && peakfr[pfr]<peakfr[pen-1]+periodinframe*1.5)
xue@1:       {
xue@1:         double mp=peaky[pfr]+copeak(F0, peakfr, pen-1, pfr, periodinframe);
xue@1:         if (mp>mpfr) impfr=pfr, mpfr=mp;
xue@1:         pfr++;
xue@1:       }
xue@1:     }
xue@1:     else
xue@1:     {
xue@1:       while (pfr<npfr)
xue@1:       {
xue@1:         if (peaky[pfr]>0) {impfr=pfr; break;}
xue@1:         pfr++;
xue@1:       }
xue@1:     }
xue@1:     if (impfr>=0) peakfr[pen++]=peakfr[impfr];
xue@1:   }
xue@1:   pst--; pfr=pst;
xue@1:   while (pfr>=0)
xue@1:   {
xue@1:     int impfr=-1;
xue@1:     while (pfr>=0 && peakfr[pfr]>peakfr[pst+1]-periodinframe*0.5) pfr--;
xue@1:     if (pfr>=0 && peakfr[pfr]>=peakfr[pst+1]-periodinframe*1.5)
xue@1:     {
xue@1:       double mpfr=peaky[pfr]+copeak(F0, peakfr, pfr, pst+1, periodinframe);
xue@1:       impfr=pfr; pfr--;
xue@1:       while (pfr>=0 && peakfr[pfr]>=peakfr[pst+1]-periodinframe*1.5)
xue@1:       {
xue@1:         double mp=peaky[pfr]+copeak(F0, peakfr, pfr, pst+1, periodinframe);
xue@1:         if (mp>mpfr) impfr=pfr, mpfr=mp;
xue@1:         pfr--;
xue@1:       }
xue@1:     }
xue@1:     else
xue@1:     {
xue@1:       while (pfr>=0)
xue@1:       {
xue@1:         if (peaky[pfr]>0) {impfr=pfr; break;}
xue@1:         pfr--;
xue@1:       }
xue@1:     }
xue@1:     if (impfr>=0) peakfr[pst--]=peakfr[impfr];
xue@1:   }
xue@1:   pst++;
xue@1:   npfr=pen-pst;
xue@1:   memmove(peakfr, &peakfr[pst], sizeof(int)*npfr); //*/
xue@1:   if (localpeaky) delete[] peaky;
xue@1: }//FindPeaks
xue@1: 
Chris@5: /**
xue@1:   function FS: Fourier series decomposition
xue@1: 
xue@1:   In: data[frst:fren]: signal to decompose
xue@1:       period: period of Fourier series decomposition
xue@1:       shift: amount of original shift of Fourier components (from frst to frst-shift)
xue@1:       K: number of Fourier components requested
xue@1:   Out: K: number of Fourier components returned
xue@1:        FXR[K], FXI[K]: real and imaginary parts of Fourier series coefficients
xue@1:        FRES[K]: decomposition residues
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void FS(int& K, double* FXR, double* FXI, double* FRES, double* data, int frst, int fren, double period, double shift)
xue@1: {
xue@1:   if (K>period/2) K=period/2;
xue@1:   int len=fren-frst+1;
xue@1:   double omg0=2*M_PI/period, ilen=1.0/len;
xue@1: 
xue@1:   //for (int fr=frst; fr<=fren; fr++) data[fr]=cos(2*M_PI*fr/period+1);
xue@1: 
xue@1:   double* res=new double[len]; memcpy(res, &data[frst], sizeof(double)*len);
xue@1:   double* rec=new double[len];
xue@1:   double resene=0; for (int i=0; i<len; i++) resene+=res[i]*res[i]; double ene=resene;
xue@1: 
xue@1:   for (int k=0; k<K; k++)
xue@1:   {
xue@1: //      double eer=0, eei=0;
xue@1:     double xr=0, xi=0, omg=k*omg0;
xue@1:     for (int fr=frst; fr<=fren; fr++)
xue@1:     {
xue@1:       double ph=omg*(fr-shift);
xue@1:       double c=cos(ph), s=sin(ph);
xue@1:       xr+=data[fr]*c, xi+=data[fr]*s;
xue@1: //        eer+=c*c, eei+=s*s;
xue@1:     }
xue@1:     xr*=ilen, xi*=ilen;
xue@1: //      if (eer>0) xr/=eer;
xue@1: //      if (eei>0) xi/=eei;
xue@1:     FXR[k]=xr, FXI[k]=xi;
xue@1:     double lresene=0;
xue@1:     for (int fr=frst; fr<=fren; fr++)
xue@1:     {
xue@1:       double ph=omg*(fr-shift);
xue@1: //        rec[fr-frst]=xr*cos(ph)+xi*sin(ph);
xue@1:       if (k==0) rec[fr-frst]=xr*cos(ph)+xi*sin(ph);
xue@1:       else rec[fr-frst]=2*(xr*cos(ph)+xi*sin(ph));
xue@1:       res[fr-frst]-=rec[fr-frst];
xue@1:       lresene+=res[fr-frst]*res[fr-frst];
xue@1:     }
xue@1:     resene=lresene;
xue@1:     FRES[k]=resene/ene;
xue@1:   }
xue@1:   delete[] rec;
xue@1:   delete[] res;
xue@1: }//FS
xue@1: 
Chris@5: /**
xue@1:   function FS_QR: Fourier series decomposition with QR orthogonalization. Since the Fourier series is
xue@1:   applied on finite-length discrate signal, the Fourier components are no longer orthogonal to each
xue@1:   other. A decreasing residue can be guaranteed by applying QR (or any other) orthogonalization process
xue@1:   to the Fourier components before decomposition.
xue@1: 
xue@1:   In: data[0:Fr]: signal to decompose
xue@1:       period: period of Fourier series decomposition
xue@1:       shift: amount of original shift of Fourier components (from 0 to -shift)
xue@1:       K: number of Fourier components requested
xue@1:   Out: K: number of Fourier components returned
xue@1:        FXR[2K-1]: Fourier series coefficients, in the order of 0r, 1i, 1r, 2i, 2r, ....
xue@1:        FRES[K]: decomposition residues, optional
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void FS_QR(int& K, double* FXR, double* data, int Fr, double period, double shift, double* FRES)
xue@1: {
xue@1:   int len=Fr+1;
xue@1:   if (K>period/2) K=period/2;
xue@1:   if (K>Fr/2) K=Fr/2;
xue@1:   if (K>len/2) K=len/2;
xue@1:   if (K<=0) return;
xue@1:   memset(FXR, 0, sizeof(double)*(2*K-1));
xue@1: 
xue@1:   double omg0=2*M_PI/period, ene, resene;
xue@1: 
xue@1:   if (FRES)
xue@1:   {
xue@1:     ene=0; for (int i=0; i<len; i++) ene+=data[i]*data[i]; resene=ene;
xue@1:   }
xue@1: 
xue@1:   /*debug only
xue@1:   double *res=new double[len*4], *rec=&res[len];
xue@1:   memcpy(res, data, sizeof(double)*len); memset(rec, 0, sizeof(double)*len);
xue@1:   double resene2=resene;//*/
xue@1: 
xue@1:   int NK=K*2-1;
xue@1: 
xue@1:   Alloc2(len, NK, A); Alloc2(len, len, Q); Alloc2(len, NK, R);
xue@1:   for (int fr=0; fr<=Fr; fr++) A[fr][0]=1;
xue@1:   for (int k=1; k<(NK+1)/2; k++)
xue@1:   {
xue@1:     double omg=k*omg0;
xue@1:     for (int fr=0; fr<=Fr; fr++){double ph=omg*(fr-shift); A[fr][2*k-1]=2*sin(ph), A[fr][2*k]=2*cos(ph);}
xue@1:   }
xue@1:   QR_householder(len, NK, A, Q, R);///////////////////
xue@1:   GIUT(NK, R);
xue@1:   for (int k=0; k<NK; k++)
xue@1:   {
xue@1:     double xr=0;
xue@1:     for (int fr=0; fr<=Fr; fr++)
xue@1:     {
xue@1:       double c=Q[fr][k];
xue@1:       xr+=data[fr]*c;
xue@1:     }
xue@1:     for (int kk=0; kk<=k; kk++) FXR[kk]+=xr*R[kk][k];
xue@1:     /*debug only
xue@1:     double lresene=0;
xue@1:     for (int fr=0; fr<=Fr; fr++)
xue@1:     {
xue@1:       rec[fr]=0; for (int kk=0; kk<=k; kk++) rec[fr]+=FXR[kk]*A[fr][kk];
xue@1:       res[fr]=data[fr]-rec[fr];
xue@1:       lresene+=res[fr]*res[fr];
xue@1:     }
xue@1:     resene2=lresene;   */
xue@1: 
xue@1:     if (FRES)
xue@1:     {
xue@1:       resene-=xr*xr;
xue@1:       if (k%2==0) FRES[k/2]=resene/ene;
xue@1:     }
xue@1:   }
xue@1:   /* debug only
xue@1:   delete[] res; //*/
xue@1:   DeAlloc2(A); DeAlloc2(Q); DeAlloc2(R);
xue@1: }//FS_QR
xue@1: 
Chris@5: /**
xue@1:   function InterpolateV: adjusts modulation rate to $rate times the current value. Since TVo is based on
xue@1:   individual cycles, this operation involves finding new cycle boundaries and recomputing other single-
xue@1:   cycle descriptors by interpolation. This is used for time stretching and cycle rate adjustment.
xue@1: 
xue@1:   In: V: a TVo object to adjust
xue@1:       newP: number of total cycles after adjustment.
xue@1:       rate: amount of adjustment.
xue@1:   Out: another TVo object with new cycle boundaries and single-cycle descriptors interpoated from V
xue@1: 
xue@1:   Returns pointer to the output TVo. This function does not affect the content in V.
xue@1: */
xue@1: TVo* InterpolateV(double newP, double rate, TVo& V)
xue@1: {
xue@1:   double Pst=-V.lp[0]/(V.lp[1]-V.lp[0]);
xue@1:   double oldP=V.P-1+(V.L-V.lp[V.P-1])/(V.lp[V.P-1]-V.lp[V.P-2])-Pst;
xue@1: 
xue@1:   TVo* newV=new TVo;
xue@1:   memcpy(newV, &V, sizeof(TVo));
xue@1: 
xue@1:   newV->F0C=0; newV->A0C=0; newV->peakfr=0; newV->Dp=0; newV->lp=0, newV->Kp=0;
xue@1:   newV->fxr=0; newV->LogASp=0; newV->F0=0; newV->F0D=0; newV->fres=0;
xue@1: 
xue@1:   int inewP=ceil(newP+Pst);
xue@1:   newV->P=inewP;
xue@1: 
xue@11:   newV->peakfr=(int*)malloc(sizeof(int)*inewP);
xue@1:   newV->AllocateP(); newV->AllocatePK();
xue@1: 
xue@1:   newV->lp[0]=V.lp[0]/rate;
xue@1:   newV->peakfr[0]=floor(newV->lp[0]+0.5);
xue@1:   for (int p=1; p<inewP; p++)
xue@1:   {
xue@1:     double rp=1+(p-1)*(oldP-2+Pst)/(newP-2+Pst);
xue@1:     int ip=floor(rp); rp=rp-ip;
xue@1:     if (ip+1>=V.P) ip=V.P-1, rp=0;
xue@1:     if (fabs(rp)<1e-16)
xue@1:     {
xue@1:       newV->lp[p]=newV->lp[p-1]+(V.lp[ip]-V.lp[ip-1])/rate;
xue@1:       newV->Dp[p]=V.Dp[ip];
xue@1:       newV->Kp[p]=V.Kp[ip];
xue@1:       memcpy(newV->fxr[p], V.fxr[ip], sizeof(double)*2*newV->Kp[p]);
xue@1:       memcpy(newV->LogASp[p], V.LogASp[ip], sizeof(double)*V.M);
xue@1:     }
xue@1:     else
xue@1:     {
xue@1:       double fv1=1.0/(V.lp[ip]-V.lp[ip-1]), fv2=1.0/(V.lp[ip+1]-V.lp[ip]);
xue@1:       double fv=rate*(fv1*(1-rp)+fv2*rp);
xue@1:       newV->lp[p]=newV->lp[p-1]+1.0/fv;
xue@1:       newV->Dp[p]=V.Dp[ip]*(1-rp)+V.Dp[ip+1]*rp;
xue@1:       int minK=V.Kp[ip]; if (minK>V.Kp[ip+1]) minK=V.Kp[ip+1];
xue@1:       for (int k=0; k<2*minK-1; k++) newV->fxr[p][k]=V.fxr[ip][k]*(1-rp)+V.fxr[ip+1][k]*rp;
xue@1:       for (int m=0; m<V.M; m++) newV->LogASp[p][m]=V.LogASp[ip][m]*(1-rp)+V.LogASp[ip+1][m]*rp;
xue@1:       newV->Kp[p]=minK;
xue@1:     }
xue@1:     newV->peakfr[p]=floor(newV->lp[p]+0.5);
xue@1:     memset(newV->fres[0], 0, sizeof(double)*newV->P*V.K);
xue@1:   }
xue@1:   int newL=newV->lp[inewP-1]+(newP+Pst-(inewP-1))*(V.lp[V.P-1]-V.lp[V.P-2])/rate+0.5;
xue@1:   newV->AllocateL(newL);
xue@1: 
xue@1:   //F0C and A0C
xue@1:   if (V.L!=newL)
xue@1:   {
xue@1:     for (int l=0; l<newL; l++)
xue@1:     {
xue@1:       double rl=1.0*l*(V.L-1)/(newL-1);
xue@1:       int il=floor(rl); rl-=il;
xue@1:       newV->F0C[l]=V.F0C[il]*(1-rl)+V.F0C[il+1]*rl;
xue@1:       newV->A0C[l]=V.A0C[il]*(1-rl)+V.A0C[il+1]*rl;
xue@1:     }
xue@1:   }
xue@1:   else
xue@1:   {
xue@1:     memcpy(newV->F0C, V.F0C, sizeof(double)*V.L);
xue@1:     memcpy(newV->A0C, V.A0C, sizeof(double)*V.L);
xue@1: //    memcpy(newV->F0, V.F0, sizeof(double)*V.L);
xue@1: //    memcpy(newV->F0D, V.F0D, sizeof(double)*V.L);
xue@1:   }
xue@1:   return newV;
xue@1: }//InterpoalteV
xue@1: 
xue@1: //vibrato rate boundaries, in Hz
xue@1: #define MAXMOD 15.0
xue@1: #define MINMOD 3.0
xue@1: 
Chris@5: /**
xue@1:   function RateAndReg: evaluates modulation rate and regularity
xue@1: 
xue@1:   In: data[frst:fren]: modulated signal, time measured in frames
xue@1:       padrate: padding rate used for inverse FFT
xue@1:       sps: sampling frequency
xue@1:       offst: hop size
xue@1:   Out: rate, regularity: rate and regularity
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void RateAndReg(double& rate, double& regularity, double* data, int frst, int fren, int padrate, double sps, double offst)
xue@1: {
xue@1:   rate=0; regularity=0;
xue@1:   //find the section used for calculating regularity and rate
xue@11:   int frlen=fren-frst+1, frlenp=1<<(int(Log2(fren))+2), frlenpp=frlenp*padrate;
xue@1:   double *lf=new double[frlenpp*4];
xue@1:   cdouble *w=(cdouble*)&lf[frlenpp];
xue@1:   cdouble *x=(cdouble*)&lf[frlenpp*2];
xue@1:   SetTwiddleFactors(frlenp, w);
xue@1: 
xue@1:   memset(lf, 0, sizeof(double)*frlenp);
xue@1:   memcpy(lf, &data[frst], sizeof(double)*frlen);
xue@1: //    for (int i=0; i<frlen; i++) lf[i]=data[frst+i]-edata;
xue@1: 
xue@1:     double* rdata=new double[50];
xue@1:     for (int i=0; i<50; i++)
xue@1:     {
xue@1:       rdata[i]=0; for (int j=0; j<fren-i; j++) rdata[i]+=data[j]*data[j+i];
xue@1:     }
xue@1: 
xue@1: 
xue@11:   RFFTC(lf, NULL, NULL, Log2(frlenp), w, x, 0);
xue@1: //  xr[0]=0;
xue@1:   for (int i=0; i<frlenp/2; i++) {x[i].x=x[i].x*x[i].x+x[i].y*x[i].y; x[i].y=0;}
xue@1:   memset(&x[frlenp/2], 0, sizeof(cdouble)*(frlenpp-frlenp/2));
xue@1: 
xue@1:   SetTwiddleFactors(frlenpp, w);
xue@11:   CIFFTR(x, Log2(frlenpp), w, lf);
xue@1: 
xue@1:   delete[] rdata;
xue@1: 
xue@1:   double minpeak=sps/(offst*MAXMOD), maxpeak=sps/(offst*MINMOD);
xue@1: 
xue@1:   //find 2nd highest peak at interval padrate
xue@1:   int imaxlf=padrate*minpeak;
xue@1:   while (imaxlf<frlenpp && imaxlf<padrate*maxpeak && lf[imaxlf]<=lf[imaxlf-padrate]) imaxlf+=padrate;
xue@1:   double maxlf=lf[imaxlf];
xue@1:   for (int i=imaxlf/padrate; i<frlen/2; i++) if (maxlf<lf[i*padrate]) maxlf=lf[i*padrate], imaxlf=i*padrate;
xue@1:   //locate the peak at interval 1
xue@1:   int ii=imaxlf; for (int i=ii-padrate+1; i<ii+padrate; i++) if (maxlf<lf[i]) maxlf=lf[i], imaxlf=i;
xue@1: 
xue@1:   double a_=0.5*(lf[imaxlf-1]+lf[imaxlf+1])-maxlf, b_=0.5*(lf[imaxlf+1]-lf[imaxlf-1]);
xue@1:   double period=imaxlf-0.5*b_/a_;
xue@1:   period=period/padrate; //period in frames
xue@1:   rate=sps/(period*offst); //rate in hz
xue@1:   regularity=(maxlf-0.25*b_*b_/a_)/lf[0]*(fren-frst)/(fren-frst-period);
xue@1: 
xue@1:   delete[] lf;
xue@1: }//RateAndReg
xue@1: 
Chris@5: /**
xue@1:   function RegularizeV: synthesize a regular vibrato from the basic descriptors.
xue@1: 
xue@1:   In: V: a TVo object hosting vibrato descriptors
xue@1:       sps: sampling rate
xue@1:       h: hop size
xue@1:   Out: HS: a harmonic sinusoid
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void RegularizeV(THS& HS, TVo& V, double sps, double h)
xue@1: {
xue@1:   int K=V.Kp[1]; for (int p=2; p<V.P; p++) if (K>V.Kp[p]) K=V.Kp[p];
xue@1:   double D=0, *fxr=new double[2*K-1];
xue@1:   memset(fxr, 0, sizeof(double)*(2*K-1));
xue@1:   for (int p=1; p<V.P; p++)
xue@1:   {
xue@1:     D+=V.Dp[p]; for (int k=0; k<2*K-1; k++) fxr[k]+=V.fxr[p][k]*V.Dp[p];
xue@1:   }
xue@1:   //uncomment the following line if re-normalization needed
xue@1:   //double alfa=fxr[0]*fxr[0]; for (int k=1; k<2*K-1; k++) alfa+=2*fxr[k]*fxr[k]; alfa=sqrt(alfa); D=alfa;
xue@1:   for (int k=0; k<2*K-1; k++) fxr[k]/=D;
xue@1:   D/=(V.P-1);
xue@1:   V.D=D;
xue@1: 
xue@1:   atom** Partials=HS.Partials;
xue@1:   for (int l=0; l<V.L; l++)
xue@1:   {
xue@1:     double theta=V.rate*l*V.h/sps;//-0.25;
xue@1:     //f0(theta)
xue@1:     double f0=fxr[0]; for (int k=1; k<K; k++) f0=f0+2*(fxr[2*k-1]*sin(2*M_PI*k*theta)+fxr[2*k]*cos(2*M_PI*k*theta));
xue@1:     V.F0D[l]=D*f0; V.F0[l]=V.F0C[l]+V.F0D[l];
xue@1:     for (int m=0; m<V.M; m++)
xue@1:     {
xue@1:       if (Partials[m][l].t<0) Partials[m][l].t=Partials[0][l].t;
xue@1:       Partials[m][l].f=(m+1)*V.F0[l];
xue@1:       if (Partials[m][l].f>0.5 || Partials[m][l].f<0) Partials[m][l].f=-1, Partials[m][l].a=0;
xue@1:       else Partials[m][l].a=V.A0C[l]*exp(V.LogAS[m]+V.LogAF[int(Partials[m][l].f*V.afres)]);
xue@1:     }
xue@1:   }
xue@1: 
xue@1:   delete[] fxr;
xue@1: }//RegularizeV
xue@1: 
Chris@5: /**
xue@1:   function QIE: computes extremum using quadratic interpolation
xue@1: 
xue@1:   In: y[-1:1]: three points to interpolate from. It is assumed y[0] is larger (or smaller) than both
xue@1:         y[-1] and y[1].
xue@1: 
xue@1:   Returns the extremum of y.
xue@1: */
xue@1: double QIE(double* y)
xue@1: {
xue@1:   double a=0.5*(y[1]+y[-1])-y[0], b=0.5*(y[1]-y[-1]);
xue@1:   return y[0]-0.25*b*b/a;
xue@1: }//QIE
xue@1: 
Chris@5: /**
xue@1:   function QIE: computes extremum using quadratic interpolation
xue@1: 
xue@1:   In: y[-1:1]: three points to interpolate from. It is assumed y[0] is larger (or smaller) than both
xue@1:         y[-1] and y[1].
xue@1:   Out: x: the argument value that yields extremum of y(x)
xue@1: 
xue@1:   Returns the extremum of y.
xue@1: */
xue@1: double QIE(double* y, double& x)
xue@1: {
xue@1:   double a=0.5*(y[1]+y[-1])-y[0], b=0.5*(y[1]-y[-1]);
xue@1:   x=-0.5*b/a;
xue@1:   return y[0]-0.25*b*b/a;
xue@1: }//QIE
xue@1: 
Chris@5: /**
xue@1:   function S_F: original source-filter estimation used in AES124.
xue@1: 
xue@1:   In: afres: filter response resolution
xue@1:       Partials[M][Fr]: HS partials
xue@1:       A0C[Fr]: partial-independent amplitude carrier
xue@1:       lp[P]: cycle boundaries
xue@1:       F0Overall: average fundamental frequency
xue@1:   Out: LogAF[afres/2]: filter response
xue@1:        LogAS[M]: source amplitudes
xue@1: 
xue@1:   Returns 0.
xue@1: */
xue@1: double S_F(int afres, double* LogAF, double* LogAS, int Fr, int M, atom** Partials, double* A0C, double* lp, int P, double F0Overall)
xue@1: {
xue@1:     //SOURCE-FILTER estimation routines
xue@1:     //derivative of log A(f), f is sampled at 1/8192 (5hz), averaged over (-15hz, 15hz)
xue@1:     double *DAF=new double[afres], *NAF=&DAF[afres/2];
xue@1:     memset(DAF, 0, sizeof(double)*afres);
xue@1:     double avgwidth1, avgwidth2;
xue@1:     for (int fr=lp[1]; fr<lp[P-1]; fr++) for (int m=0; m<M; m++)
xue@1:     {
xue@1:       if (Partials[m][fr].f>0 && Partials[m][fr-1].f>0 && Partials[m][fr+1].f>0)
xue@1:       {
xue@1:         double f_1=Partials[m][fr-1].f, f0=Partials[m][fr].f, f1=Partials[m][fr+1].f, a_1=Partials[m][fr-1].a/A0C[fr-1], a1=Partials[m][fr+1].a/A0C[fr+1];
xue@1:         if ((f0-f_1)*(f0-f1)<0)
xue@1:         {
xue@1:           double dlogaf_df=(log(a1)-log(a_1))/(f1-f_1); //calculate derivative of log
xue@1: 
xue@1:           //this derivative contributes within a frequency interval of |f1-f_1|
xue@1:           if (f1<f_1) {double tmp=f1; f1=f_1; f_1=tmp;}
xue@1:           int startbin=ceil(f_1*afres), endbin=floor(f1*afres);
xue@1: 
xue@1:             if (startbin<0) startbin=0;
xue@1:             avgwidth1=(f0-f_1)*afres, avgwidth2=(f1-f0)*afres;
xue@1:             for (int i=startbin; i<=endbin; i++)
xue@1:             {
xue@1:               double fi=i-f0*afres, w;
xue@1:               if (fi<0) w=1-fabs(fi)/avgwidth1; else w=1-fabs(fi)/avgwidth2;
xue@1:               DAF[i]+=dlogaf_df*w, NAF[i]+=w;
xue@1:             }
xue@1:         }
xue@1:       }
xue@1:     }
xue@1:     for (int i=0; i<afres/2; i++) if (NAF[i]>0) DAF[i]=DAF[i]/NAF[i]; //else NAF[i]=0;
xue@1:     int i=0, nbands=0, bandsstart[100], bandsendinc[100]; double mininaf=3;
xue@1:     while (i<afres/2)
xue@1:     {
xue@1:       while (i<afres/2 && NAF[i]<mininaf) i++; if (i>=afres/2) break;
xue@1:       bandsstart[nbands]=i;
xue@1:       while (i<afres/2 && NAF[i]>0) i++; if (i>=afres/2) {bandsendinc[nbands++]=i-1; break;}
xue@1:       if (NAF[i]<=0) while (NAF[i]<mininaf) i--;
xue@1:       bandsendinc[nbands++]=i++;
xue@1:     }
xue@1: 
xue@11:   if (nbands==0)
xue@11:   {
xue@11:     memset(LogAF, 0, sizeof(double)*afres/2);
xue@11:   }
xue@11:   else
xue@11:   {
xue@1:     //integrate DAF over the fundamental band, with AF(F0Overall)=1
xue@1:     i=F0Overall*afres;
xue@1:     double theta=F0Overall*afres-i;
xue@1:     LogAF[i]=-0.5*theta*(DAF[i]+DAF[i]*(1-theta)+DAF[i+1]*theta)/afres; i--;
xue@1:     while (i>=bandsstart[0]){LogAF[i]=LogAF[i+1]-0.5*(DAF[i]+DAF[i+1])/afres; i--;}
xue@1:     i=F0Overall*afres+1;
xue@1:     LogAF[i]=0.5*(1-theta)*(DAF[i]+DAF[i-1]*(1-theta)+DAF[i]*theta)/afres; i++;
xue@1:     while (i<=bandsendinc[0]){LogAF[i]=LogAF[i-1]+0.5*(DAF[i]+DAF[i-1])/afres; i++;}
xue@1: 
xue@1:     int band=1;
xue@1:     while (band<nbands)
xue@1:     {
xue@1:       //integrate DAF over band-th band, first ignore the blank band
xue@1:       LogAF[bandsstart[band]]=LogAF[bandsendinc[band-1]];
xue@1:       for (int i=bandsstart[band]+1; i<=bandsendinc[band]; i++) LogAF[i]=LogAF[i-1]+0.5*(DAF[i]+DAF[i-1])/afres;
xue@1:       //then look for the shift required between bands
xue@1:       double cshift=0; int ccount=0;
xue@1:       for (int fr=lp[1]; fr<lp[P-1]; fr++)
xue@1:       {
xue@1:         double f=Partials[band-1][fr].f*afres;
xue@1:         if (bandsendinc[band-1]>bandsstart[band-1] && (f<bandsstart[band-1] || f>bandsendinc[band-1])) continue;
xue@1:         int intf=floor(f);
xue@1:         double logaflast=LogAF[intf]+(DAF[intf]+0.5*(DAF[intf+1]-DAF[intf])*(f-intf))*(f-intf)/afres;
xue@1:         f=Partials[band][fr].f*afres;
xue@1:         if (bandsendinc[band]>bandsstart[band] && (f<bandsstart[band] || f>bandsendinc[band])) continue;
xue@1:         intf=floor(f);
xue@1:         double logafthis=LogAF[intf]+(DAF[intf]+0.5*(DAF[intf+1]-DAF[intf])*(f-intf))*(f-intf)/afres;
xue@1:         cshift=cshift+log(Partials[band][fr].a*(band+1)/(Partials[band-1][fr].a*band))+logaflast-logafthis;
xue@1:         ccount++;
xue@1:       }
xue@11:       if (ccount!=0) cshift/=ccount;
xue@11:       else cshift=0;
xue@1:       //apply the shift
xue@1:       for (int i=bandsstart[band]; i<=bandsendinc[band]; i++) LogAF[i]+=cshift;
xue@1:       //make up the blank band
xue@1:       int lastend=bandsendinc[band-1], thisstart=bandsstart[band];
xue@1:       for (int i=lastend+1; i<thisstart; i++) LogAF[i]=LogAF[lastend]+cshift*(i-lastend)/(thisstart-lastend);
xue@1:       band++;
xue@1:     }
xue@1: 
xue@1:     int tmpband=bandsstart[0]; for (int i=0; i<tmpband; i++) LogAF[i]=LogAF[tmpband];
xue@1:     tmpband=bandsendinc[nbands-1]; for (int i=tmpband; i<afres/2; i++) LogAF[i]=LogAF[tmpband];
xue@11:   }
xue@11:   delete[] DAF;
xue@1: 
xue@1:     //Here LogAF contains the filter model, new get the source model
xue@1:     memset(LogAS, 0, sizeof(double)*100);
xue@1:     for (int m=0; m<M; m++)
xue@1:     {
xue@1:       int ccount=0;
xue@1:       for (int fr=lp[1]; fr<lp[P-1]; fr++)
xue@1:       {
xue@1:         double f=Partials[m][fr].f*afres;
xue@1:         if (f<=0) continue;
xue@1:         int intf=floor(f);
xue@1:         LogAS[m]=LogAS[m]+log(Partials[m][fr].a/A0C[fr])-(LogAF[intf]*(intf+1-f)+LogAF[intf+1]*(f-intf));
xue@1:         ccount++;
xue@1:       }
xue@1:       if (ccount>0) LogAS[m]/=ccount;
xue@1:       else LogAS[m]=0;
xue@1:     }
xue@1:   return 0;
xue@1: }//S_F
xue@1: 
Chris@5: /**
xue@1:   function S_F_SV: slow-variation source-filter estimation used in AES126, adapted from TSF to TVo.
xue@1: 
xue@1:   In: afres: filter response resolution
xue@1:       Partials[M][Fr]: HS partials
xue@1:       A0C[Fr]: partial-independent amplitude carrier
xue@1:       lp[P]: cycle boundaries
xue@1:       F: filter response sampling interval
xue@1:       FScaleMode: set if filter sampling uses mel scale
xue@1:       theta: balancing factor of amplitude and frequency variations, needed for SV approach
xue@1:   Out: LogAF[afres/2]: filter response
xue@1:        LogAS[M]: source amplitudes
xue@1: 
xue@1:   Returns 0.
xue@1: */
xue@1: double S_F_SV(int afres, double* LogAF, double* LogAS, int Fr, int M, atom** Partials, double* A0C, double* lp, int P, double F, int FScaleMode, double theta, double Fs)
xue@1: {
xue@1:   int lp0=lp[1], lpp=lp[P-1];
xue@1:   int L=lpp-lp0;
xue@1:   Alloc2(L, M, loga); Alloc2(L, M, f);
xue@1:   double fmax=0;
xue@1:   for (int l=0; l<L; l++)
xue@1:   {
xue@1:     int fr=lp0+l;
xue@1:     for (int m=0; m<M; m++)
xue@1:     {
xue@1:       f[l][m]=Partials[m][fr].f;
xue@1:       if (FScaleMode) f[l][m]=log(1+f[l][m]*Fs/700)/log(1+Fs/700);
xue@1:       if (f[l][m]>0) loga[l][m]=log(Partials[m][fr].a);
xue@1:       else loga[l][m]=-100;
xue@1:       if (fmax<f[l][m]) fmax=f[l][m];
xue@1:     }
xue@1:   }
xue@1:   int K=floor(fmax/F);
xue@1: 
xue@1:   Alloc2(L, K+2, h); memset(h[0], 0, sizeof(double)*L*(K+2));
xue@1: 
xue@1:   SF_SV(h, L, M, K, loga, f, F, theta, 1e-6, L*(K+2));
xue@1: 
xue@1:   for (int k=0; k<K+2; k++)
xue@1:   {
xue@1:     for (int l=1; l<L; l++) h[0][k]+=h[l][k];
xue@1:     h[0][k]/=L;
xue@1:   }
xue@1: 
xue@1:   for (int i=0; i<afres/2; i++)
xue@1:   {
xue@1:     double freq=i*1.0/afres;
xue@1:     int k=floor(freq/F);
xue@1:     double f_plus=freq/F-k;
xue@1:     if (k<K+1) LogAF[i]=h[0][k]*(1-f_plus)+h[0][k+1]*f_plus;
xue@1:     else LogAF[i]=h[0][K+1];
xue@1:   }
xue@1: 
xue@1:   memset(LogAS, 0, sizeof(double)*100);
xue@1:   for (int m=0; m<M; m++)
xue@1:   {
xue@1:     int ccount=0;
xue@1:     for (int fr=lp[1]; fr<lp[P-1]; fr++)
xue@1:     {
xue@1:       double f=Partials[m][fr].f*afres;
xue@1:       if (f<=0) continue;
xue@1:       int intf=floor(f);
xue@1:       LogAS[m]=LogAS[m]+log(Partials[m][fr].a/A0C[fr])-(LogAF[intf]*(intf+1-f)+LogAF[intf+1]*(f-intf));
xue@1:       ccount++;
xue@1:     }
xue@1:     if (ccount>0) LogAS[m]/=ccount;
xue@1:     else LogAS[m]=0;
xue@1:   }
xue@1: 
xue@1:   DeAlloc2(loga); DeAlloc2(f); DeAlloc2(h);
xue@1:   return 0;
xue@1: }//S_F_SV
xue@1: 
Chris@5: /**
xue@1:   function SynthesizeV: synthesizes a harmonic sinusoid from vibrato descriptors
xue@1: 
xue@1:   In: V: a TVo object hosting vibrato descriptors
xue@1:       sps: sampling rate
xue@1:       UseK: maximal number of Fourier series components to use to construct frequency modulator
xue@1:   Out: HS: a harmonic sinusoid
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void SynthesizeV(THS* HS, TVo* V, double sps, int UseK)
xue@1: {
xue@1:   int HSt0=HS->Partials[0][0].t, HSs0=HS->Partials[0][0].s, L=V->L, M=V->M,
xue@1:       P=V->P, *Kp=V->Kp, K=V->K, afres=V->afres;
xue@1:   double **fxr=V->fxr, *lp=V->lp, **LogASp=V->LogASp, *LogAF=V->LogAF, *Dp=V->Dp,
xue@1:       *F0C=V->F0C, *F0D=V->F0D, *F0=V->F0, *A0C=V->A0C, h=V->h;
xue@1:   HS->Resize(M, L);
xue@1:   double *gamp=new double[4*P], *gamm=&gamp[P];
xue@1: 
xue@1:   Alloc2(P, K*2, newfxr);
xue@1:   int* newKp=new int[P];
xue@1:   for (int p=1; p<P; p++)
xue@1:   {
xue@1:     memcpy(newfxr[p], fxr[p], sizeof(double)*2*K);
xue@1:     newKp[p]=Kp[p]; if (UseK>=2 && newKp[p]>UseK) newKp[p]=UseK;
xue@1:     CorrectFS(newKp[p], newfxr[p]);
xue@1:   }
xue@1: 
xue@1:   double f0, f1;
xue@1:   for (int p=1; p<P-1; p++)
xue@1:   {
xue@1:     if (p==1)
xue@1:     {
xue@1:       f0=newfxr[p][0]; for (int k=1; k<newKp[p]; k++) f0=f0+2*newfxr[p][2*k];
xue@1:     }
xue@1:     else f0=f1;
xue@1:     f1=newfxr[p+1][0]; for (int k=1; k<newKp[p+1]; k++) f1=f1+2*newfxr[p+1][2*k];
xue@1:     gamp[p]=(Dp[p+1]*f1-Dp[p]*f0)/2; gamm[p+1]=-gamp[p]; //additive correction
xue@1:   }
xue@1:   gamm[1]=-gamp[1], gamp[P-1]=-gamm[P-1]; //additive correction
xue@1:   atom** Partials=HS->Partials;
xue@1:   for (int p=0; p<P; p++)
xue@1:   {
xue@1:     double l0, lrange;
xue@1:     int lst, len, ilp;
xue@1:     if (p==0) lst=0, len=ceil(lp[0]), ilp=1;
xue@1:     else if (p==P-1) lst=ceil(lp[p-1]), len=L, ilp=p-1;
xue@1:     else lst=ceil(lp[p-1]), len=ceil(lp[p]), ilp=p;
xue@1:     l0=lp[ilp-1], lrange=lp[ilp]-l0;
xue@1:     for (int l=lst; l<len; l++)
xue@1:     {
xue@1:       Partials[0][l].s=HSs0;
xue@1:       Partials[0][l].t=HSt0+h*l;
xue@1:       double theta=(l-l0)/lrange;
xue@1:       double f0=newfxr[ilp][0]; for (int k=1; k<newKp[ilp]-1; k++) f0=f0+2*(newfxr[ilp][2*k-1]*sin(2*M_PI*k*theta)+newfxr[ilp][2*k]*cos(2*M_PI*k*theta));
xue@1:       F0D[l]=f0*Dp[ilp]+gamm[ilp]*(1-theta)+gamp[ilp]*theta;   //additive correction
xue@1: //        V.F0D[l]=f0*V.Dp[lp];   //no correction
xue@1:       F0[l]=F0C[l]+F0D[l];
xue@1:       for (int m=0; m<M; m++)
xue@1:       {
xue@1:         Partials[m][l].s=HSs0;
xue@1:         Partials[m][l].t=Partials[0][l].t;
xue@1:         Partials[m][l].f=(m+1)*F0[l];
xue@1:         if (Partials[m][l].f>0.5 || Partials[m][l].f<0) Partials[m][l].a=0;
xue@1:         else if (ilp==1) Partials[m][l].a=A0C[l]*exp(LogASp[1][m]*(1-0.5*theta)+LogASp[2][m]*0.5*theta+LogAF[int(Partials[m][l].f*afres)]);
xue@1:         else if (ilp==P-1) Partials[m][l].a=A0C[l]*exp(LogASp[ilp][m]*0.5*(1+theta)+LogASp[ilp-1][m]*0.5*(1-theta)+LogAF[int(Partials[m][l].f*afres)]);
xue@1:         else
xue@1: //					  Partials[m][l].a=A0C[l]*exp(LogASp[ilp][m]*0.5+LogASp[ilp-1][m]*0.5*(1-theta)+LogASp[ilp+1][m]*0.5*theta+LogAF[int(Partials[m][l].f*afres)]);
xue@1:           Partials[m][l].a=A0C[l]*exp(LogASp[p][m]+LogAF[int(Partials[m][l].f*afres)]);
xue@1:       }
xue@1:     }
xue@1:   }
xue@1:   delete[] gamp;
xue@1:   DeAlloc2(newfxr);
xue@1:   delete[] newKp;
xue@1: //    VibratoDemoForm->Label13->Caption=L;
xue@1: }//SynthesizeV
xue@1: 
xue@1: 
xue@1: 
xue@1: 
xue@1: 
xue@1: