/**************************************************************
* Program to demonstrate the Fast Fourier Transform Procedure *
*      (frequency analysis of a discreet signal F(t))         *
*                                                             *
*       C++ version with graphic results by J-P Moreau, Paris *
*                      (www.jpmoreau.fr)                      *
*                                                             *
*           (use with tgfft.mak, graph2d.cpp, basis_r.cpp(*)  *
*               and vmblock.cpp(*)).                          *
*                                                             *
* (*) to be found in MATRICES directory.                      *
* ----------------------------------------------------------- *
* SAMPLE DATA:                                                *
*                                                             *
* Input data file (tfft.dat):                                 *
*                                                             *
* (The test signal contains 3 frequencies: 50, 250, 500 hz)   *
*                                                             *
* 1024                                                        *
* 0.00000000000000E+0000   0.00000000000000E+0000             *
* 1.95503421309917E-0004   3.53914399999776E+0001             *
* 3.91006842619834E-0004   5.95684899999760E+0001             *
* 5.86510263929974E-0004   6.54621699999552E+0001             *
* 7.82013685239669E-0004   5.24038399999845E+0001             *
* ...                      ...                                *
* 1.98826979472187E-0001   2.77372500000183E-0001             *
* 1.99022482893497E-0001  -2.43361500000174E+0000             *
* 1.99217986314807E-0001  -4.84236799999780E+0000             *
* 1.99413489736116E-0001  -6.02247899999929E+0000             *
* 1.99608993157426E-0001  -5.45615399999951E+0000             *
* 1.99804496578736E-0001  -3.22824200000105E+0000             *
* 2.00000000000045E-0001  -2.96010699999982E-0003             *
*                                                             *
* Output file (tfft.lst):                                     *
*                                                             *
*   Frequency     Value                                       *
* ------------------------                                    *
*      0.00      0.271690                                     *
*      5.00      0.546336                                     *
*      9.99      0.555560                                     *
*     14.99      0.572242                                     *
*     19.98      0.598828                                     *
*     24.98      0.640061                                     *
*   ...          ...                                          *
*   2522.53      0.020520                                     *
*   2527.53      0.020518                                     *
*   2532.52      0.020516                                     *
*   2537.52      0.020513                                     *
*   2542.51      0.020512                                     *
*   2547.51      0.020511                                     *
*   2552.50      0.020511                                     *
*                                                             *
* A graph is displayed.                                       *
**************************************************************/
#include <stdio.h>
#include <math.h>
#include <basis.h>     // ABS, MACH_EPS, REAL, WriteEnd(), WriteHead()   
#include <vmblock.h>   // dynamic allocations
#include <windows.h>

HDC  hdc;
RECT rect;
HPEN hpen, hpenOld;

//Functions used here of Graph2D.cpp
void CourbeXY(int,int,float *,double,double);
void Legendes (int,char *,char *,char *,int,int,char *,int,int,char *);

//"home made" graphic commands for hdc environment used by above functions
void DrawPixel(int ix,int iy) {
  //sorry, no other available command found
  Rectangle(hdc,rect.left+ix,rect.top+iy,
	            rect.left+ix+2,rect.top+iy+1);
}

void Swap(int *i1,int *i2) {
  int it;
  it=*i1; *i1=*i2; *i2=it;
}

void DrawLine(int ix1,int iy1,int ix2,int iy2) {
  int i,il,ix,iy;
  double dx,dy;
  if (ix2<ix1) {
	Swap(&ix1,&ix2); Swap(&iy1,&iy2);
  }
  il=(int) sqrt((ix2-ix1)*(ix2-ix1)+(iy2-iy1)*(iy2-iy1));
  dx=(double) (ix2-ix1)/il; 
  dy=(double) (iy2-iy1)/il;
  ix=ix1; iy=iy1;
  for (i=1;i<il+1; i++) {
	DrawPixel(ix,iy);
	ix=(int) (ix1+i*dx); 
	iy=(int) (iy1+i*dy);
  }
}

void OutText(int ix,int iy,char *s) {
  RECT rct;
  rct.left=ix; rct.right=ix+10*strlen(s);
  rct.top=iy-5; rct.bottom=iy+15;
  DrawText(hdc, s, -1, &rct,
                DT_SINGLELINE | DT_CENTER | DT_VCENTER);
}


// **** Begin FFT section ****

