annotate hsedit.cpp @ 13:de3961f74f30 tip

Add Linux/gcc Makefile; build fix

/*
    Harmonic sinusoidal modelling and tools

    C++ code package for harmonic sinusoidal modelling and relevant signal processing.
    Centre for Digital Music, Queen Mary, University of London.
    This file copyright 2011 Wen Xue.

    This program is free software; you can redistribute it and/or
    modify it under the terms of the GNU General Public License as
    published by the Free Software Foundation; either version 2 of the
    License, or (at your option) any later version.
*/
//---------------------------------------------------------------------------

#include "hsedit.h"
#include "splines.h"

/** \file hsedit.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