xue@11: 
xue@9: /*
xue@9:   Sample program: analysis of harmonic sinusoid in low-noise monophonic context.
xue@9: 
xue@9:   Syntax:
xue@9:   (executable path) filename <setting=value> <setting=value> ... <setting=value>
xue@9: 
xue@9:   filename: path input wave form audio file, 16-bit PCM
xue@9:   setting: optional algorithmic parameters; search "stricmp" in the source code for a complete list
xue@9:   value: non-default values given to the algorithmic parameters
xue@9: 
xue@9:   Output:
xue@9:   a *.evt file containing RIFF chunks, each starting with id "EVT\0" and hosting a harmonic sinusoid.
xue@9: 
xue@9:   The "EVT\0" chunk body contains a header ("HDR\0") chunk followed by atom ("ATM\0") chunks.
xue@9: 
xue@9:   The "HDR\0" chunk body contains 4 __int32 values: channel index, number of partials, number of
xue@9:   measurement points (frames), and a tag.
xue@9: 
xue@9:   The "ATM\0" chunk body contains an atom structure (see hs.h).
xue@9: 
xue@9: */
xue@9: 
xue@9: #include <QtCore/QCoreApplication>
xue@9: #include "arrayalloc.h"
xue@9: #include "hs.h"
xue@9: #include "matrix.h"
xue@9: #include "quickspec.h"
xue@9: #include "procedures.h"
xue@9: #include "vibrato.h"
xue@9: #include <string.h>
xue@9: #include <stdio.h>
xue@9: #include <conio.h>
xue@9: 
xue@9: //structure hosting wave header data
xue@9: struct wavehdr
xue@9: {
xue@9:   __int16 fmttag;
xue@9:   __int16 channels;
xue@9:   __int32 samplespersec;
xue@9:   __int32 bytespersec;
xue@9:   __int16 blockalign;
xue@9:   __int16 bitspersample;
xue@9: };
xue@9: 
xue@9: //tells if a wave header contains consistent data
xue@9: bool isvalidwavehdr(wavehdr* hdr)
xue@9: {
xue@9:   if (hdr->fmttag!=1) return false;
xue@9:   if (hdr->bitspersample!=8 && hdr->bitspersample!=16) return false;
xue@9:   if (hdr->channels*hdr->bitspersample!=hdr->blockalign*8) return false;
xue@9:   if (hdr->samplespersec*hdr->blockalign!=hdr->bytespersec) return false;
xue@9:   return true;
xue@9: }
xue@9: 
xue@9: //structure hosting a waveform audio
xue@9: struct waveaudio
xue@9: {
xue@9:   union
xue@9:   {
xue@9:     wavehdr hdr;
xue@9:     struct
xue@9:     {
xue@9:       __int16 fmttag;
xue@9:       __int16 channels;
xue@9:       __int32 samplespersec;
xue@9:       __int32 bytespersec;
xue@9:       __int16 blockalign;
xue@9:       __int16 bitspersample;
xue@9:     };
xue@9:   };
xue@9:   int length;
xue@9:   union
xue@9:   {
xue@9:     unsigned char* data;
xue@9:     unsigned char* data8;
xue@9:     __int16* data16;
xue@9:   };
xue@9:   waveaudio(){memset(this, 0, sizeof(waveaudio));}
xue@9:   ~waveaudio(){delete[] data;}
xue@9: };
xue@9: 
xue@9: void freewaveaudio(waveaudio* wa)
xue@9: {
xue@9:   delete[] wa->data;
xue@9:   memset(wa, 0, sizeof(waveaudio));
xue@9: }
xue@9: 
xue@9: //error messages for reading waveform audio file
xue@9: enum err_wavefile
xue@9: {
xue@9:   err_success,
xue@9:   err_RIFF_hdr,
xue@9:   err_wave_hdr,
xue@9:   err_data_hdr,
xue@9:   err_data_chunk,
xue@9:   err_io
xue@9: };
xue@9: 
xue@9: //reads wave header from file stream
xue@9: int loadwavehdr(wavehdr* hdr, FILE* fp)
xue@9: {
xue@9:   int rc, hdrlen;
xue@9:   char s[5];
xue@9: 
xue@9:   rc=fread(s, 1, 4, fp); s[4]=0;
xue@9:   if (rc<4) return err_wave_hdr;
xue@9:   if (strcmp(s, "WAVE")) return err_wave_hdr;
xue@9: 
xue@9:   rc=fread(s, 1, 4, fp); s[4]=0;
xue@9:   if (rc<4) return err_wave_hdr;
xue@9:   if (strcmp(s, "fmt ")) return err_wave_hdr;
xue@9: 
xue@9:   rc=fread(&hdrlen, 1, 4, fp);
xue@9:   if (rc<4) return err_wave_hdr;
xue@9:   if (hdrlen<16) return err_wave_hdr;
xue@9: 
xue@9:   rc=fread(hdr, 1, sizeof(wavehdr), fp);
xue@9:   if (rc<sizeof(wavehdr)) return err_wave_hdr;
xue@9:   if (!isvalidwavehdr(hdr)) return err_wave_hdr;
xue@9: 
xue@9:   if (hdrlen>16)
xue@9:   {
xue@9:     rc=fseek(fp, hdrlen-16, SEEK_CUR);
xue@9:     if (rc!=0) return err_wave_hdr;
xue@9:   }
xue@9: 
xue@9:   return err_success;
xue@9: }
xue@9: 
xue@9: //reads waveform audio, including header, from a file stream
xue@9: int loadwaveaudio(waveaudio* wa, FILE* fp)
xue@9: {
xue@9:   int rc, mainlen;
xue@9:   char s[5];
xue@9:   wavehdr hdr;
xue@9: 
xue@9:   rc=fread(s, 1, 4, fp); s[4]=0;
xue@9:   if (rc<4) return err_RIFF_hdr;
xue@9:   if (strcmp(s, "RIFF")) return err_RIFF_hdr;
xue@9: 
xue@9:   rc=fread(&mainlen, 1, 4, fp);
xue@9:   if (rc<4) return err_RIFF_hdr;
xue@9: 
xue@9:   rc=loadwavehdr(&hdr, fp);
xue@9:   if (rc!=err_success) return rc;
xue@9: 
xue@9:   rc=fread(s, 1, 4, fp); s[4]=0;
xue@9:   if (rc<4) return err_data_hdr;
xue@9:   if (strcmp(s, "data")) return err_data_hdr;
xue@9: 
xue@9:   rc=fread(&mainlen, 1, 4, fp);
xue@9:   if (rc<4) return err_data_hdr;
xue@9: 
xue@9:   if (mainlen<=0) return err_data_hdr;
xue@9: 
xue@9:   unsigned char* data=new unsigned char[mainlen];
xue@9:   rc=fread(data, 1, mainlen, fp);
xue@9:   if (rc<mainlen)
xue@9:   {
xue@9:     delete[] data;
xue@9:     return err_data_chunk;
xue@9:   }
xue@9: 
xue@9:   wa->hdr=hdr;
xue@9:   wa->length=rc/hdr.blockalign;
xue@9:   delete[] wa->data;
xue@9:   wa->data=data;
xue@9: 
xue@9:   return err_success;
xue@9: }
xue@9: 
xue@9: //reads waveform audio, including header, from a file
xue@9: int loadwavefile(waveaudio* wa, char* filename)
xue@9: {
xue@9:   FILE* fp;
xue@9:   if ((fp=fopen(filename, "rb"))==NULL) return err_io;
xue@9:   int result=loadwaveaudio(wa, fp);
xue@9:   fclose(fp);
xue@9:   return result;
xue@9: }
xue@9: 
xue@9: //returns new file path with the specific extension replacing the original one
xue@9: char* ChangeFileExt(char* FileName, char* Ext)
xue@9: {
xue@9:   char* oldext=strrchr(FileName, '.');
xue@9:   int namelen=strlen(FileName)-strlen(oldext);
xue@9:   char* dest=new char[namelen+strlen(Ext)+1];
xue@9:   memcpy(dest, FileName, namelen); dest[namelen]=0;
xue@9:   strcat(dest, Ext);
xue@9:   return dest;
xue@9: }
xue@9: 
xue@9: //ACPower for __int16 data input
xue@9: double ACPower(__int16* data, int Count)
xue@9: {
xue@9:   if (Count<=0) return 0;
xue@9:   double power=0, avg=0, tmp;
xue@9:   for (int i=0; i<Count; i++)
xue@9:   {
xue@9:     tmp=*(data++);
xue@9:     power+=tmp*tmp;
xue@9:     avg+=tmp;
xue@9:   }
xue@9:   power=(power-avg*avg/Count)/Count;
xue@9:   return power;
xue@9: }//ACPower
xue@9: 
xue@9: //double QIE(double* y, double& x){double a=0.5*(y[1]+y[-1])-y[0], b=0.5*(y[1]-y[-1]); x=-0.5*b/a; return y[0]-0.25*b*b/a;}
xue@9: 
xue@9: //correlation coefficient
xue@9: double InnerC(int N, __int16* x, __int16* y)
xue@9: {
xue@9:   double inner=0, inn1=0, inn2=0;
xue@9:   for (int n=0; n<N; n++) inner+=1.0*x[n]*y[n], inn1+=1.0*x[n]*x[n], inn2+=1.0*y[n]*y[n];
xue@9:   double inn=sqrt(inn1*inn2);
xue@9:   return (inn<=0)?0:inner/inn;
xue@9: }
xue@9: 
xue@9: //monophonic pitch estimation by autocorrelation
xue@9: double pitchautocor(__int16* data, int wid, double minf0bin, double maxf0bin, double& hsr)
xue@9: {
xue@9:   int hwid=wid/2, minT=wid/maxf0bin, maxT=wid/minf0bin;
xue@9:   double ene=InnerC(wid, data, data);
xue@9:   if (ene<=0){hsr=0; return 0;}
xue@9:   double* autocor=new double[hwid];
xue@9:   for (int i=1; i<maxT+2; i++)
xue@9:   {
xue@9:     autocor[i]=InnerC(wid-i, data, &data[i]);
xue@9:   }
xue@9: 
xue@9:   int prd=0, m=minT;
xue@9:   while(m<=maxT)
xue@9:   {
xue@9:     while (m<=maxT && (autocor[m]<0 || autocor[m]<autocor[m-1] || autocor[m]<autocor[m+1])) m++;
xue@9:     if (m<=maxT)
xue@9:     {
xue@9:       int mm=ceil(m*0.6667), mp=floor(m*1.3333), cont=0;
xue@10:       //for (int im=mm; im<=mp; im++) if (m!=im && autocor[im]>=autocor[m]) {cont=1; break;}
xue@9:       if (cont==0)
xue@9:       {
xue@9:         if (prd==0 || autocor[m]>autocor[prd]*1.05) prd=m;
xue@9:       }
xue@9:       m++;
xue@9:     }
xue@9:   }
xue@9: 
xue@9:   double pitchbin=0;
xue@9:   if (prd>=minT && prd<=maxT)
xue@9:   {
xue@9:     double mshift; hsr=QIE(&autocor[prd], mshift); pitchbin=wid/(prd+mshift);
xue@9:   }
xue@9:   else{hsr=0;}
xue@9:   delete[] autocor;
xue@9:   return pitchbin;
xue@9: }
xue@9: 
xue@9: //main function
xue@9: int main(int argc, char* argv[])
xue@9: {
xue@9:   //read audio file
xue@9:   if (argc<2){printf("Please specify input wave file."); getch(); return 0;}
xue@9:   printf("Loading wave file %s... ", argv[1]);
xue@9:   waveaudio* wa=new waveaudio;
xue@9:   int waverr=loadwavefile(wa, argv[1]); printf("[%d]\n", waverr);
xue@9:   if (waverr!=err_success){printf("Aborted: error loading wave file."); getch(); return 0;}
xue@9:   if (wa->bitspersample!=16){printf("Aborted: this program accepts 16-bit pcm only."); getch(); return 0;}
xue@9: 
xue@9:   int length=wa->length, sps=wa->samplespersec;
xue@9:   __int16* data16=wa->data16;
xue@9:   if (wa->channels>1)
xue@9:   {
xue@9:     printf("Extracting channel 0... ");
xue@9:     for (int i=1; i<length; i++) data16[i]=data16[i*wa->channels];
xue@9:     printf("[0]\n");
xue@9:   }
xue@9:   wa->data=0;
xue@9:   delete wa;
xue@9: 
xue@9:   //default settings
xue@9:   int maxhscount=100;
xue@9:   float wid_s=0.02, offstwidratio=0.5;
xue@9:   float dur_s=0.5;
xue@9:   float intv_s=0.05; //duration, in seconds, to scan for peaks
xue@9:   float minf0=55, maxf0=3520;
xue@9:   WindowType wintype=wtHann;
xue@9: 
xue@9:   NMSettings settings; memset(&settings, 0, sizeof(NMSettings));
xue@9:   settings.hB=3;          //spectral truncation half width
xue@9:   settings.maxp=50;         //maximal number of partials
xue@9:   settings.maxB=0.001;      //stiffness coefficient upper bound
xue@9:   settings.epf=1e-4;        //frequency estimation error tolerance for LSE estimation
xue@9:   settings.epf0=1;          //input frequency error bound for harmonic grouping
xue@9:   settings.delm=1.1;        //frequency error bound for harmonic grouping
xue@9:   settings.delp=1.1;        //pitch jump upper bound
xue@9: 
xue@9:   double mina=0.5;
xue@9: 
xue@9:   //user settings through commandline arguments
xue@9:   for (int i=2; i<argc; i++)
xue@9:   {
xue@9:     char* s1=strchr(argv[i], '=');
xue@9:     if (s1)
xue@9:     {
xue@9:       s1[0]=0; s1++;
xue@9:       if (!stricmp(argv[i], "hB")){settings.hB=atof(s1); continue;}
xue@9:       if (!stricmp(argv[i], "maxp")){settings.maxp=atoi(s1); continue;}
xue@9:       if (!stricmp(argv[i], "maxB")){settings.maxB=atof(s1); continue;}
xue@9:       if (!stricmp(argv[i], "epf")){settings.epf=atof(s1); continue;}
xue@9:       if (!stricmp(argv[i], "epf0")){settings.epf0=atof(s1); continue;}
xue@9:       if (!stricmp(argv[i], "delm")){settings.delm=atof(s1); continue;}
xue@9:       if (!stricmp(argv[i], "delp")){settings.delp=atof(s1); continue;}
xue@9:       if (!stricmp(argv[i], "minf0")){minf0=atof(s1); continue;}
xue@9:       if (!stricmp(argv[i], "maxf0")){maxf0=atof(s1); continue;}
xue@9:       if (!stricmp(argv[i], "wid_s")){wid_s=atof(s1); continue;}
xue@9:       if (!stricmp(argv[i], "offstwidratio")){offstwidratio=atof(s1); continue;}
xue@9:       if (!stricmp(argv[i], "dur_s")){dur_s=atof(s1); continue;}
xue@9:       if (!stricmp(argv[i], "intv_s")){intv_s=atof(s1); continue;}
xue@9:       if (!stricmp(argv[i], "maxhscount")){maxhscount=atoi(s1); continue;}
xue@9:       if (!stricmp(argv[i], "wintype"))
xue@9:       {
xue@9:         if (!stricmp(s1, "hann")){wintype=wtHann; continue;}
xue@9:         if (!stricmp(s1, "hamming")){wintype=wtHamming; continue;}
xue@9:         if (!stricmp(s1, "blackman")){wintype=wtBlackman; continue;}
xue@9:         if (!stricmp(s1, "rectangle")){wintype=wtRectangle; continue;}
xue@9:         if (!stricmp(s1, "hannsqr")){wintype=wtHannSqr; continue;}
xue@9:       }
xue@9:     }
xue@9:   }
xue@9: 
xue@9:   //
xue@9:   MList* List=new MList; List->Add(data16, 1);
xue@9: 
xue@9:   //analysis window size and hop size
xue@9:   int wid=1<<((int)floor(log2(wid_s*sps)+1)), offst=wid*offstwidratio;
xue@9: 
xue@9:   //set up the spectrogram
xue@9:   TQuickSpectrogram* Sp=new TQuickSpectrogram(0, 0, true, false, 0);
xue@9:   Sp->Data=data16; Sp->DataLength=length; Sp->BytesPerSample=2;
xue@9:   Sp->Wid=wid; Sp->Offst=offst; Sp->WinType=wintype;
xue@9:   Sp->Spec(-1); //this sets Sp->Capacity to the number of accessible frames and does nothing else
xue@9: 
xue@9:   //get number of frames
xue@9:   int Fr=Sp->Capacity;
xue@9: 
xue@9:   //the program searches every intv_fr frames over a duration of dur_fr frames for a spectral peak to start tracking.
xue@9:   int dur_fr=sps*dur_s/offst, intv_fr=sps*intv_s/offst;
xue@9:   if (dur_fr<10) dur_fr=10;
xue@9:   if (intv_fr<1) intv_fr=1;
xue@9: 
xue@9:   //minimal and maximal fundamental frequency, in bins
xue@9:   settings.minf0=minf0/sps*wid;
xue@9:   settings.maxf0=maxf0/sps*wid;
xue@9:   if (settings.minf0<settings.hB) settings.minf0=settings.hB;
xue@9: 
xue@9:   windowspec(wintype, wid, &settings.M, settings.c, &settings.iH2);
xue@9:   double *fps=new double[wid], *vps=&fps[wid/2]; List->Add(fps, 1);
xue@9:   cdouble *x=new cdouble[wid/2+1]; List->Add(x, 1);
xue@9: 
xue@9:   //output file
xue@9:   char* filename=ChangeFileExt(argv[1], ".evt"); List->Add(filename, 1);
xue@9:   TFileStream* File=new TFileStream(filename, fmWrite);
xue@9: 
xue@9:   double* frames=new double[Fr]; memset(frames, 0, sizeof(double)*Fr); List->Add(frames, 1);
xue@9:   double* framesa=new double[Fr]; memset(framesa, 0, sizeof(double)*Fr); List->Add(framesa, 1);
xue@9:   double* rsr=new double[Fr]; List->Add(rsr, 1); memset(rsr, 0, sizeof(double)*Fr);
xue@9:   double* f0s=new double[Fr]; List->Add(f0s, 1); memset(f0s, 0, sizeof(double)*Fr);
xue@9:   for (int fr=0; fr<Fr-1; fr++)
xue@9:   {
xue@9:     __int16* ldata16=&data16[offst*fr];
xue@9:     frames[fr]=ACPower(ldata16, wid);
xue@9:     framesa[fr]=1;
xue@9:  //   f0s[fr]=pitchautocor(ldata16, wid, settings.minf0, settings.maxf0, rsr[fr]); rsr[fr]=1-rsr[fr];
xue@9:   }
xue@9:   frames[Fr-1]=ACPower(&data16[(Fr-1)*offst], length-(Fr-1)*offst);
xue@9:   framesa[Fr-1]=1;
xue@9: 
xue@9:   int start=0, writecount=0;
xue@9: 
xue@9:   while (writecount<maxhscount)
xue@9:   {
xue@9:     int fr_m=start+intv_fr;
xue@9: 
xue@9:     //search for the highest spectral peak to start tracking
xue@9:     double fp_m=0, vp_m=0;
xue@9:     for (int fr=start+intv_fr; fr<start+dur_fr && fr<Fr; fr+=intv_fr)
xue@9:     {
xue@9:       cmplx<QSPEC_FORMAT>* speci=Sp->Spec(fr);
xue@9:       for (int k=0; k<=wid/2; k++) x[k]=speci[k];
xue@9:       int pc=QuickPeaks(fps, vps, wid, x, settings.M, settings.c, settings.iH2, mina, 0, wid/4);
xue@9: 
xue@9:       for(int p=0; p<pc; p++)
xue@9:         if (vps[p]>vp_m)
xue@9:         {
xue@9:           if (f0s[fr]==0 && rsr[fr]==0) f0s[fr]=pitchautocor(&data16[offst*fr], wid, settings.minf0, settings.maxf0, rsr[fr]); rsr[fr]=1-rsr[fr];
xue@9:           if (f0s[fr]>0)
xue@9:           {
xue@9:             double dpind=fps[p]/f0s[fr], pind=floor(dpind+0.5);
xue@9:             if (fabs(dpind-pind)<0.1) vp_m=vps[p], fr_m=fr, fp_m=fps[p];
xue@9:           }
xue@9:         }
xue@9:     }
xue@9: 
xue@9:     if (fp_m<=0){start=start+dur_fr; if (start>Fr) break; continue;}
xue@9: 
xue@9:     f0s[fr_m]=pitchautocor(&data16[offst*fr_m], wid, settings.minf0, settings.maxf0, rsr[fr_m]); rsr[fr_m]=1-rsr[fr_m];
xue@9: 
xue@9:     int _t=fr_m*offst+wid/2;
xue@9:     double _f=fp_m/wid;
xue@9:     int hsM, hsFr, frst=start, fren=Fr;
xue@9:     atom** hsPartials=0;
xue@9: 
xue@9:     settings.pin0=floor(fp_m/f0s[fr_m]+0.5);
xue@9:     int tag=FindNote(_t, _f, hsM, hsFr, hsPartials, frst, fren, wid, offst, Sp, settings);
xue@9: 
xue@9:     int fr1=(hsPartials[0][0].t-wid/2)/offst;
xue@9:     double consi=1;
xue@9:     for (int hsfr=0; hsfr<hsFr; hsfr++)
xue@9:     {
xue@9:       double ene=0;
xue@9:       for (int m=0; m<hsM; m++)
xue@9:       {
xue@9:         double a=(hsPartials[m][hsfr].f>0)?hsPartials[m][hsfr].a:0;
xue@9: 
xue@9:         if (hsfr==0) ene+=a*a;
xue@9:         else
xue@9:         {
xue@9:           double b=(hsPartials[m][hsfr-1].f>0)?hsPartials[m][hsfr-1].a:0;
xue@9:           ene+=(a*a+b*b+a*b)/3;
xue@9:         }
xue@9: 
xue@9:         if (hsfr==hsFr-1) ene+=a*a;
xue@9:         else
xue@9:         {
xue@9:           double b=(hsPartials[m][hsfr+1].f>0)?hsPartials[m][hsfr+1].a:0;
xue@9:           ene+=(a*a+b*b+a*b)/3;
xue@9:         }
xue@9:       }
xue@9:       framesa[fr1+hsfr]=ene;//*2;
xue@9:       if (framesa[fr1+hsfr]<frames[fr1+hsfr]) consi*=frames[fr1+hsfr]/framesa[fr1+hsfr];
xue@9:       else consi*=framesa[fr1+hsfr]/frames[fr1+hsfr];
xue@9:       if (framesa[fr1+hsfr]>frames[fr1+hsfr]) framesa[fr1+hsfr]=0;
xue@9:       else framesa[fr1+hsfr]=1-framesa[fr1+hsfr]/frames[fr1+hsfr];
xue@9:     }
xue@9:     consi=pow(consi, 1.0/hsFr);
xue@9: 
xue@9:     THS* HS=new THS;
xue@9:     HS->M=hsM; HS->Fr=hsFr; HS->Partials=hsPartials;
xue@9:     printf("%d (%.2fs), %d (%.2fs), inconsistency %f\n", start, start*offst*1.0/sps, hsFr, hsFr*offst*1.0/sps, consi-1);
xue@9:     if (consi<1.5)
xue@9:     {
xue@9:       HS->WriteToStream(File);
xue@9:       printf("Write note %d\n", writecount++);
xue@9:     }
xue@9:     delete HS;
xue@9: 
xue@9:     start=fr1+hsFr;
xue@9: 
xue@9:   }
xue@9:   delete File; delete Sp;
xue@9:   delete List;
xue@9:   printf("Completed\n"); getch(); return 0;
xue@9: }