annotate multires.cpp @ 13:de3961f74f30 tip

Add Linux/gcc Makefile; build fix

/*
    Harmonic sinusoidal modelling and tools

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

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

#include <math.h>
#include "multires.h"
#include "arrayalloc.h"
#include "procedures.h"

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