/* test_fft.c  test modified Howard's fftd */
#include <stdio.h>
#include <math.h>
#include "fftd.h"

#define Pi 3.1415926535897932384626433832795028841971
#define n 16

int main(int argc, char *argv[])
{
  double rdat[4*n];  /* real part */
  double idat[4*n];  /* imaginary part */
  double rslt[4*n];  /* cos,sin  pair terms */
  double arslt[4*n]; /* for antialias rslt */
  double AT[4*n], BT[4*n], C[4*n]; /* for convolution */
  double th, v;
  int i;
  int k = 1; /* fft, ifft */

  printf("test_fft.c test fftd.c \n");
  printf("check implies other coefficients should be near zero \n");
  for(i=0; i<n; i++) /* sin 2 th other coefficients zero */
  {
    th = (2.0*2.0*Pi)*(double)i/(double)n;
    rdat[i] = sin(th);
    idat[i] = 0.0;
    rslt[2*i] = 0.0;
    rslt[2*i+1] = 0.0;
  }
  printf("sin 2 th data \n");
  for(i=0; i<n; i++) printf("[%d] real = %g, imag = %g \n",
                            i, rdat[i], idat[i]) ;
  fftd(rslt, rdat, idat, n, 1);
  printf("sin 2 th FFT \n");
  for(i=0; i<n; i++) printf("cos %d th = %g, sin %d th = %g \n",
                            i, rslt[i*2], i, rslt[i*2+1]) ;
  printf("antialias and normalize sin 2 th FFT \n");
  for(i=1; i<n/2; i++)
  {
    rslt[i*2]  = (rslt[i*2]+rslt[(n-i)*2])/(double)n;     /* real part */
    rslt[i*2+1]=-(rslt[i*2+1]-rslt[(n-i)*2+1])/(double)n; /* conj imag part */
  }
  rslt[0]=rslt[0]/(double)n;
  rslt[1]=-rslt[1]/(double)n;
  for(i=0; i<=(n/2)-1; i++) printf("cos %d th = %g, sin %d th = %g \n",
                               i, rslt[i*2], i, rslt[i*2+1]) ;
  printf("check: sin 2 th = 1.0 \n");
  printf("\n");
  printf("\n");


  for(i=0; i<n; i++) /* cos 2 th  other coefficients zero */
  {
    th = (2.0*2.0*Pi)*(double)i/(double)n;
    rdat[i] = cos(th);
    idat[i] = 0.0;
    rslt[2*i] = 0.0;
    rslt[2*i+1] = 0.0;
  }
  printf("cos 2 th data \n");
  for(i=0; i<n; i++) printf("[%d] real = %g, imag = %g \n",
                            i, rdat[i], idat[i]) ;
  fftd(rslt, rdat, idat, n, 1);
  printf("cos 2 th FFT \n");
  for(i=0; i<n; i++) printf("cos %d th = %g, sin %d th = %g \n",
                            i, rslt[i*2], i, rslt[i*2+1]) ;
  printf("antialias and normalize cos 2 th FFT \n");
  for(i=1; i<n/2; i++)
  {
    rslt[i*2]  = (rslt[i*2]+rslt[(n-i)*2])/(double)n;     /* real part */
    rslt[i*2+1]=-(rslt[i*2+1]-rslt[(n-i)*2+1])/(double)n; /* conj imag part */
  }
  rslt[0]=rslt[0]/(double)n;
  rslt[1]=-rslt[1]/(double)n;
  for(i=0; i<=(n/2)-1; i++) printf("cos %d th = %g, sin %d th = %g \n",
                                   i, rslt[i*2], i, rslt[i*2+1]) ;
  printf("check: cos 2 th = 1.0 \n");
  printf("\n");
  printf("\n");


  for(i=0; i<n; i++) /* 1.0 + cos 2 th  other coefficients zero */
  {
    th = (2.0*2.0*Pi)*(double)i/(double)n;
    rdat[i] = 1.0 + cos(th);
    idat[i] = 0.0;
    rslt[2*i] = 0.0;
    rslt[2*i+1] = 0.0;
  }
  printf("1.0 + cos 2 th data \n");
  for(i=0; i<n; i++) printf("[%d] real = %g, imag = %g \n",
                            i, rdat[i], idat[i]) ;
  fftd(rslt, rdat, idat, n, 1);
  printf("cos 2 th FFT \n");
  for(i=0; i<n; i++) printf("cos %d th = %g, sin %d th = %g \n",
                            i, rslt[i*2], i, rslt[i*2+1]) ;
  printf("antialias and normalize 1.0 + cos 2 th FFT \n");
  for(i=1; i<n/2; i++)
  {
    rslt[i*2]  = (rslt[i*2]+rslt[(n-i)*2])/(double)n;     /* real part */
    rslt[i*2+1]=-(rslt[i*2+1]-rslt[(n-i)*2+1])/(double)n; /* conj imag part */
  }
  rslt[0]=rslt[0]/(double)n;
  rslt[1]=-rslt[1]/(double)n;
  for(i=0; i<=(n/2)-1; i++) printf("cos %d th = %g, sin %d th = %g \n",
                                   i, rslt[i*2], i, rslt[i*2+1]) ;
  printf("check: cos 0 th = 1.0  and  cos 2 th = 1.0 \n");
  printf("\n");
  printf("\n");


  for(i=0; i<n; i++) /* square wave 4/Pi(cos th - 1/3 cos 3 th + */
                     /*  1/5 cos 5 th - 1/7 cos 7 th */
  {
    v = 1.0;
    if(i>=n/4 && i<3*n/4) v = -1.0;
    if(i==n/4 || i==3*n/4) v = 0.0;
    rdat[i] = v;
    idat[i] = 0.0;
  }
  printf("square wave data, cos series \n");
  for(i=0; i<n; i++) printf("[%d] real = %g, imag = %g \n",
                            i, rdat[i], idat[i]) ;
  printf("square wave generate, cos series \n");
  for(i=0; i<n; i++)
  {
    th = 2.0*Pi*(double)i/(double)n;
    v = (4.0/Pi)*(cos(th) - cos(3.0*th)/3.0 + cos(5.0*th)/5.0 - cos(7.0*th)/7.0
                  + cos(9.0*th)/9.0 - cos(11.0*th)/11.0 + cos(13.0*th)/13.0
                  - cos(15.0*th)/15.0 + cos(17.0*th)/17.0);
    printf("[%d] real = %g \n", i, v) ;
  }
  fftd(rslt, rdat, idat, n, 1);
  printf("square wave FFT cos series \n");
  for(i=0; i<n; i++) printf("cos %d th = %g, sin %d th = %g \n",
                            i, rslt[i*2], i, rslt[i*2+1]) ;
  printf("antialias and normalize square wave FFT cos\n");
  for(i=0; i<2*n; i++)arslt[i]=rslt[i];
  for(i=1; i<n/2; i++)
  {
    arslt[i*2]  = (rslt[i*2]+rslt[(n-i)*2])/(double)n;     /* real part */
    arslt[i*2+1]=-(rslt[i*2+1]-rslt[(n-i)*2+1])/(double)n; /* conj imag part */
  }
  arslt[0]=arslt[0]/(double)n;
  arslt[1]=-arslt[1]/(double)n   ;
  for(i=0; i<=(n/2)-1; i++) printf("cos %d th = %g, sin %d th = %g \n",
                                   i, arslt[i*2], i, arslt[i*2+1]) ;
  printf("check cos 1 th = 1.27, cos 3 th = -0.42, cos 5 th = 0.25 \n");
  for(i=0; i<n; i++)
  {
    rdat[i] = rslt[2*i];
    idat[i] = rslt[2*i+1];
  }
  fftd(rslt, rdat, idat, n, -1); /* inverse FFT */
  printf("square wave inverse FFT, cos \n");
  for(i=0; i<n; i++) printf("rdat %d = %g \n", i, rslt[i*2]) ;
  printf("\n");
  printf("\n");


  for(i=0; i<n; i++) /* square wave 4/Pi(sin th + 1/3 sin 3 th + */
                     /*  1/5 sin 5 th + 1/7 sin 7 th */
  {
    v = 1.0;
    if(i>n/2) v = -1.0;
    if(i==n/2 || i==0) v = 0.0;
    rdat[i] = v;
    idat[i] = 0.0;
  }
  printf("square wave data, for sin series \n");
  for(i=0; i<n; i++) printf("[%d] real = %g, imag = %g \n",
                            i, rdat[i], idat[i]) ;
  printf("square wave generate, using sin \n");
  for(i=0; i<n; i++)
  {
    th = 2.0*Pi*(double)i/(double)n;
    v = (4.0/Pi)*(sin(th) + sin(3.0*th)/3.0 + sin(5.0*th)/5.0 + sin(7.0*th)/7.0
                  + sin(9.0*th)/9.0 + sin(11.0*th)/11.0 + sin(13.0*th)/13.0
                  + sin(15.0*th)/15.0 + sin(17.0*th)/17.0);
    printf("[%d] real = %g \n", i, v) ;
  }
  fftd(rslt, rdat, idat, n, 1);
  printf("square wave FFT sin series \n");
  for(i=0; i<n; i++) printf("cos %d th = %g, sin %d th = %g \n",
                            i, rslt[i*2], i, rslt[i*2+1]) ;
  printf("antialias and normalize square wave FFT sin\n");
  for(i=0; i<2*n; i++)arslt[i]=rslt[i];
  for(i=1; i<n/2; i++)
  {
    arslt[i*2]  = (rslt[i*2]+rslt[(n-i)*2])/(double)n;     /* real part */
    arslt[i*2+1]=-(rslt[i*2+1]-rslt[(n-i)*2+1])/(double)n; /* conj imag part */
  }
  arslt[0]=arslt[0]/(double)n;
  arslt[1]=-arslt[1]/(double)n   ;
  for(i=0; i<=(n/2)-1; i++) printf("cos %d th = %g, sin %d th = %g \n",
                                   i, arslt[i*2], i, arslt[i*2+1]) ;
  printf("check sin 1 th = 1.27, sin 3 th = 0.42, sin 5 th = 0.25 \n");
  for(i=0; i<n; i++)
  {
    rdat[i] = rslt[2*i];
    idat[i] = rslt[2*i+1];
  }
  fftd(rslt, rdat, idat, n, -1);
  printf("square wave inverse FFT, sin series \n");
  for(i=0; i<n; i++) printf("rdat %d = %g \n", i, rslt[i*2]) ;
  printf("\n");
  printf("\n");


  for(i=0; i<n; i++) /* square wave, no zero test */
  {
    v = 1.0;
    if(i>=n/2) v = -1.0;
    rdat[i] = v;
    idat[i] = 0.0;
  }
  printf("square wave data, for sin, no zero \n");
  printf("          note: cos values come in \n");
  for(i=0; i<n; i++) printf("[%d] real = %g, imag = %g \n",
                            i, rdat[i], idat[i]) ;
  fftd(rslt, rdat, idat, n, 1);
  printf("square wave FFT sin, no zero \n");
  for(i=0; i<n; i++) printf("cos %d th = %g, sin %d th = %g \n",
                            i, rslt[i*2], i, rslt[i*2+1]) ;
  printf("antialias and normalize square wave FFT sin\n");
  for(i=0; i<2*n; i++)arslt[i]=rslt[i];
  for(i=1; i<n/2; i++)
  {
    arslt[i*2]  = (rslt[i*2]+rslt[(n-i)*2])/(double)n;     /* real part */
    arslt[i*2+1]=-(rslt[i*2+1]-rslt[(n-i)*2+1])/(double)n; /* conj imag part */
  }
  arslt[0]=arslt[0]/(double)n;
  arslt[1]=-arslt[1]/(double)n   ;
  for(i=0; i<=(n/2)-1; i++) printf("cos %d th = %g, sin %d th = %g \n",
                                   i, arslt[i*2], i, arslt[i*2+1]) ;
  printf("check sin 1 th = 1.27, sin 3 th = -0.42, sin 5 th = 0.25 \n");
  for(i=0; i<n; i++)
  {
    rdat[i] = rslt[2*i];
    idat[i] = rslt[2*i+1];
  }
  fftd(rslt, rdat, idat, n, -1);
  printf("square wave inverse FFT, sin no zero \n");
  for(i=0; i<n; i++) printf("rdat %d = %g \n", i, rslt[i*2]) ;
  printf("\n");
  printf("\n");


  for(i=0; i<n; i++) /* triangle wave 8/Pi^2(cos th + 1/9 cos 3 th + */
                     /*    1/25 cos 5 th */
  {
    v = 1.0-4.0*(double)i/(double)n; /* not exactly balanced */
    if(i>=n/2) v = -3.0+4.0*(double)i/(double)n;
    rdat[i] = v;
    idat[i] = 0.0;
  }
  printf("triangle wave data, cos series \n");
  for(i=0; i<n; i++) printf("[%d] real = %g, imag = %g \n",
                            i, rdat[i], idat[i]) ;
  printf("triangle wave generate, cos series \n");
  for(i=0; i<n; i++)
  {
    th = 2.0*Pi*(double)i/(double)n;
    v = (8.0/(Pi*Pi))*(cos(th) + cos(3.0*th)/9.0 + cos(5.0*th)/25.0 +
                       cos(7.0*th)/49.0 + cos(9.0*th)/81.0 +
                       cos(11.0*th)/121.0 + cos(13.0*th)/169.0 +
                       cos(15.0*th)/225.0 + cos(17.0*th)/289.0 );
    printf("[%d] real = %g \n", i, v) ;
  }
  fftd(rslt, rdat, idat, n, 1);
  printf("triangle wave FFT, cos series \n");
  for(i=0; i<n; i++) printf("cos %d th = %g, sin %d th = %g \n",
                            i, rslt[i*2], i, rslt[i*2+1]) ;
  printf("antialias and normalize triangle wave FFT cos\n");
  for(i=0; i<2*n; i++)arslt[i]=rslt[i];
  for(i=1; i<n/2; i++)
  {
    arslt[i*2]  = (rslt[i*2]+rslt[(n-i)*2])/(double)n;     /* real part */
    arslt[i*2+1]=-(rslt[i*2+1]-rslt[(n-i)*2+1])/(double)n; /* conj imag part */
  }
  arslt[0]=arslt[0]/(double)n;
  arslt[1]=-arslt[1]/(double)n   ;
  for(i=0; i<=(n/2)-1; i++) printf("cos %d th = %g, sin %d th = %g \n",
                                   i, arslt[i*2], i, arslt[i*2+1]) ;
  printf("check cos 1 th = 0.81, cos 3 th = 0.09, cos 5 th = 0.03 \n");
  printf("\n");
  printf("\n");


  for(i=0; i<n; i++) /* triangle wave 8/Pi^2(sin th - 1/9 cos 3 th + */
                     /*    1/25 cos 5 th */
  {
    v = 4.0*(double)i/(double)n;
    if(i>n/4 && i<=3*n/4) v = 2.0-v;
    if(i>3*n/4) v = v -4.0;
    rdat[i] = v;
    idat[i] = 0.0;
  }
  printf("triangle wave data, sin series \n");
  for(i=0; i<n; i++) printf("[%d] real = %g, imag = %g \n",
                            i, rdat[i], idat[i]) ;
  printf("triangle wave generate, sin series \n");
  for(i=0; i<n; i++)
  {
    th = 2.0*Pi*(double)i/(double)n;
    v = (8.0/(Pi*Pi))*(sin(th) - sin(3.0*th)/9.0 + sin(5.0*th)/25.0 -
                       sin(7.0*th)/49.0 + sin(9.0*th)/81.0 -
                       sin(11.0*th)/121.0 + sin(13.0*th)/169.0 -
                       sin(15.0*th)/225.0 + sin(17.0*th)/289.0 );
    printf("[%d] real = %g \n", i, v) ;
  }
  fftd(rslt, rdat, idat, n, 1);
  printf("triangle wave FFT \n");
  for(i=0; i<n; i++) printf("cos %d th = %g, sin %d th = %g \n",
                            i, rslt[i*2], i, rslt[i*2+1]) ;
  printf("antialias and normalize triangle wave FFT sin\n");
  for(i=0; i<2*n; i++)arslt[i]=rslt[i];
  for(i=1; i<n/2; i++)
  {
    arslt[i*2]  = (rslt[i*2]+rslt[(n-i)*2])/(double)n;     /* real part */
    arslt[i*2+1]=-(rslt[i*2+1]-rslt[(n-i)*2+1])/(double)n; /* conj imag part */
  }
  arslt[0]=arslt[0]/(double)n;
  arslt[1]=-arslt[1]/(double)n   ;
  for(i=0; i<=(n/2)-1; i++) printf("cos %d th = %g, sin %d th = %g \n",
                                   i, arslt[i*2], i, arslt[i*2+1]) ;
  printf("check sin 1 th = 0.81, sin 3 th = -0.09, sin 5 th = 0.03 \n");
  printf("\n");
  printf("\n");


  for(i=0; i<n; i++) /* saw tooth wave 2/Pi(sin th - 1/2 sin 2 th + */
                     /*    1/3 sin 3 th - 1/4 sin 4 th */
  {
    v = 2.0*(double)i/(double)n; /* not exactly balanced */
    if(i>=n/2) v = -2.0+2.0*(double)i/(double)n;
    if(i==n/2) v = 0.0;
    rdat[i] = v;
    idat[i] = 0.0;
  }
  printf("saw tooth wave data, sin series \n");
  for(i=0; i<n; i++) printf("[%d] real = %g, imag = %g \n",
                            i, rdat[i], idat[i]) ;
  printf("sawtooth wave generate, sin series \n");
  for(i=0; i<n; i++)
  {
    th = 2.0*Pi*(double)i/(double)n;
    v =(2.0/Pi)*(sin(th) - sin(2.0*th)/2.0 + sin(3.0*th)/3.0 - sin(4.0*th)/4.0
              + sin(5.0*th)/5.0 - sin(6.0*th)/6.0 + sin(7.0*th)/7.0
              - sin(8.0*th)/8.0 + sin(9.0*th)/9.0 - sin(10.0*th)/10.0
              + sin(11.0*th)/11.0 - sin(12.0*th)/12.0 + sin(13.0*th)/13.0
              - sin(14.0*th)/14.0 + sin(15.0*th)/15.0 - sin(16.0*th)/16.0
              + sin(17.0*th)/17.0 - sin(18.0*th)/18.0 + sin(19.0*th)/19.0
              - sin(20.0*th)/20.0);
    printf("[%d] real = %g \n", i, v) ;
  }
  fftd(rslt, rdat, idat, n, 1);
  printf("saw tooth wave FFT, sin series \n");
  for(i=0; i<n; i++) printf("cos %d th = %g, sin %d th = %g \n",
                            i, rslt[i*2], i, rslt[i*2+1]) ;
  printf("antialias and normalize saw tooth wave FFT sin\n");
  for(i=0; i<2*n; i++)arslt[i]=rslt[i];
  for(i=1; i<n/2; i++)
  {
    arslt[i*2]  = (rslt[i*2]+rslt[(n-i)*2])/(double)n;     /* real part */
    arslt[i*2+1]=-(rslt[i*2+1]-rslt[(n-i)*2+1])/(double)n; /* conj imag part */
  }
  arslt[0]=arslt[0]/(double)n;
  arslt[1]=-arslt[1]/(double)n   ;
  for(i=0; i<=(n/2)-1; i++) printf("cos %d th = %g, sin %d th = %g \n",
                                   i, arslt[i*2], i, arslt[i*2+1]) ;
  printf("check sin 1 th = 0.64, sin 3 th = -0.32, sin 5 th = 0.21 \n");
  printf("\n");
  printf("\n");



  printf("\nset up convolution test \n");
  for(i=0; i<n; i++) {rdat[i+n]=0.0; idat[i+n]=0.0;}
  for(i=0; i<n; i++) /* square wave data, a */
  {
    v = 1.0;
    if(i>=n/4 && i<3*n/4) v = -1.0;
    if(i==n/4 || i==3*n/4) v = 0.0;
    rdat[i] = v;
    idat[i] = 0.0;
  }
  printf("square wave data, a , upper half zero\n");
  for(i=0; i<2*n; i++) printf("[%d] real = %g, imag = %g \n",
                            i, rdat[i], idat[i]) ;
  fftd(AT, rdat, idat, 2*n, 1);
  printf("AT FFT of sauare wave\n");
  for(i=0; i<2*n; i++) printf("cos %d th = %g, sin %d th = %g \n",
                            i, AT[2*i], i, AT[2*i+1]) ;

  for(i=0; i<n; i++) {rdat[i+n]=0.0; idat[i+n]=0.0;}
  for(i=0; i<n; i++) /* triangle wave data, b */
  {
    v = 1.0-4.0*(double)i/(double)n; /* not exactly balanced */
    if(i>=n/2) v = -3.0+4.0*(double)i/(double)n;
    rdat[i] = v;
    idat[i] = 0.0;
  }
  printf("triangle wave data, b , upper half zero \n");
  for(i=0; i<2*n; i++) printf("[%d] real = %g, imag = %g \n",
                            i, rdat[i], idat[i]) ;
  fftd(BT, rdat, idat, 2*n, 1);
  printf("BT FFT of sauare wave\n");
  for(i=0; i<2*n; i++) printf("cos %d th = %g, sin %d th = %g \n",
                            i, BT[2*i], i, BT[2*i+1]) ;
  /* form complex product, take inverse FFT */
  for(i=0; i<2*n; i++)
  {
    rdat[i] = AT[2*i]*BT[2*i] - AT[2*i+1]*BT[2*i+1]; 
    idat[i] = AT[2*i]*BT[2*i+1] + AT[2*i+1]*BT[2*i];
  }
  fftd(C, rdat, idat, 2*n, -1);
  printf("convolution of a with b \n");
  for(i=0; i<n; i++) printf("real[%d] = %g \n", i, C[2*i]) ;

  return 0;
} /* end test_fft.c */