/* 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 #include #include #include 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