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 <math.h>
xue@1: #include "multires.h"
xue@1: #include "arrayalloc.h"
xue@1: #include "procedures.h"
xue@1: 
Chris@5: /** \file multires.h */
Chris@5: 
xue@1: //---------------------------------------------------------------------------
xue@1: 
xue@1: //function xlogx(x): returns x*log(x)
xue@1: inline double xlogx(double x)
xue@1: {
xue@1:   if (x==0) return 0;
xue@1:   else return x*log(x);
xue@1: }//xlogx
xue@1: 
xue@1: //macro NORMAL_: a normalization step used for tiling
xue@1: #define NORMAL_(A, a) A=a*(A+log(a));
xue@1: //  #define NORMAL_(A, a) A=a*a*A;
xue@1: //  #define NORMAL_(A, a) A=sqrt(a)*A;
xue@1: 
Chris@5: /**
xue@1:   function DoCutSpectrogramSquare: find optimal tiling of a square. This is a recursive procedure.
xue@1: 
xue@1:   In: Specs[0][1][N], Specs[1][2][N/2], ..., Specs[log2(N)][N][1], multiresolution power spectrogram
xue@1:       e[Res]: total power of each level, e[i] equals the sum of Specs[i][][]
xue@1:       NN: maximal tile height
xue@1:   Out: cuts[N-1]: the tiling result
xue@1:        ents[Res]
xue@1: 
xue@1:   Returns the entropy of the output tiling.
xue@1: */
xue@1: double DoCutSpectrogramSquare(int* cuts, double*** Specs, double* e, int N, int NN, bool Norm, double* ents)
xue@1: {
xue@1:   double result;
xue@1:   int Res=log2(N)+1;
xue@1: 
xue@1:   if (N==1) // 1*1 only(no cuts), returns the sample function.
xue@1:   {
xue@1:     double sp00;
xue@1:     if (e[0]!=0) sp00=Specs[0][0][0]/e[0]; else sp00=0;
xue@1:     ents[0]=xlogx(sp00);
xue@1:     return ents[0];
xue@1:   }
xue@1:   else if (N==2)
xue@1:   {
xue@1:     double sp00, sp01, sp10, sp11;
xue@1:     if (e[0]!=0) sp00=Specs[0][0][0]/e[0], sp01=Specs[0][0][1]/e[0]; else sp00=sp01=0;
xue@1:     if (e[1]!=0) sp10=Specs[1][0][0]/e[1], sp11=Specs[1][1][0]/e[1]; else sp10=sp11=0;
xue@1:     double ent0=xlogx(sp00)+xlogx(sp01);
xue@1:     double ent1=xlogx(sp10)+xlogx(sp11);
xue@1:     if (ent0<ent1)
xue@1:     {
xue@1:       cuts[0]=1;
xue@1:       ents[0]=0, ents[1]=ent1;
xue@1:     }
xue@1:     else
xue@1:     {
xue@1:       cuts[0]=0;
xue@1:       ents[0]=ent0, ents[1]=0;
xue@1:     }
xue@1:   }
xue@1:   else
xue@1:   {
xue@1:     int* tmpcuts=new int[N-2];
xue@1:     int *lcuts, *rcuts;
xue@1:     double ***lSpecs, ***rSpecs, *el, *er, ent0, ent1, a;
xue@1:     double *entl0=new double[Res-1], *entr0=new double[Res-1],
xue@1:            *entl1=new double[Res-1], *entr1=new double[Res-1];
xue@1:     //vertical cuts: l->left half, r->right half
xue@1:     if (N<=NN)
xue@1:     {
xue@1:       lcuts=&cuts[1], rcuts=&cuts[N/2];
xue@1:       VSplitSpecs(N, Specs, lSpecs, rSpecs);
xue@1:       el=new double[Res-1], er=new double[Res-1];
xue@1:       memset(el, 0, sizeof(double)*(Res-1)); memset(er, 0, sizeof(double)*(Res-1));
xue@1:       if (Norm)
xue@1:       {
xue@1:         //normalization
xue@1:         for (int i=0, Fr=1, n=N/2; i<Res-1; i++, Fr*=2, n/=2)
xue@1:         for (int j=0; j<Fr; j++) for (int k=0; k<n; k++)
xue@1:           el[i]+=lSpecs[i][j][k], er[i]+=rSpecs[i][j][k];
xue@1:       }
xue@1:       else
xue@1:         for (int i=0; i<Res-1; i++) el[i]=er[i]=1;
xue@1: 
xue@1:       DoCutSpectrogramSquare(lcuts, lSpecs, el, N/2, NN, Norm, entl1);
xue@1:       DoCutSpectrogramSquare(rcuts, rSpecs, er, N/2, NN, Norm, entr1);
xue@1:                            
xue@1:       ent1=0;
xue@1: 
xue@1:       for (int i=0; i<Res-1; i++)
xue@1:       {
xue@1:         if (e[i]!=0)
xue@1:         {
xue@1:           a=el[i]/e[i]; if (a>0) {NORMAL_(entl1[i], a);} else entl1[i]=0; ent1=ent1+entl1[i];
xue@1:           a=er[i]/e[i]; if (a>0) {NORMAL_(entr1[i], a);} else entr1[i]=0; ent1=ent1+entr1[i];
xue@1:         }
xue@1:         else
xue@1:           entl1[i]=entr1[i]=0;
xue@1:       }
xue@1: 
xue@1:       DeAlloc2(lSpecs); DeAlloc2(rSpecs);
xue@1:       delete[] el; delete[] er;
xue@1: 
xue@1:     }
xue@1:     //horizontal cuts: l->lower half, r->upper half
xue@1:     lcuts=tmpcuts, rcuts=&tmpcuts[N/2-1];
xue@1:     HSplitSpecs(N, Specs, lSpecs, rSpecs);
xue@1:     el=new double[Res-1], er=new double[Res-1];
xue@1:     memset(el, 0, sizeof(double)*(Res-1)); memset(er, 0, sizeof(double)*(Res-1));
xue@1:     if (Norm)
xue@1:     {
xue@1:       //normalization
xue@1:       for (int i=0, Fr=1, n=N/2; i<Res-1; i++, Fr*=2, n/=2)
xue@1:       for (int j=0; j<Fr; j++) for (int k=0; k<n; k++)
xue@1:         el[i]+=lSpecs[i][j][k], er[i]+=rSpecs[i][j][k];
xue@1:     }
xue@1:     else
xue@1:       for (int i=0; i<Res-1; i++) el[i]=er[i]=1;
xue@1: 
xue@1:     DoCutSpectrogramSquare(lcuts, lSpecs, el, N/2, NN, Norm, entl0);
xue@1:     DoCutSpectrogramSquare(rcuts, rSpecs, er, N/2, NN, Norm, entr0);
xue@1: 
xue@1:     ent0=0;
xue@1:     
xue@1:     if (Norm)
xue@1:     for (int i=0; i<Res-1; i++)
xue@1:     {
xue@1:       if (e[i]!=0)
xue@1:       {
xue@1:         a=el[i]/e[i]; if (a>0) {NORMAL_(entl0[i], a);} else entl0[i]=0; ent0=ent0+entl0[i];
xue@1:         a=er[i]/e[i]; if (a>0) {NORMAL_(entr0[i], a);} else entr0[i]=0; ent0=ent0+entr0[i];
xue@1:       }
xue@1:       else
xue@1:         entl0[i]=entr0[i]=0;
xue@1:     }
xue@1: 
xue@1:     DeAlloc2(lSpecs); DeAlloc2(rSpecs);
xue@1:     delete[] el; delete[] er;
xue@1: 
xue@1:     if (N<=NN && ent0<ent1)
xue@1:     {
xue@1:       cuts[0]=1;
xue@1:       result=ent1;
xue@1:       for (int i=0; i<Res-1; i++)
xue@1:       {
xue@1:         ents[i+1]=entl1[i]+entr1[i];
xue@1:       }
xue@1:       ents[0]=0;
xue@1:     }
xue@1:     else
xue@1:     {
xue@1:       memcpy(&cuts[1], tmpcuts, sizeof(int)*(N-2));
xue@1:       cuts[0]=0;
xue@1:       result=ent0;
xue@1:       for (int i=0; i<Res-1; i++)
xue@1:       {
xue@1:         ents[i]=entl0[i]+entr0[i];
xue@1:       }
xue@1:       ents[Res-1]=0;
xue@1:     }
xue@1: 
xue@1:     delete[] tmpcuts; 
xue@1:     delete[] entl0; delete[] entl1; delete[] entr0; delete[] entr1;
xue@1:   }
xue@1: 
xue@1:   return result;
xue@1: }//DoCutSpectrogramSquare
xue@1: 
Chris@5: /**
xue@1:   function DoMixSpectrogramSquare: renders a composite spectrogram on a pixel grid. This is a recursive
xue@1:   procedure.
xue@1: 
xue@1:   In: Specs[0][1][N], Specs[1][2][N/2], Specs[2][4][N/4], ..., Specs[][N][1]: multiresolution power
xue@1:         spectrogram
xue@1:       cuts[N-1]: tiling
xue@1:       X, Y: dimensions of pixel grid to render
xue@1:   Out: Spec[X][Y]: pixel grid rendered to represent the given spectrograms and tiling
xue@1: 
xue@1:   No return value;
xue@1: */
xue@1: void DoMixSpectrogramSquare(double** Spec, int* cuts, double*** Specs, int N, bool Norm, int X=0, int Y=0)
xue@1: {
xue@1:   if (X==0 && Y==0) X=Y=N;
xue@1: 
xue@1:   if (N==1)
xue@1:   {
xue@1:     double value=Specs[0][0][0];//sqrt(Specs[0][0][0]);
xue@1:     value=value;
xue@1:     for (int x=0; x<X; x++) for (int y=0; y<Y; y++) Spec[x][y]=value;
xue@1:   }
xue@1:   else
xue@1:   {
xue@1:     double* e;
xue@1:     int Res;
xue@1: 
xue@1:     if (Norm)
xue@1:     {
xue@1:       //normalization
xue@1:       Res=log2(N)+1;
xue@1:       e=new double[Res];
xue@1:       memset(e, 0, sizeof(double)*Res);
xue@1:       for (int i=0, Fr=1, n=N; i<Res; i++, Fr*=2, n/=2)
xue@1:         for (int j=0; j<Fr; j++)
xue@1:           for (int k=0; k<n; k++)
xue@1:             e[i]+=Specs[i][j][k];
xue@1:       double em=e[0];
xue@1:       for (int i=1; i<Res; i++)
xue@1:       {
xue@1:         if (e[i]>em) e[i]=em/e[i];
xue@1:         else e[i]=1;
xue@1:         if (e[i]>em) em=e[i];
xue@1:       } e[0]=1;                
xue@1:       for (int i=0, Fr=1, n=N; i<Res; i++, Fr*=2, n/=2)
xue@1:       {
xue@1:         if (e[i]!=0 && e[1]!=1)
xue@1:           for (int j=0; j<Fr; j++)
xue@1:             for (int k=0; k<n; k++)
xue@1:               Specs[i][j][k]*=e[i];
xue@1:       }
xue@1:     }
xue@1: 
xue@1:     double **lSpec, **rSpec, ***lSpecs, ***rSpecs;
xue@1:     if (cuts[0]) //1: vertical split
xue@1:     {
xue@1:       VSplitSpecs(N, Specs, lSpecs, rSpecs);
xue@1:       VSplitSpec(X, Y, Spec, lSpec, rSpec);
xue@1:       DoMixSpectrogramSquare(lSpec, &cuts[1], lSpecs, N/2, Norm, X/2, Y);
xue@1:       DoMixSpectrogramSquare(rSpec, &cuts[N/2], rSpecs, N/2, Norm, X/2, Y);
xue@1:     }
xue@1:     else //0: horizontal split
xue@1:     {
xue@1:       HSplitSpecs(N, Specs, lSpecs, rSpecs);
xue@1:       HSplitSpec(X, Y, Spec, lSpec, rSpec);
xue@1:       DoMixSpectrogramSquare(lSpec, &cuts[1], lSpecs, N/2, Norm, X, Y/2);
xue@1:       DoMixSpectrogramSquare(rSpec, &cuts[N/2], rSpecs, N/2, Norm, X, Y/2);
xue@1:     }
xue@1: 
xue@1:     if (Norm)
xue@1:     {
xue@1:       for (int i=0, Fr=1, n=N; i<Res; i++, Fr*=2, n/=2)
xue@1:       {
xue@1:         if (e[i]!=0 && e[1]!=1)
xue@1:           for (int j=0; j<Fr; j++)
xue@1:             for (int k=0; k<n; k++)
xue@1:               Specs[i][j][k]/=e[i];
xue@1:       }
xue@1:       delete[] e;
xue@1:     }
xue@1: 
xue@1:     delete[] lSpec; delete[] rSpec; DeAlloc2(lSpecs); DeAlloc2(rSpecs);
xue@1:   }
xue@1: }//DoMixSpectrogramSquare
xue@1: 
Chris@5: /**
xue@1:   function DoMixSpectrogramSquare: retrieves a composite spectrogram as a vector. This is a recursive
xue@1:   procedure.
xue@1: 
xue@1:   In: Specs[0][1][N], Specs[1][2][N/2], Specs[2][4][N/4], ..., Specs[][N][1]: multiresolution power
xue@1:         spectrogram
xue@1:       cuts[N-1]: tiling
xue@1:   Out: Spec[N]: composite spectrogram sampled fron Specs according to tiling cut[]
xue@1: 
xue@1:   No return value;
xue@1: */
xue@1: void DoMixSpectrogramSquare(double* Spec, int* cuts, double*** Specs, int N, bool Norm)
xue@1: {
xue@1: //  if (X==0 && Y==0) X=Y=N;
xue@1: 
xue@1:   if (N==1)
xue@1:     Spec[0]=Specs[0][0][0];//sqrt(Specs[0][0][0]);
xue@1:   else
xue@1:   {
xue@1:     double* e;
xue@1:     int Res;
xue@1: 
xue@1:     //Norm=false;
xue@1:     if (Norm)
xue@1:     {
xue@1:       //normalization
xue@1:       Res=log2(N)+1;
xue@1:       e=new double[Res];
xue@1:       memset(e, 0, sizeof(double)*Res);
xue@1:       for (int i=0, Fr=1, n=N; i<Res; i++, Fr*=2, n/=2)
xue@1:         for (int j=0; j<Fr; j++)
xue@1:           for (int k=0; k<n; k++)
xue@1:             e[i]+=Specs[i][j][k];
xue@1:       double em=e[0];
xue@1:       for (int i=1; i<Res; i++)
xue@1:       {
xue@1:         if (e[i]>em) e[i]=em/e[i];
xue@1:         else e[i]=1;
xue@1:         if (e[i]>em) em=e[i];
xue@1:       } e[0]=1;
xue@1:       for (int i=0, Fr=1, n=N; i<Res; i++, Fr*=2, n/=2)
xue@1:       {
xue@1:         if (e[i]!=0 && e[i]!=1)
xue@1:           for (int j=0; j<Fr; j++)
xue@1:             for (int k=0; k<n; k++)
xue@1:               Specs[i][j][k]*=e[i];
xue@1:       }
xue@1:     }
xue@1: 
xue@1:     double ***lSpecs, ***rSpecs;
xue@1:     if (cuts[0]) //1: vertical split
xue@1:     {
xue@1:       VSplitSpecs(N, Specs, lSpecs, rSpecs);
xue@1:       DoMixSpectrogramSquare(Spec, &cuts[1], lSpecs, N/2, Norm);
xue@1:       DoMixSpectrogramSquare(&Spec[N/2], &cuts[N/2], rSpecs, N/2, Norm);
xue@1:     }
xue@1:     else //0: horizontal split
xue@1:     {
xue@1:       HSplitSpecs(N, Specs, lSpecs, rSpecs);
xue@1:       DoMixSpectrogramSquare(Spec, &cuts[1], lSpecs, N/2, Norm);
xue@1:       DoMixSpectrogramSquare(&Spec[N/2], &cuts[N/2], rSpecs, N/2, Norm);
xue@1:     }
xue@1: 
xue@1:     if (Norm)
xue@1:     {
xue@1:       for (int i=0, Fr=1, n=N; i<Res; i++, Fr*=2, n/=2)
xue@1:       {
xue@1:         if (e[i]!=0 && e[1]!=1)
xue@1:           for (int j=0; j<Fr; j++)
xue@1:             for (int k=0; k<n; k++)
xue@1:               Specs[i][j][k]/=e[i];
xue@1:       }
xue@1:       delete[] e;
xue@1:     }
xue@1: 
xue@1:     DeAlloc2(lSpecs); DeAlloc2(rSpecs);
xue@1:   }
xue@1: }//DoMixSpectrogramSquare
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function HSplitSpec: split a spectrogram horizontally into lower and upper halves.
xue@1: 
xue@1:   In: Spec[X][Y]: spectrogram to split
xue@1:   Out: lSpec[X][Y/2], uSpec[X][Y/2]: the two half spectrograms
xue@1: 
xue@1:   No return value. Both lSpec and uSpec are allocated anew. The caller is responsible to free these
xue@1:   buffers.
xue@1: */
xue@1: void HSplitSpec(int X, int Y, double** Spec, double**& lSpec, double**& uSpec)
xue@1: {
xue@1:   lSpec=new double*[X], uSpec=new double*[X];
xue@1:   for (int i=0; i<X; i++)
xue@1:     lSpec[i]=Spec[i], uSpec[i]=&Spec[i][Y/2];
xue@1: }//HSplitSpec
xue@1: 
Chris@5: /**
xue@1:   function HSplitSpecs: split a multiresolution spectrogram horizontally into lower and upper halves
xue@1: 
xue@1:   A full spectrogram array is given in log2(N)+1 spectrograms, with the base spec of 1*N, 1st octave of
xue@1:   2*(N/2), ..., last octave of N*1. When this array is split into two spectrogram arrays horizontally,
xue@1:   the last spec (with the highest time resolution). Each of the two new arrays is given in log2(N)
xue@1:   spectrograms.
xue@1: 
xue@1:   In: Specs[nRes+1][][]: multiresolution spectrogram
xue@1:   Out: lSpecs[nRes][][], uSpecs[nRes][][], the two half multiresolution spectrograms
xue@1: 
xue@1:   This function allocates two 2nd order arrays of double*, which the caller is responsible to free.
xue@1: */
xue@1: void HSplitSpecs(int N, double*** Specs, double***& lSpecs, double***& uSpecs)
xue@1: {
xue@1:   int nRes=log2(N); // new number of resolutions
xue@1:   lSpecs=new double**[nRes], uSpecs=new double**[nRes];
xue@1:   lSpecs[0]=new double*[nRes*N/2], uSpecs[0]=new double*[nRes*N/2]; for (int i=1; i<nRes; i++) lSpecs[i]=&lSpecs[0][i*N/2], uSpecs[i]=&uSpecs[0][i*N/2];
xue@1:   for (int i=0, Fr=1, n=N; i<nRes; i++, Fr*=2, n/=2)
xue@1:     for (int j=0; j<Fr; j++) lSpecs[i][j]=Specs[i][j], uSpecs[i][j]=&Specs[i][j][n/2];
xue@1: }//HSplitSpecs
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function MixSpectrogramSquare: find composite spectrogram from multiresolution spectrogram,output as
xue@1:   pixel grid
xue@1: 
xue@1:   In: Specs[0][1][N], Specs[1][2][N/2], ..., Specs[Res-1][N][1], multiresolution spectrogram
xue@1:   Out: Spec[N][N]: composite spectrogram as pixel grid (N time redundant)
xue@1:        cuts[N-1] (optional): tiling
xue@1: 
xue@1:   Returns entropy.
xue@1: */
xue@1: double MixSpectrogramSquare(double** Spec, double*** Specs, int N, int Res, bool Norm, bool NormMix, int* cuts=0)
xue@1: {
xue@1:   int tRes=log2(N)+1;
xue@1:   if (Res==0) Res=tRes;
xue@1:   int NN=1<<(Res-1);
xue@1:   bool createcuts=(cuts==0);
xue@1:   if (createcuts) cuts=new int[N];
xue@1:   double* e=new double[tRes], *ents=new double[tRes];
xue@1:   //normalization
xue@1:   memset(e, 0, sizeof(double)*Res);
xue@1:   for (int i=0, Fr=1, n=N; i<tRes; i++, Fr*=2, n/=2)
xue@1:     for (int j=0; j<Fr; j++)
xue@1:       for (int k=0; k<n; k++)
xue@1:         e[i]+=Specs[i][j][k];
xue@1: 
xue@1:   if (!Norm)
xue@1:   {
xue@1:     for (int i=0, Fr=1, n=N; i<tRes; i++, Fr*=2, n/=2)
xue@1:     {
xue@1:       if (e[i]!=0 && e[i]!=1)
xue@1:         for (int j=0; j<Fr; j++)
xue@1:           for (int k=0; k<n; k++)
xue@1:             Specs[i][j][k]/=e[i];
xue@1:     }
xue@1: //    for (int i=0; i<Res; i++) e[i]=1;
xue@1:   }
xue@1: 
xue@1:   double result=DoCutSpectrogramSquare(cuts, Specs, e, N, NN, Norm, ents);
xue@1:   delete[] ents;
xue@1: 
xue@1:   if (!Norm)
xue@1:   {
xue@1:     for (int i=0, Fr=1, n=N; i<tRes; i++, Fr*=2, n/=2)
xue@1:     {
xue@1:       if (e[i]!=0 && e[i]!=1)
xue@1:         for (int j=0; j<Fr; j++)
xue@1:           for (int k=0; k<n; k++)
xue@1:             Specs[i][j][k]*=e[i];
xue@1:     }
xue@1:   }
xue@1:   DoMixSpectrogramSquare(Spec, cuts, Specs, N, NormMix, N, N);
xue@1: 
xue@1:   delete[] e;
xue@1:   if (createcuts) delete[] cuts;
xue@1:   return result;
xue@1: }//MixSpectrogramSquare
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function MixSpectrogramSquare: find composite spectrogram from multiresolution spectrogram,output as
xue@1:   vectors
xue@1: 
xue@1:   In: Specs[0][1][N], Specs[1][2][N/2], ..., Specs[Res-1][N][1], multiresolution spectrogram
xue@1:   Out: spl[N-1], Spec[N]: composite spectrogram as tiling and value vectors
xue@1: 
xue@1:   Returns entropy.
xue@1: */
xue@1: double MixSpectrogramSquare(int* spl, double* Spec, double*** Specs, int N, int Res, bool Norm, bool NormMix)
xue@1: {
xue@1:   int tRes=log2(N)+1;
xue@1:   if (Res==0) Res=tRes;
xue@1:   int NN=1<<(Res-1);
xue@1: 
xue@1:   double* e=new double[tRes], *ents=new double[tRes];
xue@1: 
xue@1:   //normalization
xue@1:   memset(e, 0, sizeof(double)*Res);
xue@1:   for (int i=0, Fr=1, n=N; i<tRes; i++, Fr*=2, n/=2)
xue@1:     for (int j=0; j<Fr; j++)
xue@1:       for (int k=0; k<n; k++)
xue@1:         e[i]+=Specs[i][j][k];
xue@1: 
xue@1:   if (!Norm)
xue@1:   {
xue@1:     for (int i=0, Fr=1, n=N; i<tRes; i++, Fr*=2, n/=2)
xue@1:     {
xue@1:       if (e[i]!=0 && e[i]!=1)
xue@1:         for (int j=0; j<Fr; j++)
xue@1:           for (int k=0; k<n; k++)
xue@1:             Specs[i][j][k]/=e[i];
xue@1:     }
xue@1: //    for (int i=0; i<Res; i++) e[i]=1;
xue@1:   }
xue@1: 
xue@1:   double result=DoCutSpectrogramSquare(spl, Specs, e, N, NN, Norm, ents);
xue@1:   delete[] ents;
xue@1:   if (!Norm)
xue@1:   {
xue@1:     for (int i=0, Fr=1, n=N; i<tRes; i++, Fr*=2, n/=2)
xue@1:     {
xue@1:       if (e[i]!=0 && e[i]!=1)
xue@1:         for (int j=0; j<Fr; j++)
xue@1:           for (int k=0; k<n; k++)
xue@1:             Specs[i][j][k]*=e[i];
xue@1:     }
xue@1:   }
xue@1:   DoMixSpectrogramSquare(Spec, spl, Specs, N, NormMix);
xue@1:   return result;
xue@1: }//MixSpectrogramSquare
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function MixSpectrogramBlock: obtain the composite spectrogram of a waveform block as pixel grid.
xue@1: 
xue@1:   This method deals with the effective duration of WID/2 samples of a frame of WID samples. The maximal
xue@1:   frame width is WID, minimal frame width is wid. The spectrum with frame width WID (base) is given in
xue@1:   lSpecs[0][0], the spectra with frame width WID/2 (1st octave) is given in lSpecs[1][0] & lSpecs[1][1],
xue@1:   etc.
xue@1: 
xue@1:   The output Spec[WID/wid][WID] is a spectrogram sampled from lSpecs, with a redundancy factor WID/wid.
xue@1:   The sampling is optimized so that a maximal entropy is achieved globally. This maximal entropy is
xue@1:   returned.
xue@1: 
xue@1:   In: Specs[0][1][WID], Specs[1][2][WID/2], ..., Specs[Res-1][WID/wid][wid], multiresolution spectrogram
xue@1:   Out: Spec[WID/wid][WID], composite spectrogram as pixel grid
xue@1:        cuts[wid][WID/wid-1] (optional), tilings of the wid squares the block is divided into.
xue@1: 
xue@1:   Returns entropy.
xue@1: */
xue@1: double MixSpectrogramBlock(double** Spec, double*** Specs, int WID, int wid, bool norm, bool normmix, int** cuts=0)
xue@1: {
xue@1:   int N=WID/wid, Res=log2(WID/wid)+1;
xue@1:   double result=0, **lSpec=new double*[N], ***lSpecs=new double**[Res];
xue@1:   lSpecs[0]=new double*[Res*N]; for (int i=1; i<Res; i++) lSpecs[i]=&lSpecs[0][i*N];
xue@1: 
xue@1:   bool createcuts=(cuts==0);
xue@1:   if (createcuts)
xue@1:   {
xue@1:     cuts=new int*[wid];
xue@1:     memset(cuts, 0, sizeof(int*)*wid);
xue@1:   }
xue@1: 
xue@1:   for (int i=0; i<wid; i++)
xue@1:   {
xue@1:     for (int j=0; j<N; j++)
xue@1:       lSpec[j]=&Spec[j][i*N];
xue@1:     for (int j=0, n=N, Fr=1; j<Res; j++, n/=2, Fr*=2)
xue@1:     {
xue@1:       for (int k=0; k<Fr; k++)
xue@1:         lSpecs[j][k]=&Specs[j][k][i*n];
xue@1:     }
xue@1: 
xue@1:     int Log2N=log2(N);
xue@1:     if (i==0)
xue@1:       result+=MixSpectrogramSquare(lSpec, lSpecs, N, Log2N>2?(log2(N)-1):2, norm, normmix, cuts[i]);
xue@1:     else if (i==1)
xue@1:       result+=MixSpectrogramSquare(lSpec, lSpecs, N, Log2N>1?Log2N:2, norm, normmix, cuts[i]);
xue@1:     else
xue@1:       result+=MixSpectrogramSquare(lSpec, lSpecs, N, Log2N+1, norm, normmix, cuts[i]);
xue@1:   }
xue@1:   delete[] lSpec;
xue@1:   DeAlloc2(lSpecs);
xue@1:   if (createcuts) delete[] cuts;
xue@1:   return result;
xue@1: }//MixSpectrogramBlock
xue@1: 
Chris@5: /**
xue@1:   function MixSpectrogramBlock: obtain the composite spectrogram of a waveform block as vectors.
xue@1: 
xue@1:   In: Specs[0][1][WID], Specs[1][2][WID/2], ..., Specs[Res-1][WID/wid][wid], multiresolution spectrogram
xue@1:   Out: spl[WID], Spec[WID], composite spectrogram as tiling and value vectors. Each of the vectors is
xue@1:        made up of $wid subvectors, each subvector pair describing a N*N block of the composite
xue@1:        spectrogram.
xue@1: 
xue@1:   Returns entropy.
xue@1: */
xue@1: double MixSpectrogramBlock(int* spl, double* Spec, double*** Specs, int WID, int wid, bool norm, bool normmix)
xue@1: {
xue@1:   int N=WID/wid, Res=log2(WID/wid)+1, *lspl;
xue@1:   double result=0, *lSpec, ***lSpecs=new double**[Res];
xue@1:   lSpecs[0]=new double*[Res*N]; for (int i=1; i<Res; i++) lSpecs[i]=&lSpecs[0][i*N];
xue@1: 
xue@1:   for (int i=0; i<wid; i++)
xue@1:   {
xue@1:     lspl=&spl[i*N];
xue@1:     lSpec=&Spec[i*N];
xue@1:     for (int j=0, n=N, Fr=1; j<Res; j++, n/=2, Fr*=2)
xue@1:     {
xue@1:       for (int k=0; k<Fr; k++)
xue@1:         lSpecs[j][k]=&Specs[j][k][i*n];
xue@1:     }
xue@1:     int Log2N=log2(N);
xue@1:     /*
xue@1:     if (i==0)
xue@1:       result+=MixSpectrogramSquare(lspl, lSpec, lSpecs, N, Log2N>2?(log2(N)-1):2, norm, normmix);
xue@1:     else if (i==1)
xue@1:       result+=MixSpectrogramSquare(lspl, lSpec, lSpecs, N, Log2N>1?Log2N:2, norm, normmix);
xue@1:     else     */
xue@1:       result+=MixSpectrogramSquare(lspl, lSpec, lSpecs, N, Log2N+1, norm, normmix);
xue@1: 
xue@1:   }
xue@1:   DeAlloc2(lSpecs);
xue@1: 
xue@1:   return result;
xue@1: }//MixSpectrogramBlock
xue@1: 
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
Chris@5:   functions names as ...Block2(...) implement the same functions as the above directly without
xue@1:   explicitly dividing the multiresolution spectrogram into square blocks.
xue@1: */
xue@1: 
Chris@5: /**
xue@1:   function DoCutSpectrogramBlock2: find optimal tiling for a block
xue@1: 
xue@1:   In: Specs[R0][x0:x0+x-1][Y0:Y0+Y-1], [R0+1][2x0:2x0+2x-1][Y0/2:Y0/2+Y/2-1],...,
xue@1:         Specs[R0+?][Nx0:Nx0+Nx-1][Y0/N:Y0/N+Y/N-1], multiresolution spectrogram
xue@1:   Out: spl[Y-1], tiling of this block
xue@1: 
xue@1:   Returns entropy.
xue@1: */
xue@1: double DoCutSpectrogramBlock2(int* spl, double*** Specs, int Y, int R0, int x0, int Y0, int N, double& ene)
xue@1: {
xue@1:   double ent=0;
xue@1:   if (Y>N) //N=WID/wid, the actual square size
xue@1:   {
xue@1:     spl[0]=0;
xue@1:     double ene1, ene2;
xue@1:     ent+=DoCutSpectrogramBlock2(&spl[1], Specs, Y/2, R0, x0, Y0, N, ene1);
xue@1:     ent+=DoCutSpectrogramBlock2(&spl[Y/2], Specs, Y/2, R0, x0, Y0+Y/2, N, ene2);
xue@1:     ene=ene1+ene2;
xue@1:   }
xue@1:   else if (N==1)
xue@1:   {
xue@1:     double tmp=Specs[R0][x0][Y0];
xue@1:     ene=tmp;
xue@1:     ent=xlogx(tmp);
xue@1:   }
xue@1:   else //Y==N, the square case
xue@1:   {
xue@1:     double enel, ener, enet, eneb, entl, entr, entt, entb;
xue@1:     int* tmpspl=new int[Y];
xue@1:     entl=DoCutSpectrogramBlock2(&spl[1], Specs, Y/2, R0+1, 2*x0, Y0/2, N/2, enel);
xue@1:     entr=DoCutSpectrogramBlock2(&spl[Y/2], Specs, Y/2, R0+1, 2*x0+1, Y0/2, N/2, ener);
xue@1:     entb=DoCutSpectrogramBlock2(&tmpspl[1], Specs, Y/2, R0, x0, Y0, N/2, eneb);
xue@1:     entt=DoCutSpectrogramBlock2(&tmpspl[Y/2], Specs, Y/2, R0, x0, Y0+Y/2, N/2, enet);
xue@1:     double ene0=enet+eneb, ene1=enel+ener, ent0=entt+entb, ent1=entl+entr;
xue@1:     //normalize
xue@1:     double eneres=1-(ene0+ene1)/2, norment0, norment1;
xue@1:     //double a0=1/(ene0+eneres), a1=1/(ene1+eneres);
xue@1:     //quasi-global normalization
xue@1:     norment0=(ent0-ene0*log(ene0+eneres))/(ene0+eneres), norment1=(ent1-ene1*log(ene1+eneres))/(ene1+eneres);
xue@1:     //local normalization
xue@1:     //if (ene0>0) norment0=ent0/ene0-log(ene0); else norment0=0; if (ene1>0) norment1=ent1/ene1-log(ene1); else norment1=0;
xue@1:     if (norment1<norment0)
xue@1:     {
xue@1:       spl[0]=0;
xue@1:       ent=ent0, ene=ene0;
xue@1:       memcpy(&spl[1], &tmpspl[1], sizeof(int)*(Y-2));
xue@1:     }
xue@1:     else
xue@1:     {
xue@1:       spl[0]=1;
xue@1:       ent=ent1, ene=ene1;
xue@1:     }
xue@1:   }
xue@1:   return ent;
xue@1: }//DoCutSpectrogramBlock2
xue@1: 
Chris@5: /**
xue@1:   function DoMixSpectrogramBlock2: sampling multiresolution spectrogram according to given tiling
xue@1: 
xue@1:   In: Specs[R0][x0:x0+x-1][Y0:Y0+Y-1], [R0+1][2x0:2x0+2x-1][Y0/2:Y0/2+Y/2-1],...,
xue@1:         Specs[R0+?][Nx0:Nx0+Nx-1][Y0/N:Y0/N+Y/N-1], multiresolution spectrogram
xue@1:       spl[Y-1]; tiling of this block
xue@1:   Out: Spec[Y], composite spectrogram as value vector
xue@1: 
xue@1:   Returns 0.
xue@1: */
xue@1: double DoMixSpectrogramBlock2(int* spl, double* Spec, double*** Specs, int Y, int R0, int x0, int Y0, bool normmix, int res, double* e)
xue@1: {
xue@1:   if (Y==1)
xue@1:   {
xue@1:     Spec[0]=Specs[R0][x0][Y0]*e[0];
xue@1:   }
xue@1:   else
xue@1:   {
xue@1:     double le[32];
xue@1:     if (normmix && Y<(1<<res))
xue@1:     {
xue@1:       for (int i=0, j=1, k=Y; i<res; i++, j*=2, k/=2)
xue@1:       {
xue@1:         double lle=0;
xue@1:         for (int fr=0; fr<j; fr++) for (int n=0; n<k; n++) lle+=Specs[i+R0][x0+fr][Y0+n]*Specs[i+R0][x0+fr][Y0+n];
xue@1:         lle=sqrt(lle)*e[i];
xue@1:         if (i==0) le[0]=lle;
xue@1:         else if (lle>le[0]*2) le[i]=e[i]*le[0]*2/lle;
xue@1:         else le[i]=e[i];
xue@1:       }
xue@1:       le[0]=e[0];
xue@1:     }
xue@1:     else
xue@1:     {
xue@1:       memcpy(le, e, sizeof(double)*res);
xue@1:     }
xue@1: 
xue@1:     if (spl[0]==0)
xue@1:     {
xue@1:       int newres;
xue@1:       if (Y>=(1<<res)) newres=res;
xue@1:       else newres=res-1;
xue@1:       DoMixSpectrogramBlock2(&spl[1], Spec, Specs, Y/2, R0, x0, Y0, normmix, newres, le);
xue@1:       DoMixSpectrogramBlock2(&spl[Y/2], &Spec[Y/2], Specs, Y/2, R0, x0, Y0+Y/2, normmix, newres, le);
xue@1:     }
xue@1:     else
xue@1:     {
xue@1:       DoMixSpectrogramBlock2(&spl[1], Spec, Specs, Y/2, R0+1, x0*2, Y0/2, normmix, res-1, &le[1]);
xue@1:       DoMixSpectrogramBlock2(&spl[Y/2], &Spec[Y/2], Specs, Y/2, R0+1, x0*2+1, Y0/2, normmix, res-1, &le[1]);
xue@1:     }
xue@1:   }
xue@1:   return 0;
xue@1: }//DoMixSpectrogramBlock2
xue@1: 
Chris@5: /**
xue@1:   function MixSpectrogramBlock2: obtain the composite spectrogram of a waveform block as vectors.
xue@1: 
xue@1:   In: Specs[0][1][WID], Specs[1][2][WID/2], ..., Specs[Res-1][WID/wid][wid], multiresolution spectrogram
xue@1:   Out: spl[WID], Spec[WID], composite spectrogram as tiling and value vectors. Each of the
xue@1:        vectors is made up of $wid subvectors, each subvector pair describing a N*N block of the
xue@1:        composite spectrogram.
xue@1: 
xue@1:   Returns entropy.
xue@1: */
xue@1: double MixSpectrogramBlock2(int* spl, double* Spec, double*** Specs, int WID, int wid, bool normmix)
xue@1: {
xue@1:   double ene[32];
xue@1:   //find the total energy and normalize
xue@1:   for (int i=0, Fr=1, Wid=WID; Wid>=wid; i++, Fr*=2, Wid/=2)
xue@1:   {
xue@1:     double lene=0;
xue@1:     for (int fr=0; fr<Fr; fr++) for (int k=0; k<Wid; k++) lene+=Specs[i][fr][k]*Specs[i][fr][k];
xue@1:     ene[i]=lene;
xue@1:     if (lene!=0)
xue@1:     {
xue@1:       double ilene=1.0/lene;
xue@1:       for (int fr=0; fr<Fr; fr++) for (int k=0; k<Wid; k++) Specs[i][fr][k]=Specs[i][fr][k]*Specs[i][fr][k]*ilene;
xue@1:     }
xue@1:   }
xue@1: 
xue@1:   double result=DoCutSpectrogramBlock2(spl, Specs, WID, 0, 0, 0, WID/wid, ene[31]);
xue@1: 
xue@1:   //de-normalize
xue@1:   for (int i=0, Fr=1, Wid=WID; Wid>=wid; i++, Fr*=2, Wid/=2)
xue@1:   {
xue@1:     double lene=ene[i];
xue@1:     if (lene!=0)
xue@1:       for (int fr=0; fr<Fr; fr++) for (int k=0; k<Wid; k++) Specs[i][fr][k]=sqrt(Specs[i][fr][k]*lene);
xue@1:   }
xue@1: 
xue@1:   double e[32]; for (int i=0; i<32; i++) e[i]=1;
xue@1:   DoMixSpectrogramBlock2(spl, Spec, Specs, WID, 0, 0, 0, normmix, log2(WID/wid)+1, e);
xue@1:   return result;
xue@1: }//MixSpectrogramBlock2
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function MixSpectrogram: obtain composite spectrogram from multiresolutin spectrogram as pixel grid
xue@1: 
xue@1:   This method deals with Fr (base) frames of WID samples. Each base frame may be divided into 2 1st-
xue@1:   octave frames, 4 2nd-octave frames, ..., etc. The spectrogram calculated on base frame is given in
xue@1:   Specs[0] (Fr frames); that of 1st octave is given in Specs[1] (2*Fr frames); etc. The method resamples
xue@1:   the spectrograms of different frame width into a single spectrogram so that the entropy is maximized
xue@1:   globally.
xue@1: 
xue@1:   The output Spec is a spectrogram of apparent resolution WID at hop size wid. It is a redundant
xue@1:   representation, with equal values occupying blocks of size WID/wid.
xue@1: 
xue@1:   In: Specs[0][Fr][WID], Specs[1][Fr*2][WID/2], ..., Specs[Res-1] [Fr*(WID/wid)][wid], multiresolution
xue@1:     spectrogram
xue@1:   Out: Spec[Fr*(WID/wid)][WID], composite spectrogram as pixel grid
xue@1:        cuts[Fr][wid][N=Wid/wid], tilings of small square blocks
xue@1:   Returns 0.
xue@1: */
xue@1: double MixSpectrogram(double** Spec, double*** Specs, int Fr, int WID, int wid, bool norm, bool normmix, int*** cuts)
xue@1: {
xue@1:   //totally Fr frames of WID samples
xue@1:   //each frame is divided into wid VERTICAL parts
xue@1:   bool createcuts=(cuts==0);
xue@1:   if (createcuts)
xue@1:   {
xue@1:     cuts=new int**[Fr];
xue@1:     memset(cuts, 0, sizeof(int**)*Fr);
xue@1:   }
xue@1:   int Res=log2(WID/wid)+1;
xue@1:   double*** lSpecs=new double**[Res];
xue@1:   for (int  i=0; i<Fr; i++)
xue@1:   {
xue@1:     for (int j=0, nfr=1; j<Res; j++, nfr*=2) lSpecs[j]=&Specs[j][i*nfr];
xue@1:     MixSpectrogramBlock(&Spec[i*WID/wid], lSpecs, WID, wid, norm, normmix, cuts[i]);
xue@1:   }
xue@1: 
xue@1:   delete[] lSpecs;
xue@1:   if (createcuts) delete[] cuts;
xue@1:   return 0;
xue@1: }//MixSpectrogram
xue@1: 
Chris@5: /**
xue@1:   function MixSpectrogram: obtain composite spectrogram from multiresolutin spectrogram as vectors
xue@1: 
xue@1:   In: Specs[0][Fr][WID], Specs[1][Fr*2][WID/2], ..., Specs[Res-1] [Fr*(WID/wid)][wid], multiresolution
xue@1:         spectrogram
xue@1:   Out: spl[Fr][WID], Spec[Fr][WID], composite spectrogram as tiling and value vectors by frame.
xue@1: 
xue@1:   Returns 0.
xue@1: */
xue@1: double MixSpectrogram(int** spl, double** Spec, double*** Specs, int Fr, int WID, int wid, bool norm, bool normmix)
xue@1: {
xue@1:   //totally Fr frames of WID samples
xue@1:   //each frame is divided into wid VERTICAL parts
xue@1:   int Res=log2(WID/wid)+1;
xue@1:   double*** lSpecs=new double**[Res];
xue@1:   for (int  i=0; i<Fr; i++)
xue@1:   {
xue@1:     for (int j=0, nfr=1; j<Res; j++, nfr*=2) lSpecs[j]=&Specs[j][i*nfr];
xue@1:     MixSpectrogramBlock(spl[i], Spec[i], lSpecs, WID, wid, norm, normmix);
xue@1: //    MixSpectrogramBlock2(spl[i], Spec[i], lSpecs, WID, wid, norm);
xue@1:   }
xue@1: 
xue@1:   delete[] lSpecs;
xue@1:   return 0;
xue@1: }//MixSpectrogram
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function VSplitSpec: split a spectrogram vertically into left and right halves.
xue@1: 
xue@1:   In: Spec[X][Y]: spectrogram to split
xue@1:   Out: lSpec[X][Y/2], rSpec[X][Y/2]: the two half spectrograms
xue@1: 
xue@1:   No return value. Both lSpec and rSpec are allocated anew. The caller is responsible to free these
xue@1:   buffers.
xue@1: */
xue@1: void VSplitSpec(int X, int Y, double** Spec, double**& lSpec, double**& rSpec)
xue@1: {
xue@1:   lSpec=new double*[X/2], rSpec=new double*[X/2];
xue@1:   for(int i=0; i<X/2; i++)
xue@1:     lSpec[i]=Spec[i], rSpec[i]=Spec[i+X/2];
xue@1: }//VSplitSpec
xue@1: 
Chris@5: /**
xue@1:   function VSplitSpecs: split a multiresolution spectrogram vertically into left and right halves
xue@1: 
xue@1:   A full spectrogram array is given in log2(N)+1 spectrograms, with the base spec of 1*N, 1st octave of
xue@1:   2*(N/2), ..., last octave of N*1. When this array is split into two spectrogram arrays horizontally,
xue@1:   the last spec (with the highest time resolution). Each of the two new arrays is given in log2(N)
xue@1:   spectrograms.
xue@1: 
xue@1:   In: Specs[nRes+1][][]: multiresolution spectrogram
xue@1:   Out: lSpecs[nRes][][], rSpecs[nRes][][], the two half multiresolution spectrograms
xue@1: 
xue@1:   This function allocates two 2nd order arrays of double*, which the caller is responsible to free.
xue@1: */
xue@1: void VSplitSpecs(int N, double*** Specs, double***& lSpecs, double***& rSpecs)
xue@1: {
xue@1:   int nRes=log2(N); // new number of resolutions
xue@1:   lSpecs=new double**[nRes], rSpecs=new double**[nRes];
xue@1:   lSpecs[0]=new double*[nRes*N/2], rSpecs[0]=new double*[nRes*N/2]; for (int i=1; i<nRes; i++) lSpecs[i]=&lSpecs[0][i*N/2], rSpecs[i]=&rSpecs[0][i*N/2];
xue@1:   for (int i=0, Fr=1; i<nRes; i++, Fr*=2)
xue@1:     for (int j=0; j<Fr; j++)
xue@1:       lSpecs[i][j]=Specs[i+1][j], rSpecs[i][j]=Specs[i+1][j+Fr];
xue@1: }//VSplitSpecs
xue@1: 
xue@1: