/*

  Complete Periodic Wavelet Transform on $2^J$ Samples

  pdwt0( u[], J, h[], g[], L):
  [0] If J>0, then  do [1-] to [4]
  [1-]   Let N = 1<<J
  [1]    Allocate temp[0]=0,...,temp[N-1]=0
  [2]    Compute pcqfilter( temp[], u[], N/2, h[], g[], L):
  [3]    For i=0 to N-1, let u[i] = temp[i]
  [4]    Compute pdwt0( u[], J-1, h[], g[], L )
  [5] Return


  Reconstruction of $2^J$ Samples from Complete Periodic Wavelet Expansion

  ipdwt0( u[], J, h[], g[], L):
  [0] If J>0, then  do [1-] through [4]
  [1-]   Let N = 1<<J
  [1]    Allocate temp[0]=0,...,temp[N-1]=0
  [2]    Compute ipdwt0( u[], J-1, h[], g[], L )
  [3]    Compute ipcqfilter( temp[], u[], N/2, h[], g[], L )
  [4]    For i=0 to N-1, let u[i] = temp[i]
  [5] Return

*/

#include <assert.h>
#include <stdio.h>
#include <stdlib.h>

#include "pcqfilt.c"

void
pdwt0(float u[], int J, const float h[], const float g[], int L)
{
  if(J>0) {
      float *temp;
      int i, N;
      N = 1<<J;
      temp = (float *)calloc(N, sizeof(float)); assert(temp);
      pcqfilter(temp, u, N/2, h,g,L);
      for(i=0;i<N;i++) u[i]=temp[i];
      free(temp);
      pdwt0(u, J-1, h,g,L);
  }
  return;
}

void
ipdwt0(float u[], int J, const float h[], const float g[], int L)
{
  if(J>0) {
      float *temp;
      int i, N;
      N = 1<<J;
      ipdwt0(u, J-1, h,g,L);
      temp = (float *)calloc(N, sizeof(float)); assert(temp);
      ipcqfilter(temp, u, N/2, h,g,L);
      for(i=0;i<N;i++) u[i]=temp[i];
      free(temp);
  }
  return;
}


int
main(void)
{
  /* Daubechies 4 filter coefficients: */
  const float h[4] = {  0.48296291314453416, 0.83651630373780794,
			0.22414386804201339, -0.12940952255126037};
  const float g[4] = { -0.12940952255126037, -0.22414386804201339,
		       0.83651630373780794, -0.48296291314453416};
  const int L=4, J=4;
  float u[16];
  int n;

  for(n=0; n<16; n++) u[n]=(float)(n*n-16*n+41);	/* test signal */

  printf("   u[%d:%d]:", 0,15);
  for(n=0; n<16; n++) printf(" %10g", u[n]);
  putchar('\n');

  pdwt0(u,J,h,g,L);

  printf("  Wu[%d:%d]:", 0,15);
  for(n=0; n<16; n++) printf(" %f", u[n]);
  putchar('\n');

  ipdwt0(u,J,h,g,L);

  printf("iWWu[%d:%d]:", 0,15);
  for(n=0; n<16; n++) printf(" %10g", u[n]);
  putchar('\n');

  return 0;
}