dfstest.c

这是数字信号处理方面的一些源码

/**********************************************************************

DFSTEST.C - 离散傅里叶级数合成连续周期信号

gen_freq   产生信号的频域数据
cal_time   计算信号的时域波形
draw_image 输入数据绘图子程序

***********************************************************************/

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

/* function prototypes and structures for fft functions */

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

#define    PI	(4.0*atan(1.0))

void gen_freq(COMPLEX *x, int m ,double a, double w0,int type);
void cal_time(COMPLEX *x, int m, 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       *amp;
  double       a,w0;
  COMPLEX      *fsamp,*tsamp;

  a = 1;         /* a : 信号幅度 */
  w0 = 2.0*PI*1; /* w0: 信号的频率(弧度) */
  m = 10;        /* m : DFS 长度 */
  length = 512;  /* length: 显示样点数据长度 */
  type = 3;      /* type: 对应于表2.2 中的信号类型(1,2,...,6) */

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

  /* Input frequency DFS sampling data for processing */
  gen_freq(fsamp,m,a,w0,type);

  for (i = 0; i < 2*m+1; i++) {
     amp[i]= sqrt(fsamp[i].real*fsamp[i].real+fsamp[i].imag*fsamp[i].imag);
  }

  strcpy(title1,"The DFS Data Sequency N=");
  strcat(title1,itoa(m,tmp1,10));
  strcpy(title2,"The Magnitude of DFS");
  strcpy(tmp1,"-N");
  strcpy(tmp2,"N");
  draw_image(amp,2*m,title1,title2,tmp1,tmp2,1);

/* Calculate the time domain signal wave of the input DFS*/
  printf("Waitting for the calculation...\n");
  cal_time(fsamp,m,tsamp,length);

  for (i = 0 ; i < length ; i++) {
    amp[i]   = tsamp[i].real;
  }
  strcpy(title1,"The Signal In Time Domain N=");
  strcat(title1,itoa(m,tmp1,10));
  strcpy(title2,"The Time Signal Magnitude");
  strcpy(tmp1,"0");
  gcvt(4.0*PI/w0,6,tmp2);
  strcat(tmp2,"(s)");
  draw_image(amp,length,title1,title2,tmp1,tmp2,0);

  free(fsamp);
  free(tsamp);
  free(amp);
}

/**************************************************************************
gen_freq - 信号的频率域数据产生子程序

输入参数:
	x : 复数指针,作为输出频域信号的存放区
	m : 整数,输出频域信号的长度(length=2^m+1)
	a : 信号的幅度, 双精度浮点数
	w0: 信号基波的频率(弧度),双精度浮点数
      type: 产生信号的类型, 对应于表2.2(1,2,...,6)
输出参数:
	产生的信号存放在 x 所指定的区域

void gen_freq(COMPLEX *x, int m ,double a, double w0,int type)
*************************************************************************/

void gen_freq(COMPLEX *x, int m, double a, double w0,int type)
{
    int i,j,le;

    le = 2*m + 1;
/* calculate the F-Domain X(nw0) values recursively */
    switch(type) {
       case 1:  i= -m;
		for (j=0; j<le; j++) {
		   if(abs(i%2) == 0) {
		      (x+j)->real = 2.0*a/(PI*(1.0-i*i));
		      (x+j)->imag = 0;
		   }
		   else {
		     (x+j)->real = 0;
		     (x+j)->imag = 0;
		   }
		   i++;
		}
		break;
       case 2:  i= -m;
		for (j=0; j<le; j++) {
		   if(abs(i%2) == 0) {
		      (x+j)->real = a/(PI*(1.0-i*i));
		      (x+j)->imag = 0;
		   }
		   else if(abs(i) != 1) {
		     (x+j)->real = 0;
		     (x+j)->imag = 0;
		   }
		   else {
		     (x+j)->real = 0;
		     (x+j)->imag = -1.0*i*a/4.0;
		   }
		   i++;
		}
		break;
       case 3:  i= -m;
		for (j=0; j<le; j++) {
		   if(abs(i)%2 == 0) {
		      (x+j)->real = 0;
		      (x+j)->imag = 0;
		   }
		   else if( abs(i>>1)%2 == 0 ) {
		     (x+j)->real = 2.0*a/i/PI;
		     (x+j)->imag = 0;
		   }
		   else {
		     (x+j)->real = -2.0*a/i/PI;
		     (x+j)->imag = 0;
		   }
		   i++;
		}
		break;
       case 4:  i= -m;
		for (j=0; j<le; j++) {
		   if(abs(i%2) == 0) {
		      (x+j)->real = 0;
		      (x+j)->imag = 0;
		   }
		   else {
		     (x+j)->real = 4.0*a*cos(i*PI/2)/(PI*PI*i*i);
		     (x+j)->imag = -4.0*a*sin(i*PI/2)/(PI*PI*i*i);
		   }
		   i++;
		}
		break;
       case 5:  i= -m;
		for (j=0; j<le; j++) {
		   if(i == 0) {
		      (x+j)->real = 0;
		      (x+j)->imag = 0;
		   }
		   else {
		     (x+j)->real = a*pow(-1,i+1)*cos(PI/2)/(i*PI);
		     (x+j)->imag = -a*pow(-1,i+1)*sin(PI/2)/(i*PI);
		   }
		   i++;
		}
		break;
       case 6:  for (j=0; j<le; j++) {
		  (x+j)->real = w0/2.0/PI;
		  (x+j)->imag = 0;
		}
		break;
       default: for (j=0; j<le; j++) {
		  (x+j)->real = 0;
		  (x+j)->imag = 0;
		}
		break;
    }
}

/***********************************************************************

cal_time - 计算信号的时间域波形子程序

输入参数:
	 x : 输入频率域的数据指针,复数指针;
	 m : 输入数据的长度, 整数;
    Dataout: 输出数据指针, 复数指针;
	 N : 输出数据长度;
输出参数:
	 输出数据在 Dataout 中, 长度为 N;
	 无输出参数.

void cal_time(COMPLEX *x, int m, COMPLEX *Dataout, int N)
*************************************************************************/

void cal_time(COMPLEX *x, int m, COMPLEX *Dataout, int N)
{
    int i,j,k;
    COMPLEX *Dinptr,*Doutptr;
    double arg;

/* Perform the DFT calculation */

    arg = 4*PI/N;
    printf("\n");
    Doutptr=Dataout;
    for (i=0; i<N; i++)  {
       Doutptr->real = 0;
       Doutptr->imag = 0;
       Doutptr++;
    }

    Dinptr = x;
    for (i=-m ; i<=m ; i++) {
	Doutptr=Dataout;
	for (j=0; j<N; j++) {
	   Doutptr->real += Dinptr->real*cos(arg*i*j)-Dinptr->imag*sin(arg*i*j);
	   Doutptr->imag += Dinptr->imag*cos(arg*i*j)+Dinptr->real*sin(arg*i*j);
	   Doutptr++;
	}
	printf("*");
	Dinptr++;
    }
    printf("\n");
}

/************************************************************************
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();
}