/**********************************************************************
DIF_FFT.C - DIF FFT 測試程序

fft        基2 DIF FFT 子程序
ifft       基2 DIF IFFT 子程序
draw_image 繪圖子程序
***********************************************************************/

#include    <math.h>
#include    <stdlib.h>
#include    <stdio.h>
#include    <string.h>
#include    <conio.h>
#include    <graphics.h>

/* COMPLEX STRUCTURE */
typedef struct {
    float real, imag;
} COMPLEX;

#define    PI	(4.0*atan(1.0))

 void fft(COMPLEX *,int);
 void ifft(COMPLEX *,int);
 void draw_image(double *x,int m,char *title1,char *title2,
		char *xdis1,char *xdis2,int dis_type);

/********************************************************/
void main(void)
{
  int          i,length,m,j;
  char         title[80],tmp[20];
  double       *amp;
  double       a,tempflt;
  COMPLEX      *samp;

  m=8;
  length = 1<<m;
  amp = (double *) calloc(length+1,sizeof(double));
  samp = (COMPLEX *) calloc(length+1,sizeof(COMPLEX));
	if(!samp) {
	    printf("\nUnable to allocate complex array for fft\n");
	    exit(1);
	}

  /* Input sampling data for processing */

  printf("Waitting for sampling data...");
  for (i = 0; i < length; i++) {
      amp[i] = cos(2.0*PI*(double)i*2.0/length);
      printf("*");
  }
  for (i = 0; i < length; i++) {
     samp[i].real = amp[i]*cos(16.0*PI*(double)i*2.0/length)*amp[i];
     amp[i]= samp[i].real;
     printf("*");
  }

  strcpy(title,"The Sampling Signal Data");
  draw_image(amp,length,title,"The Magnitude","0",itoa(length,tmp,10),0);

/* Find the spectrum of the data */
  printf("Waitting for the FFT calculation...\n");
  fft(samp,m);

/* calculate magnitude */
  tempflt=0;
  for (i = 0 ; i < length ; i++) {
    tempflt  = samp[i].real * samp[i].real + samp[i].imag * samp[i].imag;
    amp[i]   = tempflt;
  }

  strcpy(title,"The Signal Freq Magnitude Result");
  draw_image(amp,length,title,"The Magnitude","0",itoa(length,tmp,10),0);

  printf("Waitting for the IFFT calculation...\n");
  ifft(samp,m);
  strcpy(title,"The IFFT Result");
  for (i=0;i<length;i++) amp[i]=samp[i].real;
  draw_image(amp,length,title,"The Magnitude","0",itoa(length,tmp,10),0);

  free(samp);
  free(amp);
}

/************************************************************************
fft - 基2 DIF FFT 子程序

輸入參數:
	 COMPLEX *x : FFT 輸入和輸出數據區指針;
	      int m : FFT 長度 ( length = 2^m );
輸出參數:
	 輸出數據放在 x 所指的輸入數據區.
	 無輸出參數.

void fft(COMPLEX *x, int m)
*************************************************************************/

