Fast discrete biorthogonal CDF 9/7 wavelet forward and inverse transform (lifting implementation)

dwt97.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
 */

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

#define FSHIFT  8
#define FMUL    ((double)(1 << FSHIFT))
#define FMUL_I  ((int)(1 << FSHIFT))



#define print_matrix_type(t,fmt_t,sz) \
static void print_matrix_##sz##t( t *matrix) { \
   int n = sqrt(sz); \
   for (int i=0; i < sz; i++) { \
      printf(fmt_t, matrix[i]); \
         if (i && i%n == (n-1) || i==(sz-1)) \
            printf("\n"); \
         else \
            printf(" "); \
   } \
}

print_matrix_type(double,"%1.2f",64)
print_matrix_type(int,"%d",64)

#define print_matrix_8x8_double(a) print_matrix_64double(a)
#define print_matrix_8x8t_double(t,a) printf("%s\n",t); print_matrix_8x8_double(a); printf("\n");


#define print_matrix_8x8_int(a) print_matrix_64int(a)
#define print_matrix_8x8t_int(t,a) printf("%s\n",t); print_matrix_8x8_int(a); printf("\n");


void fwt97_int(int* x,int n) {
  int a;
  int i;
  alignas(32) int tempbank[n];

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

  print_matrix_8x8t_int("i after predict 1", x);

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

  print_matrix_8x8t_int("i after update 1", x);

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

  print_matrix_8x8t_int("i after predict 2", x);

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

  print_matrix_8x8t_int("i after Update 2", x);

  // Scale
  const int scale_dn = round(1/1.149604398*FMUL);
  const int scale_up = round(1.149604398*FMUL);
  for (i=0;i<n;i++) {
    x[i] = x[i] * ((i%2) ? scale_dn : scale_up) / FMUL_I;
  }
  // Pack
  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];

  print_matrix_8x8t_int("i after pack", x);
}

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

  // Unpack
  alignas(32) int tempbank[n];

  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
  for (i=0;i<n;i++) {
    x[i] = x[i] * (int)round( (i%2) ? 1.149604398*FMUL : 1./1.149604398*FMUL ) / FMUL_I;
  }
  // Undo update 2
  a=round(-0.4435068522*FMUL);
  for (i=2;i<n;i+=2) {
    x[i]+=a*(x[i-1]+x[i+1]) / FMUL_I;
  }
  x[0]+=2*a*x[1] / FMUL_I;

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

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

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

}



/**
 *  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;
  alignas(32) double tempbank[n];

  // 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];

  print_matrix_8x8t_double("i after Predict 1", x);

  // 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];
  print_matrix_8x8t_double("i after Update 1", x);

  // 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];
  print_matrix_8x8t_double("i after Predict 2", x);

  // 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];
  print_matrix_8x8t_double("i after Update 2", x);

  // Scale
  a=1/1.149604398;
  for (i=0;i<n;i++) {
    x[i] = x[i] * ( (i%2) ? 1/1.149604398 : 1.149604398 );
  }

  // Pack
  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];
  print_matrix_8x8t_double("i after pack", x);
}

/**
 *  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
  alignas(32) double tempbank[n];

  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
  for (i=0;i<n;i++) {
    x[i] = x[i] * ( (i%2) ? 1.149604398 : 1./1.149604398 );
  }

  // 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];
}


#define fwt97_double fwt97
#define iwt97_double iwt97

#define cdw97_type(t) \
static int cdw97_##t() {         \
  t x[64];            \
  int i;              \
  /* Makes a fancy cubic signal */ \
  for (i=0;i<64;i++) x[i]=i-32;    \
                                   \
  /* Prints original sigal x */    \
  print_matrix_8x8t_##t("Original signal:", x); \
                                    \
  /* Do the forward 9/7 transform */ \
  fwt97_##t(x,64); \
                  \
  /* Prints the wavelet coefficients */ \
  print_matrix_8x8t_##t("Wavelets coefficients:",x); \
                                        \
  /* Do the inverse 9/7 transform */ \
  iwt97_##t(x,64); \
\
  /* Prints the reconstructed signal */ \
  print_matrix_8x8t_##t("Reconstructed signal:",x); \
} \

cdw97_type(double)
cdw97_type(int)

int main() {
   cdw97_double();
   printf("-- integer --\n");
   cdw97_int();
}