/* convolution.c  test ffts use for convolution */
/*                16 point, but must antialias using complex */
/* definition
 * f(t) = h(t) o g(t) = integral, T over range of t, of h(T)*g(t-T) dT
 *                                symmetric about 0
 * f(t) = h(t) 0 g(t) = sum, T over range range of t, of h(T)*g(t-(max(t)/2-T))
 *                           t from 0 to N-1 discrete
 *                           max(t)=N-1
 */

#include <stdio.h>
#include <math.h>
#include "fftc.h"

void print16c(char msg[], float a[])
{
  int i;
 
  printf("\n %s \n", msg);
  for(i=0; i<32; i=i+2) printf("%s[%d] =  %g, %g\n",
                                msg, i, a[i], a[i+1]);
}

void print16(char msg[], float a[])
{
  int i;
 
  printf("\n %s \n", msg);
  for(i=0; i<16; i++) printf("%s[%d] =  %g \n",
                               msg, i, a[i]);
}

int main(int argc, char * argv[])
{
  float ant[32]={0.00001, 0.0,  0.00005, 0.0,  0.0003, 0.0,
                 0.002,   0.0,  0.02,    0.0,  0.4,    0.0,
                 0.8,     0.0,  1.0,     0.0,
                 1.0,     0.0,  0.8,     0.0,  0.4,    0.0,
                 0.02,    0.0,  0.002,   0.0,  0.0003, 0.0,
                 0.00005, 0.0,  0.00001, 0.0};
  float tgt[32]={0.0, 0.0,  0.0, 0.0,  0.0, 0.0,  0.0, 0.0,
                 1.0, 0.0,  1.0, 0.0,  1.0, 0.0,  1.0, 0.0,
                 1.0, 0.0,  1.0, 0.0,  0.0, 0.0,  0.0, 0.0,
                 0.0, 0.0,  0.0, 0.0,  0.0, 0.0,  0.0, 0.0};
  float A[32], B[32], C[32];
  float aa[16] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0,
                  1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
  float bb[16] = {0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0,
                  0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0};
  float cc[16]; /* "long" way */
  float Amean, Bmean, Csum, rCsum;
  int i, j;
  int N;

  printf("convolution.c \n");

  N = 32; /* number of floats in 16 point complex fft */
  for(i=0; i<N; i++) {A[i] = ant[i]; B[i] = tgt[i];}
  print16c("A", A);
  fft16(A);
  print16c("fft(A)", A);
  print16c("B", B);
  fft16(B);
  print16c("fft(B)", B);
  /* complex multiply */
  for(i=0; i<N; i=i+2) {C[i]   = A[i]*B[i]-A[i+1]*B[i+1];
                        C[i+1] = A[i]*B[i+1]+A[i+1]*B[i];}
  print16c("C", C);
  ifft16(C);
  print16c("ifft(C)", C);

  /* normalize before fft, make mean zero, complex part still zero */
  for(i=0; i<N; i++) {A[i] = ant[i]; B[i] = tgt[i];}
  Amean = 0.0;
  Bmean = 0.0;
  for(i=0; i<N; i=i+2) {Amean += A[i]; Bmean += B[i];}
  Amean = Amean/16.0;
  Bmean = Bmean/16.0;
  for(i=0; i<N; i=i+2) {A[i]= A[i]-Amean; B[i] = B[i]-Bmean;}
  print16c("rmeanA", A);
  fft16(A);
  print16c("fft(rmeanA)", A);
  print16c("rmeanB", B);
  fft16(B);
  print16c("fft(rmeanB)", B);
  /* complex multiply */
  for(i=0; i<N; i=i+2) {C[i]   = A[i]*B[i]-A[i+1]*B[i+1];
                        C[i+1] = A[i]*B[i+1]+A[i+1]*B[i];}
  print16c("rmeanC", C);
  ifft16(C);
  print16c("fft(rmeanC)", C);
  Csum = 0.0;
  for(i=0; i<N; i=i+2) Csum = Csum + C[i]; /* just real */
  printf("rCsum=%g\n",rCsum/16.0);

  for(i=0; i<16; i++) {A[2*i] = aa[i]; B[2*i] = bb[i];
                        A[2*i+1] = 0.0; B[2*i+1] = 0.0;}
  print16c("aa", A);
  fft16(A);
  print16c("fft(aa)", A);
  print16c("bb", B);
  fft16(B);
  print16c("fft(bb)", B);
  /* complex multiply */
  for(i=0; i<N; i=i+2) {C[i]   = A[i]*B[i]-A[i+1]*B[i+1];
                        C[i+1] = A[i]*B[i+1]+A[i+1]*B[i];}
  print16c("cc", C);
  ifft16(C);
  print16c("ifft(cc)", C);

  /* numerical convolution */
  print16("aa", aa);
  print16("bb", bb);
  for(i=0; i<16; i++)
  {
    cc[i] = 0.0;
    for(j=0; j<16; j++) /* must center (7-j) is Tau */
    {
      if(i-(7-j)<0 || i-(7-j)>=16) continue;
      cc[i] = cc[i] + aa[j]*bb[i-(7-j)];
    }
  }
  print16("cc", cc);
  return 0;
} /* end convolution.c */