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 "hsedit.h"
xue@1: #include "splines.h"
xue@1: 
Chris@5: /** \file hsedit.h */
Chris@5: 
xue@1: //---------------------------------------------------------------------------
xue@1: 
Chris@5: /**
xue@1:   function DeFM: frequency de-modulation
xue@1: 
xue@1:   In: peakfr[npfr]: segmentation into FM cycles, peakfr[0]=0, peakfr[npfr-1]=Fr-1
xue@1:       a1[Fr], f1[Fr]: sequence of amplitudes and frequencies
xue@1:       arec[Fr]: amplitude-based weights for frequency averaging
xue@1:   Out: a2[Fr], f2[Fr]: amplitude and frequency sequance after demodulation
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void DeFM(double* a2, double* f2, double* a1, double* f1, double* arec, int npfr, int* peakfr)
xue@1: {
xue@1:   double *frs=new double[npfr*12], *a=&frs[npfr], *f=&frs[npfr*2],
xue@1:         *aa=&frs[npfr*3], *ab=&frs[npfr*4], *ac=&frs[npfr*5], *ad=&frs[npfr*6],
xue@1:         *fa=&frs[npfr*7], *fb=&frs[npfr*8], *fc=&frs[npfr*9], *fd=&frs[npfr*10];
xue@1:   a[0]=a1[0], f[0]=f1[0], frs[0]=peakfr[0];
xue@1:   for (int i=1; i<npfr-1; i++)
xue@1:   {
xue@1:     a[i]=f[i]=frs[i]=0; double lrec=0;
xue@1:     for (int fr=peakfr[i-1]; fr<peakfr[i+1]; fr++)
xue@1:       a[i]+=a1[fr]*a1[fr], f[i]+=f1[fr]*arec[fr], frs[i]+=fr*arec[fr], lrec+=arec[fr];
xue@1:     a[i]=sqrt(a[i]/(peakfr[i+1]-peakfr[i-1])), f[i]/=lrec, frs[i]/=lrec;
xue@1:   }
xue@1:   a[npfr-1]=a1[peakfr[npfr-1]], f[npfr-1]=f1[peakfr[npfr-1]], frs[npfr-1]=peakfr[npfr-1];
xue@1:   CubicSpline(npfr-1, aa, ab, ac, ad, frs, a, 1, 1, a2);
xue@1:   CubicSpline(npfr-1, fa, fb, fc, fd, frs, f, 1, 1, f2);
xue@1:   delete[] frs;
xue@1: }//DeFM
xue@1: 
Chris@5: /**
xue@1:   function DFMSeg: segments HS frames into FM cycles
xue@1: 
xue@1:   In: partials[M][Fr]: HS partials
xue@1:   Out: peakfr[npfr]: segmentation, peakfr[0]=0, peakfr[npfr-1]=Fr-1.
xue@1:        arec[Fr]: total amplitudes of frames
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void DFMSeg(double* arec, int& npfr, int* peakfr, int M, int Fr, atom** partials)
xue@1: {
xue@1:   double *frec=new double[Fr];
xue@1:   memset(arec, 0, sizeof(double)*Fr); memset(frec, 0, sizeof(double)*Fr);
xue@1:   for (int m=0; m<M; m++) for (int fr=0; fr<Fr; fr++) {double la=partials[m][fr].a; la=la*la; arec[fr]+=la; frec[fr]+=partials[m][fr].f/(m+1)*la;}
xue@1:   for (int fr=0; fr<Fr; fr++) frec[fr]=frec[fr]/arec[fr];
xue@1:   peakfr[0]=0; npfr=1;
xue@1:   for (int fr=1; fr<Fr-1; fr++)
xue@1:   {
xue@1:     if ((frec[fr]<frec[fr-1] && frec[fr]<frec[fr+1]) || (frec[fr]>frec[fr-1] && frec[fr]>frec[fr+1]))
xue@1:     {
xue@1:       peakfr[npfr]=fr;
xue@1:       if (peakfr[npfr]-peakfr[npfr-1]>2) npfr++;
xue@1:     }
xue@1:   }
xue@1:   peakfr[npfr++]=Fr-1;
xue@1:   delete[] frec;
xue@1: }//DFMSeg
xue@1: 	
Chris@5: /**
xue@1:   function HSAM: harmonic sinusoid amplitude modulation
xue@1: 
xue@1:   In: SrcHS: source harmonic sinusoid
xue@1:       dep: modulation depth
xue@1:       fre: modulator frequency
xue@1:       ph: modulator phase
xue@1:   Out: HS: destination harmonic sinusoid
xue@1: 
xue@1:   No reutrn value.
xue@1: */
xue@1: void HSAM(THS* HS, THS* SrcHS, double dep, double fre, double ph)
xue@1: {
xue@1:   double omg=M_PI*2*fre;
xue@1:   for (int m=0; m<HS->M; m++)
xue@1:     for (int fr=0; fr<HS->Fr; fr++)
xue@1:       HS->Partials[m][fr].a=SrcHS->Partials[m][fr].a*(1+dep*cos(omg*SrcHS->Partials[m][fr].t+ph));
xue@1: }//HSAM
xue@1: 
Chris@5: /**
xue@1:   function HSFM: harmonic sinusoid frequency modulation
xue@1: 
xue@1:   In: SrcHS: source harmonic sinusoid
xue@1:       a: modulation extent, in semitones
xue@1:       fre: modulator frequency
xue@1:       ph: modulator phase
xue@1:   Out: HS: destination harmonic sinusoid
xue@1: 
xue@1:   No reutrn value.
xue@1: */
xue@1: void HSFM(THS* HS, THS* SrcHS, double a, double freq, double ph)
xue@1: {
xue@1:   double omg=M_PI*2*freq, pa=pow(2, a/12.0)-1;
xue@1:   for (int m=0; m<HS->M; m++)
xue@1:     for (int fr=0; fr<HS->Fr; fr++)
xue@1:       HS->Partials[m][fr].f=SrcHS->Partials[m][fr].f*(1+pa*cos(omg*SrcHS->Partials[m][fr].t+ph));
xue@1: }//HSFM
xue@1: 
Chris@5: /**
xue@1:   function HSFM_SF: harmonic sinusoid frequency modulation with source-filter model
xue@1: 
xue@1:   In: SrcHS: source harmonic sinusoid
xue@1:       a: modulation extent, in semitones
xue@1:       fre: modulator frequency
xue@1:       ph: modulator phase
xue@1:       SF: source-filter model
xue@1:   Out: HS: destination harmonic sinusoid
xue@1: 
xue@1:   No reutrn value.
xue@1: */
xue@1: void HSFM_SF(THS* HS, THS* SrcHS, double a, double freq, double ph, TSF* SF)
xue@1: {
xue@1:   double omg=M_PI*2*freq, pa=pow(2, a/12.0)-1;
xue@1:   for (int m=0; m<HS->M; m++) for (int fr=0; fr<HS->Fr; fr++)
xue@1:   {
xue@1:     double f0=SrcHS->Partials[m][fr].f;
xue@1:     double f1=f0*(1+pa*cos(omg*SrcHS->Partials[m][fr].t+ph));
xue@1:     HS->Partials[m][fr].f=f1;
xue@1:     HS->Partials[m][fr].a=SrcHS->Partials[m][fr].a*exp(SF->LogAF(f1)-SF->LogAF(f0));
xue@1:   }
xue@1: }//HSFM_SF
xue@1: 
Chris@5: /**
xue@1:   function: HSPitchShift: harmonic sinusoid pitch shifting
xue@1: 
xue@1:   In: SrcHS: source harmonic sinusoid
xue@1:       ps12: amount of pitch shift, in semitones
xue@1:   Out: HS: destination harmonic sinusoid
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void HSPitchShift(THS* HS, THS* SrcHS, double ps12)
xue@1: {
xue@1:   double pa=pow(2, ps12/12.0);
xue@1:   for (int m=0; m<HS->M; m++) for (int fr=0; fr<HS->Fr; fr++) HS->Partials[m][fr].f=SrcHS->Partials[m][fr].f*pa;
xue@1: }//HSPitchShift
xue@1: 
Chris@5: /**
xue@1:   function ReFM: frequency re-modulation
xue@1: 
xue@1:   In: partials[M][Fr]: HS partials
xue@1:       amount: relative modulation depth after remodulation
xue@1:       rate: relateive modulation rate after remodulation
xue@1:       SF: a source-filter model, optional
xue@1:   Out: partials2[M][Fr]: remodulated HS partials. Must be allocated before calling.
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void ReFM(int M, int Fr, atom** partials, atom** partials2, double amount, double rate, TSF* SF)
xue@1: {
xue@1:   double *arec=new double[Fr]; int *peakfr=new int[Fr], npfr;
xue@1:   DFMSeg(arec, npfr, peakfr, M, Fr, partials);
xue@1: 
xue@1:   double *a1=new double[Fr*8];
xue@1:   double *f1=&a1[Fr], *a2=&a1[Fr*3], *f2=&a1[Fr*4], *da=&a1[Fr*5], *df=&a1[Fr*6];
xue@1: 
xue@1:   for (int m=0; m<M; m++)
xue@1:   {
xue@1:     atom *part=partials[m], *part2=partials2[m]; bool fzero=false;
xue@1:     for (int fr=0; fr<Fr; fr++)
xue@1:     {
xue@1:       if (part[fr].f<=0){fzero=true; break;}
xue@1:       a1[fr]=part[fr].a*2;
xue@1:       f1[fr]=part[fr].f;
xue@1:     }
xue@1:     if (fzero){part2[0].f=0; break;}
xue@1:     DeFM(a2, f2, a1, f1, arec, npfr, peakfr);
xue@1:     for (int i=0; i<Fr; i++) da[i]=a1[i]-a2[i], df[i]=f1[i]-f2[i];
xue@1:     for (int fr=0; fr<Fr; fr++)
xue@1:     {
xue@1:       double frd=fr/rate; int dfrd=floor(frd); frd-=dfrd;
xue@1:       double lda=0, ldf=0;
xue@1:       if (dfrd<Fr-1) lda=da[dfrd]*(1-frd)+da[dfrd+1]*frd, ldf=df[dfrd]*(1-frd)+df[dfrd+1]*frd;
xue@1:       else if (dfrd==Fr-1) lda=da[dfrd]*(1-frd), ldf=df[dfrd]*(1-frd);
xue@1:       part2[fr].f=f2[fr]=f2[fr]+ldf*amount;
xue@1:       if (SF) part2[fr].a=part[fr].a*exp(SF->LogAF(part2[fr].f)-SF->LogAF(part[fr].f));
xue@1:       else part2[fr].a=(a2[fr]+lda*amount)*0.5;
xue@1:     }
xue@1:   }
xue@1:   delete[] a1;
xue@1:   delete[] arec; delete[] peakfr;
xue@1: }//ReFM
xue@1: 
xue@1: 
xue@1: 
xue@1: