/*
   Lifting Implementation of the Haar Filter Transform on $2q$ Samples

   haarlift( u[], q, dq ):
   [1] For n=0 to q-1, replace u[(2*n+1)*dq] -=  u[(2*n)*dq]
   [2] For n=0 to q-1, replace u[(2*n)*dq] += 0.5 * u[(2*n+1)*dq]
   [3] For n=0 to q-1, replace u[(2*n+1)*dq] /= sqrt(2.0)
   [4] For n=0 to q-1, replace u[(2*n)*dq] *= sqrt(2.0)


   Invert the Haar Filter Transform to $2q$ samples

   ihaarlift( u[], q, dq ):
   [1] For n=0 to q-1, replace u[(2*n)*dq] /= sqrt(2.0)
   [2] For n=0 to q-1, replace u[(2*n+1)*dq] *= sqrt(2.0)
   [3] For n=0 to q-1, replace u[(2*n)*dq] -= 0.5 * u[(2*n+1)*dq]
   [4] For n=0 to q-1, replace u[(2*n+1)*dq] +=  u[(2*n)*dq]



   Complete Discrete Haar Transform by Lifting on $2^J$ Samples

   ldht0( u[], J, dq ):
   [0] If J>0, then  do [1] and [2]
   [1]    Compute haarlift( u[], (1<<J)/2, dq ):
   [2]    Compute ldht0( u[], J-1, 2*dq )
   [3] Return

   Reconstruction by Lifting from a Complete Discrete Haar Transform
   on $2^J$ Samples 

   ildht0( u[], J, dq ):
   [0] If J>0, then  do [1] and [2]
   [1]    Compute ildht0( u[], J-1, 2*dq )
   [2]    Compute ihaarlift( u[], (1<<J)/2, dq ):
   [4] Return

*/

#include <math.h>
#include <stdio.h>

void
haarlift( float u[], int q, int dq)
{
  int n;
  for(n=0; n<q; n++) u[(2*n+1)*dq] -=  u[(2*n)*dq];
  for(n=0; n<q; n++) u[(2*n)*dq] += 0.5 * u[(2*n+1)*dq];
  for(n=0; n<q; n++) u[(2*n+1)*dq] /= sqrt(2.0);
  for(n=0; n<q; n++) u[(2*n)*dq] *= sqrt(2.0);
}

void
ihaarlift( float u[], int q, int dq)
{
  int n;
  for(n=0; n<q; n++) u[(2*n)*dq] /= sqrt(2.0);
  for(n=0; n<q; n++) u[(2*n+1)*dq] *= sqrt(2.0);
  for(n=0; n<q; n++) u[(2*n)*dq] -= 0.5 * u[(2*n+1)*dq];
  for(n=0; n<q; n++) u[(2*n+1)*dq] +=  u[(2*n)*dq];
}

void
ldht0(float u[], int J, int dq )
{
  if(J>0) {
    haarlift( u, (1<<J)/2, dq );
    ldht0( u, J-1, 2*dq );
  }
  return;
}

void
ildht0(float u[], int J, int dq )
{
 if(J>0) {
   ildht0( u, J-1, 2*dq );
   ihaarlift( u, (1<<J)/2, dq );
 }
 return;
}

int
main(void)
{
  float in[8], in1[8]={0}, save[8]={0};
  int n, N=8, J=3;

  for(n=0; n<N; n++)
    in1[n]=save[n]=in[n]=(n-N/2.0)*(n-N/2.0)-N*N/8.0;	/* test signal */


  printf("   u[%d:%d]:", 0,N-1);
  for(n=0; n<N; n++) printf(" %f", in[n]);
  putchar('\n');

  haarlift(in, N/2, 1);

  printf("  Fu[%d:%d]:", 0,N-1);
  for(n=0; n<N; n++) printf(" %f", in[n]);
  putchar('\n');

  ihaarlift(in, N/2, 1);

  printf("F*Fu[%d:%d]:", 0,N-1);
  for(n=0; n<N; n++) printf(" %f", in[n]);
  putchar('\n');

  printf("diff[%d:%d]:", 0,N-1);
  for(n=0; n<N; n++) printf(" %f", in[n]-save[n]);
  putchar('\n');

  ldht0( in1, J, 1 );

  printf("  Hu[%d:%d]:", 0,N-1);
  for(n=0; n<N; n++) printf(" %f", in1[n]);
  putchar('\n');

  ildht0( in1, J, 1 );

  printf("H*Hu[%d:%d]:", 0,N-1);
  for(n=0; n<N; n++) printf(" %f", in1[n]);
  putchar('\n');

  printf("diff[%d:%d]:", 0,N-1);
  for(n=0; n<N; n++) printf(" %f", in1[n]-save[n]);
  putchar('\n');

  return 0;
}