FILE  *f_in,*f_out;
int   i,ip,n1,ndata;
REAL df,dt,f,temp,tbegin,tend,y;
float *signal;         // pointer to real vector
REAL  **A;             // pointer to matrix (complex vector)  
void  *vmblock;

 //complex numbers utility functions

 //returns real part of a complex number z
 double Re(double z[2])
 {
   return z[0];
 }

 //returns imaginary part of a complex number z
 double Im(double z[2])
 {
   return z[1];
 }

 //adds up to complex numbers (z1+z2)
 void Sum(double z1[2],double z2[2],double z1_plus_z2[2])
 {
   z1_plus_z2[0]=z1[0]+z2[0];
   z1_plus_z2[1]=z1[1]+z2[1];
 }

 //z1-z2
 void InvSum(double z1[2],double z2[2],double z1_z2[2])
 {
   z1_z2[0]=z1[0]-z2[0];
   z1_z2[1]=z1[1]-z2[1];
 }

 //z1*z2
 void Prod(double z1[2],double z2[2],double z1xz2[2])
 {
   double a1,a2,b1,b2;
   a1= Re(z1);
   b1= Im(z1);
   a2= Re(z2);
   b2= Im(z2);
   z1xz2[0]=a1*a2-b1*b2;
   z1xz2[1]=a1*b2+a2*b1;
 }

 /**************************************************************
 *             FAST FOURIER TRANSFORM PROCEDURE                *
 * ----------------------------------------------------------- *
 * This procedure calculates the fast Fourier transform of a   *
 * real function sampled in N points 0,1,....,N-1. N must be a *
 * power of two (2 power M).                                   *
 *  T being the sampling duration, the maximum frequency in    *
 * the signal can't be greater than fc = 1/(2T). The resulting *
 * specter H(k) is discreet and contains the frequencies:      *
 *         fn = k/(NT) with k = -N/2,.. -1,0,1,2,..,N/2.       *
 *         H(0) corresponds to null fréquency.                 *
 *         H(N/2) corresponds to fc frequency.                 * 
 * ----------------------------------------------------------- *
 * INPUTS:                                                     *
 *                                                             *
 *     A[i][2]  complex vector of size N, the real part of     *
 *              which contains the N sampled points of real    *
 *              signal to analyse (time spacing is constant).  *
 *                                                             *
 *          M   integer such as N=2^M                          *
 *                                                             *
 * OUTPUTS:                                                    *
 *                                                             *
 *     A[i][2]  complex vector of size N, the vector modulus   *
 *              contain the frequencies of input signal, the   *
 *              vector angles contain the corresponding phases *
 *              (not used here).                               *
 *                                                             *
 *                                     J-P Moreau/J-P Dumont   *
 **************************************************************/
void FFT(REAL **A, int M)  {
    
    const double pi=3.1415926535;
    double U[2],W[2],T[2];            //complex numbers
    int N,NV2,NM1,J,I,IP,K,L,LE,LE1;
    double Phi;

    N=(2 << (M-1));
    NV2=(N >> 1);
    NM1=N-1;
    J=1;
    for (I=1;I<=NM1;I++) {
      if (I<J) {
	    T[0]=A[J-1][0]; 
	    T[1]=A[J-1][1];
	    A[J-1][0]=A[I-1][0]; 
	    A[J-1][1]=A[I-1][1];
	    A[I-1][0]=T[0]; 
	    A[I-1][1]=T[1];
      }
      K=NV2;
      while (K<J) {
	    J=J-K;
	    K=(K >> 1);
      }
      J=J+K;
    }

    LE=1;
    for (L=1;L<=M;L++) {
      LE1=LE;
      LE=(LE << 1);
      U[0]=1.0;
      U[1]=0.0;
      Phi= pi/LE1;
      W[0] = cos(Phi);
      W[1] = sin(Phi);

      for (J=1;J<=LE1;J++)  {
	    I=J-LE;
	    while (I< N-LE) {
	      I=I+LE;
	      IP=I+LE1;
	      Prod(A[IP-1],U,T);
	      InvSum(A[I-1],T,A[IP-1]);
	      Sum(A[I-1],T,A[I-1]);
	    }
	    Prod(U,W,U);
      }
    }
    for (I=0;I<N;I++)  {
      A[I][0]=A[I][0]/N;
      A[I][1]=A[I][1]/N;
    }

} //FFT()


