diff hsedit.cpp @ 1:6422640a802f

first upload

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