/**********************************************************************
FSAMPLE.C -  取樣混迭演示程序

dft          離散傅里葉變換函數(shù)
idft         IDFT 函數(shù)
draw_image   繪圖子程序
***********************************************************************/

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

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

#define    PI	(4.0*atan(1.0))

 void idft(COMPLEX *Datain, COMPLEX *Dataout, int N);
 void dft(COMPLEX *Datain, COMPLEX *Dataout, int N);
 void draw_image(double *x,int m,char *title1,char *title2,
		 char *xdis1,char *xdis2,int dis_type);

/********************************************************/
void main(void)
{
  int          i,j,k,length,m,type;
  char         title1[80],title2[80],tmp1[20],tmp2[20];
  double       far *amp;
  COMPLEX      far *fsamp, far *tsamp, far *ssamp;

  length = 256; /* length is the display sample size */

  amp = (double *) farcalloc(length+1,sizeof(double));
  ssamp = (COMPLEX *) farcalloc(length+1,sizeof(COMPLEX));
  fsamp = (COMPLEX *) farcalloc(length+1,sizeof(COMPLEX));
  tsamp = (COMPLEX *) farcalloc(length+1,sizeof(COMPLEX));
	if(!tsamp) {
	    printf("\nUnable to allocate complex array for calculation\n");
	    exit(1);
	}

  /* Input frequency sampling data for processing */
  for (i = 0; i < length; i++) {
     if(i<length/14) fsamp[i].real = 24.0-24.0*14.0*i/length;
     else if(i>length-length/14) fsamp[i].real = 24.0*14.0*i/length-24.0*14.0+24.0;
     else fsamp[i].real = 0;
     fsamp[i].imag=0;
  }

  printf("Generate the Signal. Waitting for the calculation...\n");
  idft(fsamp, ssamp, length);
  for(i = 0; i< length; i++)
     amp[i]=ssamp[i].real;

  strcpy(title1,"The Signal Data Sequency");
  strcpy(title2,"The Magnitude of Signal");
  strcpy(tmp1,"0");
  gcvt((float)length/2/70,5,tmp2);
  strcat(tmp2,"(s) t");
  draw_image(amp,length,title1,title2,tmp1,tmp2,0);

  /* Resample with fs = 10 Hz */
  for (i = 0 ; i < length/7 ; i++) {
    tsamp[i].real = ssamp[i*7].real;
    tsamp[i].imag = ssamp[i*7].imag;
  }

  printf("Waitting for the calculation...\n");
  dft(tsamp, fsamp, (int)(length/7));

  for(i = 0; i< length/7; i++)
     amp[i]=fsamp[i].real;

  strcpy(title1,"The Signal Spectrum fs=10Hz");
  strcpy(title2,"The Spectrum Magnitude");
  strcpy(tmp1,"0");
  strcpy(tmp2,"2*PI");
  draw_image(amp,(int)(length/7),title1,title2,tmp1,tmp2,0);

  /* Resample with fs = 7 Hz */
  for (i = 0 ; i < length/10 ; i++) {
    tsamp[i].real = ssamp[i*10].real;
    tsamp[i].imag = ssamp[i*10].imag;
  }

  printf("Waitting for the calculation...\n");
  dft(tsamp, fsamp, (int)(length/10));

  for(i = 0; i< length/10; i++)
     amp[i]=fsamp[i].real;

  strcpy(title1,"The Signal Spectrum fs=7Hz");
  strcpy(title2,"The Spectrum Magnitude");
  strcpy(tmp1,"0");
  strcpy(tmp2,"2*PI");
  draw_image(amp,(int)(length/10),title1,title2,tmp1,tmp2,0);

  /* Resample with fs = 14 Hz */
  for (i = 0 ; i < length/5 ; i++) {
    tsamp[i].real = ssamp[i*5].real;
    tsamp[i].imag = ssamp[i*5].imag;
  }

  printf("Waitting for the calculation...\n");
  dft(tsamp, fsamp, (int)(length/5));

  for(i = 0; i< length/5; i++)
     amp[i]=fsamp[i].real;

  strcpy(title1,"The Signal Spectrum fs=14Hz");
  strcpy(title2,"The Spectrum Magnitude");
  strcpy(tmp1,"0");
  strcpy(tmp2,"2*PI");
  draw_image(amp,(int)(length/5),title1,title2,tmp1,tmp2,0);

  farfree(fsamp);
  farfree(ssamp);
  farfree(tsamp);
  farfree(amp);
}

