diff multires.cpp @ 1:6422640a802f

first upload

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