pse.c

本文件包含三个小程序

    static int ready = 0;       /* flag to indicated stored value */
    static double gstore;        /* place to store other value */
    double v1,v2,r,fac,gaus;
    double uniform();

/* make two numbers if none stored */
    if(ready == 0) {
	do {
	    v1 = 2.*uniform();
	    v2 = 2.*uniform();
	    r = v1*v1 + v2*v2;
	} while(r > 1.0);       /* make radius less than 1 */

/* remap v1 and v2 to two Gaussian numbers */
	fac = sqrt(-2.*log(r)/r);
	gstore = v1*fac;        /* store one */
	gaus = v2*fac;          /* return one */
	ready = 1;              /* set ready flag */
    }

    else {
	ready = 0;      /* reset ready flag for next pair */
	gaus = gstore;  /* return the stored one */
    }

    return(gaus);
}

/**************************************************************************
uniform - 产生零均值,均匀分布的随机数(-0.5 到 0.5 之间)

函数返回一个零均值,均匀分布的 double型随机数.
无参数, 调用 rand()函数.

double uniform()
*************************************************************************/

/* this define gives the biggest value rand() will return and
is used only if not defined in a header file (non-ANSI
implementations of C may require this) */

#ifndef RAND_MAX
#define RAND_MAX 32767
#endif

double uniform(void)
{
    return((double)(rand() & RAND_MAX) / RAND_MAX - 0.5);
}

/*************************************************************************
ham - Hamming 窗函数

输入参数:
	  COMPLEX  *x : 输入复数数组指针;
	   int	    n : 输入数组长度.
输出参数:
	  加窗后结果放在输入数组中返回,
	  无输出参数.

void ham(COMPLEX *x, int n)
*************************************************************************/

void ham(COMPLEX *x,int n)
{
  int    i;
  double ham,factor;

  factor = 8.0*atan(1.0)/(n-1);
  for (i = 0 ; i < n ; i++){
    ham = 0.54 - 0.46*cos(factor*i);
    x->real *= ham;
    x->imag *= ham;
    x++;
  }
}

/*************************************************************************
han - Hanning 窗函数

输入参数:
	  COMPLEX  *x : 输入复数数组指针;
	   int	    n : 输入数组长度.
输出参数:
	  加窗后结果放在输入数组中返回,
	  无输出参数.

void han(COMPLEX *x, int n)
*************************************************************************/

void han(COMPLEX *x,int n)
{
  int    i;
  double factor,han;

  factor = 8.0*atan(1.0)/(n-1);
  for (i = 0 ; i < n ; i++){
    han = 0.5 - 0.5*cos(factor*i);
    x->real *= han;
    x->imag *= han;
    x++;
  }
}

/*************************************************************************
triang - triangle 窗函数

输入参数:
	  COMPLEX  *x : 输入复数数组指针;
	   int	    n : 输入数组长度.
输出参数:
	  加窗后结果放在输入数组中返回,
	  无输出参数.

void triang(COMPLEX *,int n)
*************************************************************************/

void triang(COMPLEX *x,int n)
{
  int    i;
  double tri,a;
  a = 2.0/(n-1);

  for (i = 0 ; i <= (n-1)/2 ; i++) {
    tri = i*a;
    x->real *= tri;
    x->imag *= tri;
    x++;
  }
  for ( ; i < n ; i++) {
    tri = 2.0 - i*a;
    x->real *= tri;
    x->imag *= tri;
    x++;
  }
}

/*************************************************************************
black - Blackman 窗函数

输入参数:
	  COMPLEX  *x : 输入复数数组指针;
	   int	    n : 输入数组长度.
输出参数:
	  加窗后结果放在输入数组中返回,
	  无输出参数.

void black(COMPLEX *x, int n)
*************************************************************************/

void black(COMPLEX *x,int n)
{
  int    i;
  double black,factor;

  factor = 8.0*atan(1.0)/(n-1);
  for (i=0; i<n; ++i){
    black = 0.42 - 0.5*cos(factor*i) + 0.08*cos(2*factor*i);
    x->real *= black;
    x->imag *= black;
    x++;
  }
}

/*************************************************************************
harris - 4 term Blackman-Harris 窗函数

输入参数:
	  COMPLEX  *x : 输入复数数组指针;
	   int	    n : 输入数组长度.
输出参数:
	  加窗后结果放在输入数组中返回,
	  无输出参数.

void harris(COMPLEX *x, int n)
*************************************************************************/

void harris(  COMPLEX *x,int n)
{
  int    i;
  double harris,factor,arg;

  factor = 8.0*atan(1.0)/n;
  for (i=0; i<n; ++i){
    arg = factor * i;
    harris = 0.35875 - 0.48829*cos(arg) + 0.14128*cos(2*arg)
	       - 0.01168*cos(3*arg);
    x->real *= harris;
    x->imag *= harris;
    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();
}