fastcor.c

本文件包含三个小程序

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++;
    }
}

/**************************************************************************
log2 - 基2求对数函数

函数返回以2为底的输入整数的对数值的整数.
如果对数值位于两个整数之间,则返回较大值.

int log2(unsigned int x)
*************************************************************************/

int log2(x)
    unsigned int x;
{
    unsigned int mask,i;

    if(x == 0) return(-1);      /* zero is an error, return -1 */

    x--;                        /* get the max index, x-1 */

    for(mask = 1 , i = 0 ; ; mask *= 2 , i++) {
	if(x == 0) return(i);   /* return log2 if all zero */
	x = x & (~mask);        /* AND off a bit */
    }
}

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