void fft(COMPLEX *x,int m)
{
    static COMPLEX *w;           /* used to store the w complex array */
    static int mstore = 0;       /* stores m for future reference */
    static int n = 1;            /* length of fft stored for future */

    COMPLEX u,temp,tm;
    COMPLEX *xi,*xip,*xj,*wptr;

    int i,j,k,l,le,windex;

    double arg,w_real,w_imag,wrecur_real,wrecur_imag,wtemp_real;

    if(m != mstore) {

/* free previously allocated storage and set new m */

	if(mstore != 0) free(w);
	mstore = m;
	if(m == 0) return;       /* if m=0 then done */

/* n = 2**m = fft length */

	n = 1 << m;
	le = n/2;

/* allocate the storage for w */

	w = (COMPLEX *) calloc(le-1,sizeof(COMPLEX));
	if(!w) {
	    printf("\nUnable to allocate complex W array\n");
	    exit(1);
	}

/* calculate the w values recursively */

	arg = 4.0*atan(1.0)/le;         /* PI/le calculation */
	wrecur_real = w_real = cos(arg);
	wrecur_imag = w_imag = -sin(arg);
	xj = w;
	for (j = 1 ; j < le ; j++) {
	    xj->real = (float)wrecur_real;
	    xj->imag = (float)wrecur_imag;
	    xj++;
	    wtemp_real = wrecur_real*w_real - wrecur_imag*w_imag;
	    wrecur_imag = wrecur_real*w_imag + wrecur_imag*w_real;
	    wrecur_real = wtemp_real;
	}
    }

/* start fft */

    le = n;
    windex = 1;
    for (l = 0 ; l < m ; l++) {
	le = le/2;

/* first iteration with no multiplies */

	for(i = 0 ; i < n ; i = i + 2*le) {
	    xi = x + i;
	    xip = xi + le;
	    temp.real = xi->real + xip->real;
	    temp.imag = xi->imag + xip->imag;
	    xip->real = xi->real - xip->real;
	    xip->imag = xi->imag - xip->imag;
	    *xi = temp;
	}

/* remaining iterations use stored w */

	wptr = w + windex - 1;
	for (j = 1 ; j < le ; j++) {
	    u = *wptr;
	    for (i = j ; i < n ; i = i + 2*le) {
		xi = x + i;
		xip = xi + le;
		temp.real = xi->real + xip->real;
		temp.imag = xi->imag + xip->imag;
		tm.real = xi->real - xip->real;
		tm.imag = xi->imag - xip->imag;
		xip->real = tm.real*u.real - tm.imag*u.imag;
		xip->imag = tm.real*u.imag + tm.imag*u.real;
		*xi = temp;
	    }
	    wptr = wptr + windex;
	}
	windex = 2*windex;
    }

/* rearrange data by bit reversing */

    j = 0;
    for (i = 1 ; i < (n-1) ; i++) {
	k = n/2;
	while(k <= j) {
	    j = j - k;
	    k = k/2;
	}
	j = j + k;
	if (i < j) {
	    xi = x + i;
	    xj = x + j;
	    temp = *xj;
	    *xj = *xi;
	    *xi = temp;
	}
    }
}

/************************************************************************
ifft - 基2 DIF IFFT 子程序

輸入參數:
	 COMPLEX *x : IFFT 輸入和輸出數據區指針;
	      int m : IFFT 長度 ( length = 2^m );
輸出參數:
	 輸出數據放在 x 所指的輸入數據區.
	 無輸出參數.

void ifft(COMPLEX *x, int m)
*************************************************************************/

