/*
  Whole-Sample Symmetric Lifting: Prediction Step

  wslpredict( u[], N, dq, coeff ):
  [0] Let i=dq
  [1] While i<N-2*dq, do [2] to [3]
  [2]    Accumulate  u[i] += coeff*(u[i-dq] + u[i+dq])
  [3]    Increment  i += 2*dq  
  [4] If i+dq<N, then accumulate u[i] += coeff*(u[i-dq]+u[i+dq])
  [5] Else accumulate u[i] += 2*coeff*u[i-dq]
  [6] Return


  Whole-Sample Symmetric Lifting: Updating Step

  wslupdate( u[], N, dq, coeff ):
  [0] Accumulate u[0] += 2*coeff*u[dq]
  [1] Let i=2*dq
  [2] While i<N-2*dq, do [2] to [3]
  [3]    Accumulate  u[i] += coeff*(u[i-dq] + u[i+dq])
  [4]    Increment  i += 2*dq  
  [5] If i<N, then do [6] to [7]
  [6]    If i+dq<N, then accumulate u[i] += coeff*(u[i-dq]+u[i+dq])
  [7]    Else accumulate u[i] += 2*coeff*u[i-dq]
  [8] Return


  Lifting: Normalization Step

  lnormalize( u[], N, dq, coeff ):
  [0] Let i=0
  [1] While i<N-2*dq, do [2] to [4]
  [2]    Replace  u[i] *= coeff
  [3]    Replace  u[i+dq] /= coeff
  [4]    Increment  i += 2*dq  
  [5] Replace  u[i] *= coeff
  [6] If i+dq<N, then replace u[i+dq] /= coeff
  [7] Return

*/

void
wslpredict(float u[], int N, int dq, float coeff)
{
  int i;
  i = dq;
  while(i<N-2*dq) {
    u[i] += coeff*(u[i-dq]+u[i+dq]);
    i += 2*dq;
  }
  if(i+dq<N) u[i] += coeff*(u[i-dq]+u[i+dq]);
  else u[i] += 2*coeff*u[i-dq];
  return;
}

void
wslupdate(float u[], int N, int dq, float coeff)
{
  int i;
  u[0] += 2*coeff*u[dq];
  i = 2*dq;
  while(i<N-2*dq) {
    u[i] += coeff*(u[i-dq]+u[i+dq]);
    i += 2*dq;
  }
  if(i<N)
    if(i+dq<N) u[i] += coeff*(u[i-dq]+u[i+dq]);
    else u[i] += 2*coeff*u[i-dq];
  return;
}

void
lnormalize(float u[], int N, int dq, float coeff)
{
  int i;
  i = 0;
  while(i<N-2*dq) {
    u[i] *= coeff;
    u[i+dq] /= coeff;
    i += 2*dq;
  }
  u[i] *= coeff;
  if(i+dq<N) u[i] /= coeff;
  return;
}

/*

  Whole-Sample Symmetric Lifting: 4,2 Biorthogonal Filter Transform

  wsl42filter( u[], N, dq ):
  [0] Compute wslpredict ( u[], N, dq, alpha42 )
  [1] Compute wslupdate ( u[], N, dq, beta42 )
  [2] Compute wslpredict ( u[], N, dq, gamma42 )
  [3] Compute lnormalize ( u[], N, dq, zeta42 )  
  [4] Return

  Whole-Sample Symmetric Lifting: Inverse 4,2 Biorthogonal Filter Transform

  wsl42ifilter( u[], N, dq ):
  [0] Compute lnormalize ( u[], N, dq, 1/zeta42 )
  [1] Compute wslpredict ( u[], N, dq, -gamma42 )
  [2] Compute wslupdate ( u[], N, dq, -beta42 )
  [3] Compute wslpredict ( u[], N, dq, -alpha42 )
  [4] Return

*/

/*  Lifting coefficients for 4,2-biorthogonal DWT (exact values) */
const float	alpha42 = { (float)  0.2500 },	/*  1/4  */
		beta42  = { (float)  1.0000 },	/*   1   */
		gamma42 = { (float) -0.1875 },	/* -3/16 */
		zeta42  = { (float)  0.5000 };	/*  1/2  */

void
wsl42filter(float u[], int N, int dq )
{
  wslpredict ( u, N, dq, alpha42 );
  wslupdate  ( u, N, dq, beta42 );
  wslpredict ( u, N, dq, gamma42 );
  lnormalize ( u, N, dq, zeta42 );
  return;
}

void
wsl42ifilter(float u[], int N, int dq )
{
  lnormalize ( u, N, dq, 1/zeta42 );
  wslpredict  ( u, N, dq, -gamma42 );
  wslupdate ( u, N, dq, -beta42 );
  wslpredict  ( u, N, dq, -alpha42 );
  return;
}


