Algorithm Alley 
by Steven Pigeon


Listing One
#include <imath.h>
#include <stdlib.h>
#include <math.h>
/***************************************************************************/
/* Those two routines implement the Mallat decimation algorithm. First     */
/* phase computes the means, second phase the differences. All isO(2n),    */
/* which is rather efficient.                                              */
/* keeps ~20 bits of precision (out of 23), or a precision of about 1e-7   */
/***************************************************************************/
void FMT(const float source[], float destination[], int length)
 { 
  float *t;
  int i,d,j;
  t = (float *)malloc(2*length*sizeof(float));
  for (i=0;i<length;i++) t[i]=source[i];
  for (i=0,d=length;i<2*length-2;i+=2,d++) t[d] = t[i]+t[i+1];
  destination[length-1]=t[2*length-2];
  for ( i=2*length-2, j=2*length-3, d=length-2; 
        i>length-1; 
        i--, j-=2, d--) 
    destination[d] = t[j]-t[i]/2.0f;
  free(t);
 }
/*******************************/
void iFMT(const float source[], float destination[], int length)
 {
  float *t; 
  int i,d,p;
  t = (float *)malloc(2*length*sizeof(float));
  t[2*length-2]=source[length-1];
  for ( p=length-2, i=2*length-4 , d=2*length-2; 
       i>=0; i-=2, d--, p--)
   {
    float t0 = t[d]/2;
    float s0 = source[p];
    t[i]   = t0 - s0;
    t[i+1] = t0 + s0;
   }
  for (i=0;i<length;i++) destination[i]=t[i];
  free(t);
 }
/***************************************************************************/
/* These routines compute the fast (normalized) Haar transform             */
/* Works of course with a vector of any lenght > 3.                        */
/* Lonely coefficients are left as is, since we can't pair them correctly, */
/* but are put in the result vector. The precision is also about ~20 bits  */
/* (out of 23)                                                             */
/***************************************************************************/

void FHT(const float source[], float destination[], int length)
 {
  int i,l;
  float *t = (float *)malloc(length*sizeof(float));
  for (i=0;i<length;i++) destination[i]=source[i];
  l=length;
  for (i=Log2(length-1);i;i--)
   {
    int j,k;
    int l2 = l/2;
    for(j=0;j<l;j++) t[j]=destination[j];
    for (k=0,j=0;k<l2;k++,j+=2)
     {
      destination[k]    = (t[j]+t[j+1])*0.7071067812f;
      destination[k+l2] = (t[j]-t[j+1])*0.7071067812f;
     }
    l /= 2 ;  // for the leftovers
   }
  free(t);
 }
/*******************************/
void iFHT(const float source[], float destination[], int length)
 {
  int i,j;
  int lim[32];
  
  int    k  = Log2(length-1);
  float *t  = (float *)malloc( length * sizeof(float));

  if (( 1 << k) != length) k--;

  for ( j=length, i=k; i; i--, j/=2 ) lim[i]=j;
  lim[0]=1;

  for ( i=0; i<length; i++ ) destination[i]=source[i];

  for(i=0;i<k;i++)
   {
    int d;
    int l=lim[i];
    int l2=lim[i+1];

    for(j=0; j<l2; j++) t[j]=destination[j];

    for( d=0, j=0; j<l; j++, d+=2 )
     {
      destination[d]   =(t[j]+t[j+l])*0.7071067812f;
      destination[d+1] =(t[j]-t[j+l])*0.7071067812f;
     }
   }
  free(t);
 }






3