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

#define PI 3.141592653589
#define SWAP(a,b) tempr=(a); (a)=(b); (b)=tempr

void pause()
{
  int k;
  printf("\npaused, enter 1 to continue...\n");
  scanf("%d",&k);
}

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 bitreversal(x,N)   /*  bit reverses a given vector x[] */
float x[];
int N;
{ int i,j,k, m;
  float tempr;
  m=(int)log2f((float)N);
  for (i=0; i<N; ++i) {       /* bit reversal */
    j=0;
    for (k=0; k<m; ++k)
      j=(j << 1) | (1 & (i >> k));
    if (j < i) {
      SWAP(x[i],x[j]);
    }
  }
}

fft(xr,xi,N,inverse)
     float xr[],xi[];
     int N,inverse;
{ int mmax,step,i,j,k,j1,n,nu,m;
  float arg,s,c,tempr,tempi,v;

  //  N=pow(2.0,(float)m);
  m=(int)log2f((float)N);
  /*
  if (inverse) printf("IFFT of size %d (%d)\n",N,m);
  else printf("FFT of size %d (%d)\n",N,m);
  printf("\n before reversal...\n");
  for (i=0; i<N; i++)  printf("%3.0f ",xr[i]); printf("\n");
  */
  bitreversal(xr,N); bitreversal(xi,N);   /* bit reversal of input data */
  /*
  printf("\n after reversal...\n");
  for (i=0; i<N; i++)  printf("%3.0f ",xr[i]); printf("\n");
  */
  for (i=0; i<m; i++) {         /* for m=log N stages */
    n=pow(2.0,(float)i);                 /* length of a section */
    k=0;
    while (k<N-1) {             /* for N components in each stage */
      for (j=0; j<n; j++) {     /* for each section */
	arg=-j*PI/n;
	if (inverse) arg=-arg;
	c=cos(arg);
	s=sin(arg);
	j1=k+j;
	tempr=xr[j1+n]*c-xi[j1+n]*s;
	tempi=xi[j1+n]*c+xr[j1+n]*s;
	xr[j1+n]=xr[j1]-tempr;
	xi[j1+n]=xi[j1]-tempi;
	xr[j1]=xr[j1]+tempr;
	xi[j1]=xi[j1]+tempi;
      }
      k+=2*n;
    }
  }
  if (inverse) {
    v=1;
    for (i=0; i<N; i++) 
      { xr[i]=xr[i]/v; xi[i]=xi[i]/v; }
  }
}

fht(x,N,inverse)
     float *x;
     int N,inverse;
{ int i,j,n,m,n1;
  float *tmp,u,v,s,s1;

  tmp=(float *)malloc(N*sizeof(float));
  m=(int)log2f((float)N);
  u=1.0/sqrt((float)N);
  bitreversal(x,N);
  n=N; s=pow(sqrt(2.0),n/2);
  for (i=m-1; i >= 0; i--) {     // for m=log N stages 
    n=n/2; s=s/sqrt(2.0);
    printf("n=%d %d, s=%.3f %.3f\n",n,n1,s,s1);
    for (j=0; j<n; j++) {    
      tmp[j]=x[j]-x[j+n];     // second half 
      x[j]=x[j]+x[j+n];       // first half
      //      printf("x[%d]=%.3f ",j,x[j]);
    }
    bitreversal(tmp,n);
    //    for (j=0; j<n; j++) printf("tmp[%d]=%.3f ",j,tmp[j]);
    //    printf("\n");
    for (j=0; j<n; j++) {
      //      x[n+j]=tmp[j]*s/u;
      //      x[n+j]=tmp[j]*s/u;
      x[n+j]=tmp[j]*s;
      //      printf("x[%d]=%.3f ",n+j,x[n+j]);
    }
    //    printf("\n");
  }
  for (i=0; i<N; i++) x[i]*=u;
}


dht(x,N,inverse)
     float *x;
     int N,inverse;
{
  int i,n;
  float *y;
  
  y=(float *)malloc(N*sizeof(float));
  if (inverse) {
    n=1;
    while(n<N) {
      for (i=0; i<n; i++) {
	y[2*i]  =(x[i]+x[i+n])/sqrt(2.0);
	y[2*i+1]=(x[i]-x[i+n])/sqrt(2.0);
      }
      for (i=0; i<n*2; i++) x[i]=y[i];
      n=n*2;
    }
  }
  else {
    n=N;
    while(n>1) {
      n=n/2;
      for (i=0; i<n; i++) {
	y[i]  =(x[2*i]+x[2*i+1])/sqrt(2.0);
	y[i+n]=(x[2*i]-x[2*i+1])/sqrt(2.0);
      }
      for (i=0; i<n*2; i++) x[i]=y[i];
    }
  }
  free(y);
}


