/*
  Mallat's Discrete Wavelet Transform

  dwt( u[], x, y, J, h[], g[], L):
  [0]  If J=0 then for n = x to y, print u[n]
  [1]  Else do [2] through [9]
  [2]     Let x1 = ceiling((1+x-L)/2), let y1 = floor(y/2)
  [3]     For n = x1 to y1 do [4] through [8]
  [4]        Let s[n] = 0, let d[n] = 0 
  [5]        For k = 0 to L-1 do [6] and [7]
  [6]           Accumulate s[n] += h[k]*u[k+2*n] 
  [7]           Accumulate d[n] += g[k]*u[k+2*n] 
  [8]        Print d[n]
  [9]     Compute dwt( s[], x1, y1, J-1, h[], g[], L) 
  [10] Return

  */

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

void
dwt( float u[], int x, int y, int J, const float h[], const float g[], int L)
{
  int k, n, x1, y1;
  float *s, *d, *s0, *d0;
  if(J==0) {
    printf("s_{J-0}[%d:%d]", x,y);
    for(n=x; n<=y; n++) printf(" %f", u[n]);
    putchar('\n');
  }
  else {
    x1 = (int)ceil((float)(1+x-L)/2.0);
    y1 = (int)floor(y/2.0);
    s0 = (float *)malloc((1+y1-x1)*sizeof(float)); assert(s0); s=s0-x1;
    d0 = (float *)malloc((1+y1-x1)*sizeof(float)); assert(d0); d=d0-x1;
    printf("d_{J-%d}[%d:%d]:", J-1, x1, y1);
    for(n=x1; n<=y1; n++ ) {
      s[n] = d[n] = 0;
      for(k=0; k<L; k++) {
	s[n] += h[k]*u[k+2*n];
	d[n] += g[k]*u[k+2*n];
      }
      printf(" %f", d[n]);
    }
    putchar('\n');
    dwt( s, x1, y1, J-1, h, g, L);
    free(s0); free(d0);
  }
  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[17];
  int n, x=0, y=17-1;

  for(n=x; n<=y; n++) u[n]=(float)(n*n-16*n+41);	/* test signal */
  dwt(u,x,y,J,h,g,L);
  return 0;
}