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: #ifndef hsH
xue@1: #define hsH
xue@1: 
Chris@5: /**
Chris@5:   \file hs.h - harmonic sinusoid model
xue@1: 
xue@1:   Further reading: Wen X. and M. Sandler, "Sinusoid modeling in a harmonic context," in Proc. DAFx'07, Bordeaux, 2007.
xue@1: */
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@2: #include <string.h>
xue@1: #include "arrayalloc.h"
xue@1: #include "align8.h"
xue@1: #include "fft.h"
Chris@2: #include "quickspec.h"
xue@11: 
xue@11: #ifndef __BORLANDC__
Chris@2: #include "tstream.h"
xue@11: #else
xue@11: #include <Classes.hpp>
xue@11: #endif
xue@1: 
xue@1: #define ATOM_LOCALANCHOR 1
xue@1: #define HS_CONSTF 1
xue@1: #define MAX_CAND 64
xue@1: #define STIFF_B_MAX 0.01
xue@1: 
xue@1: enum atomtype //type flags of an HS atom
xue@1: {
xue@1:   atAnchor,   //"anchor" is an atom whose validity is taken for granted, e.g. user input
xue@1:   atPeak,     //an atom whose frequency is an actual spectral peak
xue@1:   atInfered,  //an atom which has no spectral peak and whose frequency is inferred from other atoms
xue@1:   atMuted,    //an atom which is being muted
xue@1:   atBuried    //an atom which is buried in background noise
xue@1: };
xue@1: 
xue@1: struct atom   //an atom of a harmonic sinusoid
xue@1: {
xue@1:   double t;   //time
xue@1:   double s;   //scale
xue@1:   double f;   //digital frequency
xue@1:   double a;   //amplitude
xue@1:   double p;   //phase angle
xue@1:   int pin;    //partial index
xue@1:   atomtype type;  //atom type
xue@1:   int tags;   //additional info on atom type
xue@1:   float r1;
xue@1: };
xue@1: 
xue@1: struct candid //candidate atom
xue@1: {
xue@1:   int p;      //partial index
xue@1:   double f;   //frequency
xue@1:   double df;  //delta freqdel folk
xue@1:   double s;   //score in partial-sum
xue@1:   int prev;   //index of last partial f-df
xue@1:   int type;   //0: f0, 1: local maximum, 2: not a maxixum
xue@1:   double ms;
xue@1: };
xue@1: 
xue@1: struct dsparams1      //argument structure for calling ds1()
xue@1: {
xue@1:   double s;           //score
xue@1:   double lastene;     //energy of last frame
xue@1:   double currentacce; //energy of this frame, currently accumulated
xue@1:   int p;              //partial index
xue@1:   int lastP;          //number of efficient partials in last frame
xue@1:   double* lastvfp;    //amplitude estimates of partials at last frame
xue@1: };
xue@1: 
xue@1: struct NMResults      //note match output structure
xue@1: {
xue@1:   double* fp;         //frequencies, in bins
xue@1:   double* vfp;        //amplitudes
xue@1:   double* pfp;        //phase angles
xue@1:   atomtype* ptype;    //atom types
xue@1: };
xue@1: 
xue@1: struct NMSettings     //note match algorithm settings
xue@1: {
xue@1:   double c[4];        //cosine-family window specifiers
xue@1:   int M;              //ditto
xue@1:   double iH2;         //ditto
xue@1:   double hB;          //spectral truncation half width
xue@1:   int maxp;           //maximal number of partials
xue@1:   double maxB;        //stiffness coefficient upper bound
xue@1:   double epf;         //frequency estimation error tolerance for LSE estimation
xue@1:   double epf0;        //input frequency error bound for harmonic grouping
xue@1:   double delm;        //frequency error bound for harmonic grouping
xue@1:   double delp;        //pitch jump upper bound
xue@1:   double minf0;       //minimal fundamental frequency
xue@1:   double maxf0;       //maximal fundamental frequency
xue@1:   int pin0;           //input partial index
xue@1:   bool forcepin0;     //force the peak nearest to input frequency as an atom
xue@1:   bool pin0asanchor;  //mark atom found near the input frequency as anchor
xue@1:   int pcount;         //number of "pinned" (anchor) partials
xue@1:   int* pin;           //partial indices of pinned partials
xue@1:   double* pinf;       //frequencies of pinned partials
xue@1:   int* pinfr;         //frame indices of pinned partials, used in multi-frame constant-pitch note match
xue@1: };
xue@1: 
xue@1: 
Chris@5: /**
xue@1:   stiffcandid is the harmonic atom class used internally by harmonic grouping and tracking routines. Literally
xue@1:   it means "candidate harmonic atoms on a stiff string model". The stiff string is the main harmonic model used
xue@1:   by me for describing frequency relations between partials of a harmonic sound.
xue@1: 
xue@1:   stiffcandid is superceded by TTempAtom class.
xue@1: */
xue@1: 
xue@1: class stiffcandid
xue@1: {
xue@1: public:
xue@1:   int P;  //number of partial estimates located
xue@1:   int* p;   //partial indices of located partials, sizeof(int)*P
xue@1:   double* f;  //frequencies of located partials, sizeof(double)*P
xue@1:   double s; //score in partial-sum
xue@1:   //{N; F, G}: the feasible polygonal area of (F, G) given the P partials, where F=F0*F0, G=F0*B
xue@1:   int N;
xue@1:   double* F;  //sizeof(double)*N
xue@1:   double* G;  //sizeof(double)*N
xue@1: 
xue@1:   stiffcandid(int Wid); //create empty with minimal info
xue@1:   stiffcandid(int Wid, int ap, double af, double errf); //create empty with pitch range
xue@1:   stiffcandid(int Wid, int ap, double af, double errf, double ds); //create with 1 atom
xue@1:   stiffcandid(stiffcandid* prev, int ap, double af, double errf, double ds); //create by updating with new atom
xue@1:   ~stiffcandid();
xue@1: };
xue@1: 
xue@1: 
Chris@5: /**
xue@1:   THS is the data structure hosting a harmonic sinusoid representation. Its key members include the number
xue@1:   of partials, number of frames, and all its atoms (sinusoid parameters of all partials at measurement points).
xue@1: */
xue@1: class THS
xue@1: {
xue@1: public:
xue@1:   int Channel;      //channel id: THS describes a harmonic sinusoid in single channel
xue@1:   int M;            //number of partials
xue@1:   int Fr;           //number of frames
xue@1:   atom** Partials;  //atoms, [M][Fr]
xue@1: 
xue@1:   int* BufM[128];   //a buffer for algorithmic use
xue@1: 
xue@1:   int isconstf;     //constant frequencies flag
xue@1: 
xue@1:   double** startamp;//onset amplifiers, optional
xue@1:   int st_start;     //position of the first onset amplifying point
xue@1:   int st_offst;     //interval of onset amplifying points
xue@1:   int st_count;     //number of onset amplifying points
xue@1: 
xue@1:   //constructors and destructor
xue@1:   THS();
xue@1:   THS(int aM, int aFr);
xue@1:   THS(THS* HS);
xue@1:   THS(THS* HS, double start, double end);
xue@1:   ~THS();
xue@1: 
xue@1:   int StartPos();   //start position
xue@1:   int EndPos();     //end position
xue@1:   int EndPosEx();   //extended end position
xue@1:   int StdOffst();   //hop size
xue@1: 
xue@1:   void ClearBufM(); //free buffers registered in BufM[]
xue@1:   void Resize(int aM, int aFr, bool copydata=false, bool forcealloc=false); //change size (M or Fr)
xue@1: 
xue@1:   //I/O routines
xue@1:   int WriteHdrToStream(TStream* Stream);
xue@1:   int ReadHdrFromStream(TStream* Stream);
xue@1:   int WriteAtomToStream(TStream* Stream, atom* atm);
xue@1:   int ReadAtomFromStream(TStream* Stream, atom* atm);
xue@1:   int WriteToStream(TStream* Stream);
xue@1:   int ReadFromStream(TStream* Stream);
xue@1: };
xue@1: 
xue@1: 
Chris@5: /**
xue@1:   TPolygon is a polygon class. This class itself does not enforce convexness. However, when used for solving
xue@1:   the stiff string model, all polygons are convex, as they are the results of current a first convex polygon
xue@1:   by straight lines.
xue@1: 
xue@1:   For the convenience of computation, the sequence of vertices are arranged so that they come in clockwise
xue@1:   order starting from the leftmost (smallest x coordinate) vertex; in case two vertices are both leftmost
xue@1:   the upper (larger y coordinate) one is placed at the start and the lower one is placed at the end.
xue@1: */
xue@1: 
xue@1: class TPolygon
xue@1: {
xue@1: public:
xue@1:   int N;      //number of vertices (equals number of edges)
xue@1:   double* X;  //x-coordinates of vertices
xue@1:   double* Y;  //y-coordinates of vertices
xue@1:   TPolygon(int cap);
xue@1:   TPolygon(int cap, TPolygon* R);
xue@1:   ~TPolygon();
xue@1: }; 
xue@1: 
xue@1: 
Chris@5: /**
xue@1:   TTempAtom is an atom class within the stiff-stirng harmonic atom context used internally by harmonic
xue@1:   grouping and tracking routines. Literally it means "a temporary atom", and truly it hosts most info
xue@1:   one expects of a potential atom. TTempAtom replaces stiffcandid class in harmonic sinusoid group and
xue@1:   tracking algorithms, due to the advanges that it carries atom information (frequency, amplitude, etc.)
xue@1:   as well as harmonic atom information (link to other atom, the F-G polygon), and is therefore both an
xue@1:   individual atom and a harmonic atom containing all atoms tracked down through the linked list Prev.
xue@1: */
xue@1: 
xue@1: class TTempAtom
xue@1: {
xue@1: public:
xue@1:   int pind;     //partial index
xue@1:   double f;     //frequency
xue@1:   double a;     //amplitude
xue@1:   double s;     //scale
xue@1:   double rsr;   //residue-sinusoid ratio
xue@1:   TTempAtom* Prev;  //previous atom
xue@1:   TPolygon *R;  //F-G polygon for harmonic group
xue@1:   union {double acce; int tag[2];};
xue@1: 
xue@1:   TTempAtom(double af, double ef, double maxB); //create empty with frequency range
xue@1:   TTempAtom(int apin, double af, double ef, double maxB); //create empty with frequency range
xue@1:   TTempAtom(TPolygon* AR, double delf1, double delf2, double minf); //create empty with extended R
xue@1:   TTempAtom(int apind, double af, double ef, double aa, double as, double maxB);  //create with one partial
xue@1:   TTempAtom(TTempAtom* APrev, bool DupR); //create duplicate
xue@1:   TTempAtom(TTempAtom* APrev, int apind, double af, double ef, double aa, double as, bool updateR=true); //duplicate and add one partial
xue@1:   ~TTempAtom();
xue@1: };
xue@1: 
xue@1: //--internal function--------------------------------------------------------
xue@1: double ds0(double, void*); //a score function, internal use only
xue@1: 
xue@1: //--general polygon routines-------------------------------------------------
xue@1: void areaandcentroid(double& A, double& cx, double& cy, int N, double* x, double* y); //compute area and centroid
xue@1: void cutcvpoly(int& N, double* x, double* y, double A, double B, double C, bool protect=false); //sever polygon by line
xue@1: double maximalminimum(double& x, double& y, int N, double* sx, double* sy); //maximum inscribed circle
xue@1: 
xue@1: //--F-G polygon routines-----------------------------------------------------
xue@1: void CutR(TPolygon* R, int apind, double af, double ef, bool protect=false);
xue@1: void ExBStiff(double& Bmin, double& Bmax, int N, double* F, double* G);
xue@1: void ExFmStiff(double& Fmin, double& Fmax, int m, int N, double* F, double* G);
xue@1: void ExtendR(TPolygon* R, double delf1, double delf2, double minf);
xue@1: void InitializeR(TPolygon* R, double af, double ef, double maxB);
xue@1: void InitializeR(TPolygon* R, int apind, double af, double ef, double maxB);
xue@1: 
xue@1: //--internal structure conversion routines-----------------------------------
xue@11: void AtomsToPartials(int k, atom* part, int& M, int& Fr, atom**& partials, int offst);
xue@1: int NMResultToAtoms(int M, atom* HP, int t, int wid, NMResults results);
xue@1: int NMResultToPartials(int M, int fr, atom** Partials, int t, int wid, NMResults results);
xue@1: 
xue@1: //--batch sinusoid estimation routines---------------------------------------
xue@1: double PeakShapeC(double f, int Fr, int N, cdouble** x, int B, int M, double* c, double iH2);
xue@1: double PeakShapeC(double f, int Fr, int N, cfloat** x, int B, int M, double* c, double iH2);
xue@1: int QuickPeaks(double* f, double* a, int N, cdouble* x, int M, double* c, double iH2, double mina, int binst=-1, int binen=-1, int B=5, double* rsr=0);
xue@1: int QuickPeaks(double* f, double* a, int Fr, int N, cdouble** x, int fr0, int r0, int M, double* c, double iH2, double mina, int binst=-1, int binen=-1, int B=5, double* rsr=0);
xue@1: 
xue@1: //--harmonic atom detection (harmonic grouping) routines---------------------
xue@1: double NoteMatchStiff3(TPolygon* R, double &f0, double& B, int pc, double* fps, double* vps, int Fr, cdouble** x, int N, int offst, NMSettings* settings, NMResults* results, int lastp, double* lastvfp, double (*computes)(double a, void* params)=ds0, int forceinputlocalfr=-1); //basic grouping without rsr
xue@1: double NoteMatchStiff3(TPolygon* R, double &f0, double& B, int pc, double* fps, double* vps, double* rsr, int Fr, cdouble** x, int N, int offst, NMSettings* settings, NMResults* results, int lastp, double* lastvfp, double (*computes)(double a, void* params)=ds0, int forceinputlocalfr=-1); //basic grouping with rsr
xue@1: double NoteMatchStiff3(TPolygon* R, double &f0, double& B, int Fr, cdouble** x, int N, int offst, NMSettings* settings, NMResults* results, int lastp, double* lastvfp, double (*deltas)(double a, void* params)=ds0, bool forceinputlocalfr=false, int startfr=-1, int validfrrange=0); //basic grouping with rsr - wrapper
xue@1: double NoteMatchStiff3(TPolygon* R, int peak0, int pin0, cdouble* x, int pc, double* fps, double* vps, int N, NMSettings* settings, double* vfp, int** pitchind, int newpc); //grouping with given pitch, single-frame
xue@1: 
xue@1: //--harmonic sinusoid tracking routines--------------------------------------
xue@1: int FindNote(int _t, double _f, int& M, int& Fr, atom**& partials, int frst, int fren, int wid, int offst, TQuickSpectrogram* Spec, NMSettings settings); //harmonic sinusoid tracking (forward and backward tracking)
xue@1: int FindNoteConst(int _t, double _f, int& M, int& Fr, atom**& partials, int frst, int fren, int wid, int offst, TQuickSpectrogram* Spec, NMSettings settings, double brake); //constant-pitch harmonic sinusoid tracking
xue@1: int FindNoteF(atom* part, double& starts, TPolygon* R, int startp, double* startvfp, int frst, int fren, int wid, int offst, TQuickSpectrogram* Spec, NMSettings settings, double brake); //forward harmonic sinusoid tracking
xue@1: int FindNoteFB(int frst, TPolygon* Rst, double* vfpst, int fren, TPolygon* Ren, double* vfpen, int M, atom** partials, int wid, int offst, TQuickSpectrogram* Spec, NMSettings settings); //forward-backward harmonic sinusoid tracking
xue@1: void NoteMatchStiff3FB(int& pitchcount, TPolygon**& R, double**& vfp, double*& sc, int* newpitches, int* prev, int pc, double* fps, double* vps, cdouble* x, int wid, int maxpitch, NMSettings* settings); //single DP step of FindNoteFB()
xue@1: 
xue@1: //--harmonic sinusoid synthesis routines-------------------------------------
xue@1: double* SynthesisHS(int pm, int pfr, atom** partials, int& dst, int& den, bool* terminatetag=0);
xue@11: double* SynthesisHS(THS* HS, int& dst, int& den, bool* terminatetag=0);
xue@1: double* SynthesisHS2(int M, int Fr, atom** partials, int& dst, int& den, bool* terminatetag=0);
xue@1: double* SynthesisHSp(int pm, int pfr, atom** partials, int& dst, int& den, double** startamp=0, int st_start=0, int st_offst=0, int st_count=0);
xue@1: double* SynthesisHSp(THS* HS, int& dst, int& den);
xue@1: 
xue@1: //--other functions----------------------------------------------------------
xue@1: void ReEstHS1(THS* HS, __int16* Data16); //reestimation of HS
xue@1: 
xue@1: #endif