/***********************************************************************
dft -  離散傅里葉正變換子程序
輸入?yún)?shù):
	  COMPLEX *Datain : 輸入數(shù)據(jù)區(qū)指針;
	  COMPLEX *Dataout: 輸出數(shù)據(jù)區(qū)指針;
		   int  N : 數(shù)據(jù)長度;
輸出參數(shù):
	  輸出數(shù)據(jù)存放在 Dataout 所指的數(shù)據(jù)區(qū);
	  無輸出參數(shù).

void dft(COMPLEX *Datain, COMPLEX *Dataout, int N)
************************************************************************/
void dft(COMPLEX *Datain, COMPLEX *Dataout, int N)
{
    int i,k,n,p;
    static int nstore = 0;      /* store N for future use */
    static COMPLEX *cf;         /* coefficient storage */
    COMPLEX *cfptr,*Dinptr;
    double arg;

/* Create the coefficients if N has changed */

    if(N != nstore) {
	if(nstore != 0) free((char *) cf);    /* free previous */

	cf = (COMPLEX  *) calloc(N, sizeof(COMPLEX));
	if (!cf) {
	    printf("\nUnable to allocate memory for coefficients.\n");
	    exit(1);
	}

	arg = 8.0*atan(1.0)/N;
	for (i=0 ; i<N ; i++) {
	    cf[i].real = (float)cos(arg*i);
	    cf[i].imag = -(float)sin(arg*i);
	}
    }

/* Perform the DFT calculation */

    printf("\n");
    for (k=0 ; k<N ; k++) {

	Dinptr = Datain;
	Dataout->real = Dinptr->real;
	Dataout->imag = Dinptr->imag;
	Dinptr++;
	for (n=1; n<N; n++) {

	p = (int)((long)n*k % N);
	    cfptr = cf + p;         /* pointer to cf modulo N */

	    Dataout->real += Dinptr->real * cfptr->real
			     - Dinptr->imag * cfptr->imag;

	    Dataout->imag += Dinptr->real * cfptr->imag
			     + Dinptr->imag * cfptr->real;
	    Dinptr++;
	}
	if (k % 32 == 31) printf("*");
	Dataout++;          /* next output */
    }
    printf("\n");
}

/***********************************************************************
idft - 離散傅里葉反變換子程序
輸入?yún)?shù):
	  COMPLEX *Datain : 輸入數(shù)據(jù)區(qū)指針;
	  COMPLEX *Dataout: 輸出數(shù)據(jù)區(qū)指針;
		   int  N : 數(shù)據(jù)長度;
輸出參數(shù):
	  輸出數(shù)據(jù)存放在 Dataout 所指的數(shù)據(jù)區(qū);
	  無輸出參數(shù).

void idft(COMPLEX *Datain, COMPLEX *Dataout, int N)
************************************************************************/
void idft(COMPLEX *Datain, COMPLEX *Dataout, int N)
{
    int i,k,n,p;
    static int nstore = 0;      /* store N for future use */
    static COMPLEX *cf;         /* coefficient storage */
    COMPLEX *cfptr,*Dinptr;
    double arg;

/* Create the coefficients if N has changed */

    if(N != nstore) {
        if(nstore != 0) free((char *) cf);    /* free previous */

        cf = (COMPLEX  *) calloc(N, sizeof(COMPLEX));
        if (cf == 0) {
            printf("\nUnable to allocate memory for coefficients.\n");
            exit(1);
        }

/* scale stored values by 1/N */
        arg = 8.0*atan(1.0)/N;
	for (i=0 ; i<N ; i++) {
            cf[i].real = (float)(cos(arg*i)/(double)N);
            cf[i].imag = (float)(sin(arg*i)/(double)N);
        }
    }

/* Perform the DFT calculation */

    printf("\n");
    for (k=0 ; k<N ; k++) {

        Dinptr = Datain;
        Dataout->real = Dinptr->real * cf[0].real;
        Dataout->imag = Dinptr->imag * cf[0].real;
        Dinptr++;
        for (n=1; n<N; n++) {

        p = (int)((long)n*k % N);
            cfptr = cf + p;         /* pointer to cf modulo N */

	    Dataout->real += Dinptr->real * cfptr->real
                             - Dinptr->imag * cfptr->imag;

            Dataout->imag += Dinptr->real * cfptr->imag
                             + Dinptr->imag * cfptr->real;
            Dinptr++;
        }
        if (k % 32 == 31) printf("*");
	Dataout++;          /* next output */
    }
    printf("\n");
}

/************************************************************************
draw_image - 將輸入數(shù)據(jù)的幅度畫出圖形。該函數(shù)可自動(dòng)調(diào)整顯示的比例, 使圖形
	     充滿整個(gè)屏幕。

輸入?yún)?shù): double *x    -   輸入數(shù)據(jù)序列的指針;
	  int m        -   輸入數(shù)據(jù)序列的長度;
	  char *title1 -   顯示圖形的上標(biāo)題字符串指針;
	  char *xdis1  -   X 坐標(biāo)左邊顯示標(biāo)題字符串指針.
	  char *title2 -   顯示圖形的左標(biāo)題字符串指針.
	  char *xdis2  -   X 坐標(biāo)右邊顯示標(biāo)題字符串指針.
	  int dis_type -   顯示類型, 0:連線 1:直線.
輸出參數(shù): 無
*************************************************************************/
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();
}