void Test_fft()  {
  extern MaxX, MaxY;
  //client drawing zone in pixels
  MaxX=(rect.right-rect.left);
  MaxY=(rect.bottom-rect.top);
  // open input and output file
  f_in=fopen("tfft.dat","r");
  f_out=fopen("tfft.lst","w");
  //read number of input signal points in input file
  fscanf(f_in,"%d",&ndata);
  //take nearest power of two
  if (ndata > 2048) {ndata=2048; ip=11;}
  if (ndata < 2048 && ndata>1023) {ndata=1024; ip=10;}
  if (ndata < 1024 && ndata>511) {ndata=512; ip=9;}
  if (ndata <  512 && ndata>255) {ndata=256; ip=8;}
  if (ndata <  256 && ndata>127) {ndata=128; ip=7;}
  if (ndata <  128 && ndata> 63) {ndata=64; ip=6;}
  if (ndata <   64 && ndata> 31) {ndata=32; ip=5;}
  if (ndata <   32) {
    fprintf(f_out," Error: number of points too small (<32) !\n");
    fclose(f_out);
    printf(" Results in file tfft.lst (error).\n");
  }
  //Allocate memory space for A and signal
  vmblock = vminit();
  A      = (REAL **) vmalloc(vmblock, MATRIX, ndata, 2);
  signal = (float *) vmalloc(vmblock, VEKTOR, ndata, 0);
  if (! vmcomplete(vmblock)) {
    LogError ("Memory full!", 0, __FILE__, __LINE__);
    return;
  }	  
  
  //read ndata couples T(i), Y(i) in input data file
  for (i=0; i<ndata; i++) {
    fscanf(f_in,"%lf %lf",&temp,&y);
    signal[i] = (float) y;
    if (i==0)  tbegin=temp;
    if (i==ndata-1)  tend=temp;
  }

  //Draw input signal F(t)
  CourbeXY(ndata,5,signal,tbegin,tend);
  Legendes (5,"  INPUT SIGNAL (TFFT.DAT)  ","Time s","Acc. m/s2",
	        MaxX-200,50,"1024 points",0,0,"");

  // Put input signal in real part of complex vector A
  for (i=0; i<ndata; i++) {
    A[i][0]=signal[i];
    A[i][1]=0.0;
  }
  // call FFT procedure
  FFT(A,ip);
  // get frequencies in signal real vector
  n1=ndata/2;
  for (i=0; i<n1; i++) {
    temp= sqrt(A[i][0]*A[i][0]+A[i][1]*A[i][1]);
    if (i>0)  temp=2*temp;
    signal[i]=(float) temp;
  }

  // calculate sampling time range dt
  dt=(tend-tbegin)/(ndata-1);
  // calculate frequency step df
  df= 1.0/(ndata*dt);
  // print frequency specter to output file
  f=0.0;
  fprintf(f_out,"   Frequency     Value   \n");
  fprintf(f_out," ------------------------\n");
  for (i=0; i<n1; i++) {
    fprintf(f_out,"  %8.2f  %12.6f\n",f,signal[i]);
    f+=df;
  }
  //restore last frequency
  f-=df;
  fclose(f_out);
  //Draw graph of Fourier analysis (left quarter)
  CourbeXY(n1/4,3,signal,0,f/4);
  free(vmblock);
  Legendes (3,"  FOURIER ANALYSIS  ","freq. hz","Acc. m/s2",
	        MaxX/2-175,(MaxY/2)+50,"Numerical results in file tfft.lst.",0,0,"");
}

// **** End FFT section ****


//Handle window Paint and Destroy messages
long FAR PASCAL /*_export*/ WndProc(HWND hwnd,
    UINT message, UINT wParam, LONG lParam) 
{
   
    PAINTSTRUCT ps;

    switch (message) {
        case WM_PAINT:
	    hdc = BeginPaint(hwnd, &ps);
	    GetClientRect(hwnd, &rect);
            hpen = CreatePen(PS_SOLID, 1, RGB(0, 0, 255));
            hpenOld = SelectObject(hdc,hpen);

            Test_fft();
               
            SelectObject(hdc,hpenOld);
            DeleteObject(hpen);
            EndPaint(hwnd, &ps);
            return 0;
        case WM_DESTROY:
            PostQuitMessage(0);
            return 0;
    }
    return DefWindowProc(hwnd, message, wParam, lParam);
}

// main program
int PASCAL WinMain(HANDLE hInstance, HANDLE hPrevInstance,
    LPSTR lpszCmdParam, int nCmdShow)
{
    static char szAppName[] = "Tgfft";
    HWND        hwnd;
    MSG         msg;
    WNDCLASS    wndclass;

    if (!hPrevInstance) {
        wndclass.style         = CS_HREDRAW | CS_VREDRAW;
        wndclass.lpfnWndProc   = WndProc;
        wndclass.cbClsExtra    = 0;
        wndclass.cbWndExtra    = 0;
        wndclass.hInstance     = hInstance;
        wndclass.hIcon         = LoadIcon(NULL, IDI_APPLICATION);
        wndclass.hCursor       = LoadCursor(NULL, IDC_ARROW);
        wndclass.hbrBackground = GetStockObject(WHITE_BRUSH);
        wndclass.lpszMenuName  = NULL;
        wndclass.lpszClassName = szAppName;

        RegisterClass(&wndclass);
    }

    hwnd = CreateWindow(szAppName,  // window class name
        " FAST FOURIER TRANSFORM",  // window caption
        WS_OVERLAPPEDWINDOW,        // window style                       
        CW_USEDEFAULT,              // initial x position
        CW_USEDEFAULT,              // initial y position
        CW_USEDEFAULT,              // initial x size
        CW_USEDEFAULT,              // initial y size
        NULL,                       // parent window handle
        NULL,                       // window menu handle
        hInstance,                  // program instance handle
        NULL);                      // creation parameters

   ShowWindow(hwnd, nCmdShow);
   UpdateWindow(hwnd);

   while (GetMessage(&msg, NULL, 0, 0)) {
       TranslateMessage(&msg);
       DispatchMessage(&msg);
   }
   return msg.wParam;
}

//end of file tgfft.cpp