/*

  Whole-Sample Symmetric Lifting 4,2-Biorthogonal Discrete Wavelet Transform

  wsl42dwt( u[], N, dq, J ):
  [0] If J>0, then  do [1] and [2]
  [1]    Compute wsl42filter( u[], N, dq )
  [2]    Compute wsl42dwt( u[], N, 2*dq, J-1 )
  [3] Return

  Whole-Sample Symmetric Lifting 4,2-Biorthogonal Inverse Discrete Wavelet Transform

  wsl42idwt( u[], N, dq, J ):
  [0] If J>0, then  do [1] and [2]
  [1]    Compute wsl42idwt( u[], N, 2*dq, J-1 )
  [2]    Compute wsl42ifilter( u[], N, dq )
  [3] Return

*/

void
wsl42dwt( float u[], int N, int dq, int J)
{
  if(J>0) {
    wsl42filter( u, N, dq);
    wsl42dwt( u, N, 2*dq, J-1);
  }
  return;
}

void
wsl42idwt( float u[], int N, int dq, int J)
{
  if(J>0) {
    wsl42idwt( u, N, 2*dq, J-1);
    wsl42ifilter( u, N, dq);
  }
  return;
}

/*

  Whole-Sample Symmetric Lifting: 9,7 Biorthogonal Filter Transform

  wsl97filter( u[], N, dq ):
  [0] Compute wslpredict ( u[], N, dq, alpha97 )
  [1] Compute wslupdate ( u[], N, dq, beta97 )
  [2] Compute wslpredict ( u[], N, dq, gamma97 )
  [3] Compute wslupdate ( u[], N, dq, delta97 )
  [4] Compute lnormalize ( u[], N, dq, zeta97 )  
  [5] Return

  Whole-Sample Symmetric Lifting: Inverse 9,7 Biorthogonal Filter Transform

  wsl97ifilter( u[], N, dq ):
  [0] Compute lnormalize ( u[], N, dq, 1/zeta97 )
  [2] Compute wslupdate ( u[], N, dq, -delta97 )
  [1] Compute wslpredict ( u[], N, dq, -gamma97 )
  [2] Compute wslupdate ( u[], N, dq, -beta97 )
  [3] Compute wslpredict ( u[], N, dq, -alpha97 )
  [4] Return

*/

/*  Lifting coefficients for 9,7-biorthogonal DWT (approximate values) */
const float	alpha97 = { (float) -1.5861343420599235584283154513374 },
		beta97  = { (float) -0.0529801185729614146241295675035 },
		gamma97 = { (float)  0.882911075530933295919790099003  },
		delta97	= { (float)  0.4435068520439711521156042151689 },
		zeta97	= { (float)  1.149604398860241159795075642191  };

void
wsl97filter(float u[], int N, int dq )
{
  wslpredict ( u, N, dq, alpha97 );
  wslupdate  ( u, N, dq, beta97 );
  wslpredict ( u, N, dq, gamma97 );
  wslupdate  ( u, N, dq, delta97 );
  lnormalize ( u, N, dq, zeta97 );
  return;
}

void
wsl97ifilter(float u[], int N, int dq )
{
  lnormalize ( u, N, dq, 1/zeta97 );
  wslupdate  ( u, N, dq, -delta97 );
  wslpredict ( u, N, dq, -gamma97 );
  wslupdate  ( u, N, dq, -beta97 );
  wslpredict ( u, N, dq, -alpha97 );
  return;
}


/*

  Whole-Sample Symmetric Lifting 9,7-Biorthogonal Discrete Wavelet Transform

  wsl97dwt( u[], N, dq, J ):
  [0] If J>0, then  do [1] and [2]
  [1]    Compute wsl97filter( u[], N, dq )
  [2]    Compute wsl97dwt( u[], N, 2*dq, J-1 )
  [3] Return

  Whole-Sample Symmetric Lifting 9,7-Biorthogonal Inverse Discrete Wavelet Transform

  wsl97idwt( u[], N, dq, J ):
  [0] If J>0, then  do [1] and [2]
  [1]    Compute wsl97idwt( u[], N, 2*dq, J-1 )
  [2]    Compute wsl97ifilter( u[], N, dq )
  [3] Return

*/

void
wsl97dwt( float u[], int N, int dq, int J)
{
  if(J>0) {
    wsl97filter( u, N, dq);
    wsl97dwt( u, N, 2*dq, J-1);
  }
  return;
}

void
wsl97idwt( float u[], int N, int dq, int J)
{
  if(J>0) {
    wsl97idwt( u, N, 2*dq, J-1);
    wsl97ifilter( u, N, dq);
  }
  return;
}