void ifft(x,m)
    COMPLEX *x;
    int m;
{
    static COMPLEX *w;           /* used to store the w complex array */
    static int mstore = 0;       /* stores m for future reference */
    static int n = 1;            /* length of ifft stored for future */

    COMPLEX u,temp,tm;
    COMPLEX *xi,*xip,*xj,*wptr;

    int i,j,k,l,le,windex;

    double arg,w_real,w_imag,wrecur_real,wrecur_imag,wtemp_real;
    float scale;

    if(m != mstore) {

/* free previously allocated storage and set new m */

	if(mstore != 0) free(w);
	mstore = m;
	if(m == 0) return;       /* if m=0 then done */

/* n = 2**m = inverse fft length */

	n = 1 << m;
	le = n/2;

/* allocate the storage for w */

	w = (COMPLEX *) calloc(le-1,sizeof(COMPLEX));
	if(!w) {
	    printf("\nUnable to allocate complex W array\n");
	    exit(1);
	}

/* calculate the w values recursively */

	arg = 4.0*atan(1.0)/le;         /* PI/le calculation */
	wrecur_real = w_real = cos(arg);
	wrecur_imag = w_imag = sin(arg);  /* opposite sign from fft */
	xj = w;
	for (j = 1 ; j < le ; j++) {
	    xj->real = (float)wrecur_real;
	    xj->imag = (float)wrecur_imag;
	    xj++;
	    wtemp_real = wrecur_real*w_real - wrecur_imag*w_imag;
	    wrecur_imag = wrecur_real*w_imag + wrecur_imag*w_real;
	    wrecur_real = wtemp_real;
	}
    }

/* start inverse fft */

    le = n;
    windex = 1;
    for (l = 0 ; l < m ; l++) {
	le = le/2;

/* first iteration with no multiplies */

	for(i = 0 ; i < n ; i = i + 2*le) {
	    xi = x + i;
	    xip = xi + le;
	    temp.real = xi->real + xip->real;
	    temp.imag = xi->imag + xip->imag;
	    xip->real = xi->real - xip->real;
	    xip->imag = xi->imag - xip->imag;
	    *xi = temp;
	}

/* remaining iterations use stored w */

	wptr = w + windex - 1;
	for (j = 1 ; j < le ; j++) {
	    u = *wptr;
	    for (i = j ; i < n ; i = i + 2*le) {
		xi = x + i;
		xip = xi + le;
		temp.real = xi->real + xip->real;
		temp.imag = xi->imag + xip->imag;
		tm.real = xi->real - xip->real;
		tm.imag = xi->imag - xip->imag;
		xip->real = tm.real*u.real - tm.imag*u.imag;
		xip->imag = tm.real*u.imag + tm.imag*u.real;
		*xi = temp;
	    }
	    wptr = wptr + windex;
	}
	windex = 2*windex;
    }

/* rearrange data by bit reversing */

    j = 0;
    for (i = 1 ; i < (n-1) ; i++) {
	k = n/2;
	while(k <= j) {
	    j = j - k;
	    k = k/2;
	}
	j = j + k;
	if (i < j) {
	    xi = x + i;
	    xj = x + j;
	    temp = *xj;
	    *xj = *xi;
	    *xi = temp;
	}
    }

/* scale all results by 1/n */
    scale = (float)(1.0/n);
    for(i = 0 ; i < n ; i++) {
	x->real = scale*x->real;
	x->imag = scale*x->imag;
	x++;
    }
}

