/* This is a sample program of 1D Walsh-Hadamard transform (WHT).
   Compile it with "cc mywht.c -lm -o mywht". 
   You can choose either Hadamard order or Sequency (Walsh) order by 
   setting the switch "sequency" in the function wht(x,m,sequency,inverse).
   
   You can use the wht function in this program for your project. I wrote
   this wht function following closely what I showed you in class (and also
   the handout). Do not treat it as a black box. Read it to understand it.
*/
   
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>

#include "array_utility.h"

int howo(i,m)    /*  converts a sequency index i to Haramard index j  */
int i,m;
{ int j,k;
/*  j=i>>1;  i=i^j;  /* converting binary to Gray code */
  i=i^(i>>1);
  j=0;
  for (k=0; k<m; ++k)
    j=(j << 1) | (1 & (i >> k));    /* bit-reversal */
  return j;
}

void wht(x,N,sequency,inverse)
     float *x;
     int N,inverse,sequency;
{ int mmax,step,i,j,k,j1,m,n,nu;
  float w;

  m=log((float)N)/log(2.0);
  if (sequency) {
    for (i=0; i<N; ++i) {       // bit reversal of data
      j=0;
      for (k=0; k<m; ++k)
	j=(j << 1) | (1 & (i >> k));
      if (j < i) {
	w=x[i]; x[i]=x[j]; x[j]=w;
      }
    }
  }
  for (i=0; i<m; i++) {         // for log2(N) stages 
    n=pow(2,m-1-i);             // length of section 
    k=0;
    while (k<N-1) {             // before reaching end
      for (j=0; j<n; j++) {     // for each section 
	j1=k+j;
	w=x[j1]+x[j1+n];
	x[j1+n]=x[j1]-x[j1+n];
	x[j1]=w;
      }
      k+=2*n;
    }
  }
  w=1.0/sqrt((float)N);
  for (i=0; i<N; i++) x[i]=x[i]*w;

  /*
  for (i=0; i<N; i++)
    printf("%4.1f ",x[i]);
  printf("\n");
  
  if (sequency) {     // converting to sequency (Walsh) order 
    for (i=0; i<N; i++) { 
      j=howo(i,m);  
      printf("(%d %d) ",i,j);
      y[i]=x[j];
      //      if (j<i) { w=x[i]; x[i]=x[j]; x[j]=w; }
    }
    for (i=0; i<N; i++) x[i]=y[i];
  }
  for (i=0; i<N; i++)
    printf("%4.1f ",x[i]);
  printf("\n");
  */
}

main()
{
  int i,j,k, m, N;
  float x[8];
  printf("Ordering: (0) Hadamard, (1) Sequency (Walsh): ");
  scanf("%d",&k);

  N=8;
  m=log((float)N)/log(2.0);
  //  n=pow(2,m);
  for (i=0; i<N; i++) { x[i]=0;}
  for (i=2; i<5; i++) { x[i]=i;}

  printf("\nx in time domain:\n");
  for (i=0; i<N; i++) {
    printf("%4.1f ",x[i]);
  }
  printf("\n");
  
  wht(x,N,k,0);         /* forward WHT */

  printf("X in sequency domain:\n");
  for (i=0; i<N; i++) {
    printf("%4.1f ",x[i]);
  }
  printf("\n");

/* This is where you should insert your filtering codes. */

  wht(x,N,k,1);         /* inverse WHT */

  printf("x back to time domain:\n");
  for (i=0; i<N; i++) {
    printf("%4.1f ",x[i]);
  }
  printf("\n");
}