Wavelet.c

/*Template:C

• dwt97.c - Fast discrete biorthogonal CDF 9/7 wavelet forward and inverse transform (lifting implementation)
• This code is provided "as is" and is given for educational purposes.
• 2006 - Gregoire Pau - gregoire.pau@ebi.ac.uk - http://www.ebi.ac.uk/~gpau/misc/dwt97.c
• /

1. include <stdio.h>
2. include <stdlib.h>

double *tempbank=0;

/*

• fwt97 - Forward biorthogonal 9/7 wavelet transform (lifting implementation)
• x is an input signal, which will be replaced by its output transform.
• n is the length of the signal, and must be a power of 2.
• The first half part of the output signal contains the approximation coefficients.
• The second half part contains the detail coefficients (aka. the wavelets coefficients).
• See also iwt97.
• /

void fwt97(double* x,int n) { double a; int i;

// Predict 1 a=-1.586134342; for (i=1;i<n-2;i+=2) x[i]+=a*(x[i-1]+x[i+1]); x[n-1]+=2*a*x[n-2];

// Update 1 a=-0.05298011854; for (i=2;i<n;i+=2) x[i]+=a*(x[i-1]+x[i+1]); x[0]+=2*a*x[1];

// Predict 2 a=0.8829110762; for (i=1;i<n-2;i+=2) x[i]+=a*(x[i-1]+x[i+1]); x[n-1]+=2*a*x[n-2];

// Update 2 a=0.4435068522; for (i=2;i<n;i+=2) x[i]+=a*(x[i-1]+x[i+1]); x[0]+=2*a*x[1];

// Scale a=1/1.149604398; for (i=0;i<n;i++) if (i%2) x[i]*=a; else x[i]/=a;

// Pack if (tempbank==0) tempbank=(double *)malloc(n*sizeof(double)); for (i=0;i<n;i++) if (i%2==0) tempbank[i/2]=x[i]; else tempbank[n/2+i/2]=x[i]; for (i=0;i<n;i++) x[i]=tempbank[i]; }

/*

• iwt97 - Inverse biorthogonal 9/7 wavelet transform
• This is the inverse of fwt97 so that iwt97(fwt97(x,n),n)=x for every signal x of length n.
• See also fwt97.
• /

void iwt97(double* x,int n) { double a; int i;

// Unpack if (tempbank==0) tempbank=(double *)malloc(n*sizeof(double)); for (i=0;i<n/2;i++) { tempbank[i*2]=x[i]; tempbank[i*2+1]=x[i+n/2]; } for (i=0;i<n;i++) x[i]=tempbank[i];

// Undo scale a=1.149604398; for (i=0;i<n;i++) if (i%2) x[i]*=a; else x[i]/=a;

// Undo update 2 a=-0.4435068522; for (i=2;i<n;i+=2) x[i]+=a*(x[i-1]+x[i+1]); x[0]+=2*a*x[1];

// Undo predict 2 a=-0.8829110762; for (i=1;i<n-2;i+=2) x[i]+=a*(x[i-1]+x[i+1]); x[n-1]+=2*a*x[n-2];

// Undo update 1 a=0.05298011854; for (i=2;i<n;i+=2) x[i]+=a*(x[i-1]+x[i+1]); x[0]+=2*a*x[1];

// Undo predict 1 a=1.586134342; for (i=1;i<n-2;i+=2) x[i]+=a*(x[i-1]+x[i+1]); x[n-1]+=2*a*x[n-2]; }

int main() { double x[32]; int i;

// Makes a fancy cubic signal for (i=0;i<32;i++) x[i]=5+i+0.4*i*i-0.02*i*i*i;

// Prints original sigal x printf("Original signal:\n"); for (i=0;i<32;i++) printf("x[%d]=%f\n",i,x[i]); printf("\n");

// Do the forward 9/7 transform fwt97(x,32);

// Prints the wavelet coefficients printf("Wavelets coefficients:\n"); for (i=0;i<32;i++) printf("wc[%d]=%f\n",i,x[i]); printf("\n");

// Do the inverse 9/7 transform iwt97(x,32);

// Prints the reconstructed signal printf("Reconstructed signal:\n"); for (i=0;i<32;i++) printf("xx[%d]=%f\n",i,x[i]); }