/************************************************************************
draw_image - 將輸入數據的幅度畫出圖形。該函數可自動調整顯示的比例, 使圖形
	     充滿整個屏幕。

輸入參數: double *x    -   輸入數據序列的指針;
	  int m        -   輸入數據序列的長度;
	  char *title1 -   顯示圖形的上標題字符串指針;
	  char *xdis1  -   X 坐標左邊顯示標題字符串指針.
	  char *title2 -   顯示圖形的左標題字符串指針.
	  char *xdis2  -   X 坐標右邊顯示標題字符串指針.
	  int dis_type -   顯示類型, 0:連線 1:直線.
輸出參數: 無
*************************************************************************/
void draw_image(double *x,int m,char *title1,char *title2,
		char *xdis1,char *xdis2,int dis_type)
{
 int gdriver=DETECT, gmode,errorcode;
 int i,scx,scy,y0,signa,signb;
 int style, userpat;
 int start_x=40,start_y=40,end_x=10,end_y=60;
 long tlen;
 double ys,xs,ym;
 char dis[40];
 /*initializes the graphics mode */
 initgraph(&gdriver,&gmode,"");
 errorcode=graphresult();
 if (errorcode != grOk) {
    printf("Graphics error: %s\n",grapherrormsg(errorcode));
    printf("Press any key to halt!\n");
    getch();
    exit(1);
 }
 scx=getmaxx();
 scy=getmaxy();
 ym=1.e-90;
 signa=0;
 signb=0;

 for(i=0;i<m;i++) {
    if ((*(x+i)>0)&&(*(x+i)>ym))  ym = *(x+i);
    if ((*(x+i)<0)&&(- *(x+i)>ym))  ym = - *(x+i);
 }
 for(i=0;i<m;i++)  {
    if (*(x+i)>fabs(ym/20)) signa=1;
    if (*(x+i)<-fabs(ym/20)) signb=1;
 }
 if ((signa==1)&&(signb==1)) ys=(double)((scy - start_y - end_y)>>1)/ym;
 else ys=(double)((scy - start_y - end_y)/ym);
 xs=(double)(scx - start_x - end_x)/m;
 y0=((scy - start_y - end_y)>>1)+start_y;

 /* draw the frame */

 setcolor(LIGHTGREEN);
 rectangle(start_x-1,start_y-20,scx-end_x+1,scy-end_y+20);

 setcolor(DARKGRAY);
 /* select the line style */
 style=DASHED_LINE;
 userpat = 1;
 setlinestyle(style, userpat, 1);
 /* a user defined line pattern */
 /* binary: "0000000000000001"  */
 for(i=0;i<=10;i++)
    line(start_x,start_y+(scy-start_y-end_y)*i/10,scx-end_x,start_y+(scy-start_y-end_y)*i/10);
 for(i=0;i<=10;i++)
    line(start_x+(scx-start_x-end_x)*i/10,start_y,start_x+(scx-start_x-end_x)*i/10,scy-end_y);
 setcolor(GREEN);
 style=SOLID_LINE;
 userpat = 1;
 setlinestyle(style, userpat, 1);
 rectangle(start_x,start_y,scx-end_x,scy-end_y);
 setcolor(YELLOW);
 for(i=0;i<=10;i++)
    line(start_x,start_y+(scy-start_y-end_y)*i/10,start_x+5,start_y+(scy-start_y-end_y)*i/10);
 for(i=0;i<=10;i++)
    line(start_x+(scx-start_x-end_x)*i/10,scy-end_y+15,start_x+(scx-start_x-end_x)*i/10,scy-end_y+20);
 settextstyle(DEFAULT_FONT,HORIZ_DIR,1);
 setcolor(YELLOW);
 if((signa==1)&&(signb==0)) {
    strcpy(dis,"0");
    outtextxy(start_x+2,scy-end_y+4,dis);
    gcvt(ym,5,dis);
    outtextxy(start_x+1,start_y-10,dis);
    outtextxy(start_x-10,scy-end_y+24,xdis1);
    outtextxy(scx-2-strlen(xdis2)*8,scy-end_y+24,xdis2);
 }
 else if((signb==1)&&(signa==0)) {
    strcpy(dis,"0");
    outtextxy(start_x+2,start_y-10,dis);
    gcvt(ym,5,dis);
    outtextxy(start_x+2,scy-end_y+4,"-");
    outtextxy(start_x+10,scy-end_y+4,dis);
    outtextxy(start_x-10,scy-end_y+24,xdis1);
    outtextxy(scx-2-strlen(xdis2)*8,scy-end_y+24,xdis2);
 }
 else {
    line(start_x,y0,scx-end_x,y0);
    strcpy(dis,"0");
    outtextxy(start_x-10,y0,dis);
    gcvt(ym,5,dis);
    outtextxy(start_x+2,start_y-10,dis);
    outtextxy(start_x+2,scy-end_y+4,"-");
    outtextxy(start_x+10,scy-end_y+4,dis);
    outtextxy(start_x-10,scy-end_y+24,xdis1);
    outtextxy(scx-2-strlen(xdis2)*8,scy-end_y+24,xdis2);
 }
 strcpy(dis,"Press any key to continue...");
 setcolor(LIGHTRED);
 outtextxy((scx-28*8)>>1,scy-16,dis);

 settextstyle(DEFAULT_FONT,HORIZ_DIR,2);
 tlen=strlen(title1);
 if ((tlen<<4)<scx) {
    setcolor(LIGHTGREEN);
    outtextxy((start_x+scx-end_x-(tlen<<4))>>1,start_y-40,title1);
 }

 settextstyle(DEFAULT_FONT,VERT_DIR,1);
 tlen=strlen(title2);
 if ((tlen<<4)<scy) {
    setcolor(LIGHTGREEN);
    outtextxy(start_x-20,(scy-end_y-(tlen<<3))>>1,title2);
 }
 /*draw the amplitude image*/
 setcolor(WHITE);
 if((signa==1)&&(signb==0)) y0=scy-end_y;
 else if((signb==1)&&(signa==0)) y0=start_y;
 if (dis_type == 0) {
    for(i=0;i<m-1;i++)
      line(xs*i+start_x,y0-*(x+i)*ys,xs*(i+1)+start_x,y0-*(x+i+1)*ys);
 }
 else if (dis_type == 1) {
    for(i=0;i<=m;i++)
      line(xs*i+start_x,y0-*(x+i)*ys,xs*i+start_x,y0);
 }
 getch();
 closegraph();
}