fdct(x,N,inverse)
     float *x;
     int N,inverse;
{
  int m,n,N2=N/2;
  float a,u,v,w, *yr,*yi;

  w=PI/2/N; 
  a=sqrt(2.0/N);
  yr=(float *)malloc(N*sizeof(float));
  yi=(float *)malloc(N*sizeof(float));

  if (inverse) {
    for (n=0; n<N; n++)
      x[n]=x[n]*a;
    x[0]=x[0]/sqrt(2.0);
    for (n=0; n<N; n++) {
      yr[n]=x[n]*cos(n*w);
      yi[n]=x[n]*sin(n*w);
    }
  }
  else {
    for (m=0; m<N2; m++) {
      yr[m]=x[2*m];
      yr[N-1-m]=x[2*m+1];
      yi[m]=yi[N2+m]=0;
    } 
  }

  fft(yr,yi,N,inverse);

  if (inverse) {
    for (m=0; m<N2; m++) {
      x[2*m]=yr[m];
      x[2*m+1]=yr[N-1-m];
    }
  }
  else {
    for (n=0; n<N; n++) 
      x[n]=cos(n*w)*yr[n]+sin(n*w)*yi[n];
    for (n=0; n<N; n++)
      x[n]=x[n]*a;
    x[0]=x[0]/sqrt(2.0);
  }
  free(yr); free(yi);
}


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

  m=log((float)N)/log(2.0);
  y=(float *)malloc(N*sizeof(float));

  for (i=0; i<m; i++) {         // for log 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;
	t=x[j1]+x[j1+n];
	x[j1+n]=x[j1]-x[j1+n];
	x[j1]=t;
      }
      k+=2*n;
    }
  }
  w=1.0/sqrt((float)N);
  for (i=0; i<N; i++) x[i]=x[i]*w;
  if (sequency)     // converting to sequency (Walsh) order 
    {  
      for (i=0; i<N; i++) { j=howo(i,m);  y[i]=x[j]; }
      for (i=0; i<N; i++) x[i]=y[i];
    }
  free(y);
}

dct(x,N,inverse)
float *x;
int N,inverse;
{ int mmax,i,j;
  float *y,u,u0,w,t;

  w=0.5*PI/N;  
  u0=sqrtf(1.0/N); u=sqrtf(2.0/N); 
  y=(float *)malloc(N*sizeof(float));
  if (inverse)
    for (i=0; i<N; i++) {
      y[i]=u0*x[0];
      for (j=1; j<N; j++)
	y[i]+=u*x[j]*cosf((2*i+1)*j*w);
    }
  else
    for (i=0; i<N; i++) { 
      y[i]=0;
      for (j=0; j<N; j++)
	y[i]+=x[j]*cosf((2*j+1)*i*w);
      if (i) y[i]*=u; else y[0]*=u0;
    }
  for (i=0; i<N; i++) x[i]=y[i];  
  free(y);
}


main()
{
  int i,j,k, m, n, N;
  float *x, *y, *z;
  float u,v,w, av,au,pw,qw;

 label:
  printf("N=");  scanf("%d",&N);
  //  printf("DCT (0) or fast DCT (1) "); scanf("%d",&k);

  x = (float *) malloc(N*sizeof(float));
  y = (float *) malloc(N*sizeof(float));

  for (m=2; m<5; m++) x[m]=m;

  printf("Time signal before transform:\n");
  for (m=0; m<N; m++) printf("%6.2f ",x[m]); printf("\n");

  //  fht(x,N);
  dht(x,N,0);

  printf("\nSpectrum after transform:\n");
  for (n=0; n<N; n++) printf("%6.2f ",x[n]);
  printf("\n");

  dht(x,N,1);

  printf("\nTime signal after inverse transform:\n");
  for (n=0; n<N; n++) printf("%6.2f ",x[n]);
  printf("\n");

  pause();

  qw=0; 
  for (n=0; n<N; n++) { qw+=x[n]*x[n]; }
  au=x[0]*sqrtf(N);
  qw=sqrt(qw);
  printf("Frequency domain: Power=%.2f, Average=%.2f\n",qw,au);

  //  if (k) fdct(x,N,1);  else dct(x,N,1);

  wht(x,N,1,1);

  printf("\nTime signal back from DCT:\n");
  for (m=0; m<N; m++) printf("%6.2f ",x[m]); printf("\n");
  v=0;
  for (m=0; m<N; m++) 
    v=v+(x[m]-y[m])*(x[m]-y[m]);
  printf("Should be zero: %.2f\n",v);
  
  free(x);

  if (1) goto label;

}