rad4fft.c

DSP信号处理源码,包括数字信号处理课程中的基本源程序。

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

rad4fft    基4 DIT FFT
rad4ifft   基4 DIT 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 rad4fft(COMPLEX *,int);
 void rad4ifft(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=4;
  length = pow(4,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");
  rad4fft(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");
  rad4ifft(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);
}

/************************************************************************
rad4fft - 基4 DIT FFT 函数

输入参数:
	 COMPLEX *x : FFT 输入和输出数据区指针;
	      int m : FFT 长度 ( length = 2^m );
输出参数:
	 输出数据放在 x 所指的输入数据区.
	 无输出参数.

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

void rad4fft(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 u0,u1,u2,u3,temp,temp0,temp1,temp2,temp3;
    COMPLEX *xi,*xj,*xip,*xip0,*xip1,*xip2,*xip3,tm;
    COMPLEX *wptr0,*wptr1,*wptr2,*wptr3;

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

    double arg;

    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 = 4**m = fft length */

	n = pow(4,m);
	le = n*3/4;

/* allocate the storage for w */

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

/* calculate the w values recursively */

	arg = 2.0*4.0*atan(1.0)/n;         /* 2*PI/n calculation */
	xj = w;
	for (j = 0 ; j < le ; j++) {
	    xj->real = cos(arg*j);
	    xj->imag = -sin(arg*j);
	    xj++;
	}
    }

/* rearrange data by 4 reversing */

    for (i=0; i<n; i++) {
       j = 0;
       l = i;
       for (k = 0; k < m; ++k) {
	  j = (j<<2)|(l%4);
	  l = l/4;
       }
       if (i < j) {
	  xi = x + i;
	  xj = x + j;
	  temp = *xj;
	  *xj = *xi;
	  *xi = temp;
       }
    }

/* start fft */
    for (l = 0 ; l < m ; l++) {
	le = pow(4,l);  /* number of butterfly per part */
	windex = pow(4,(m-l-1));
	for (j = 0 ; j < le ; j++) {
	    wptr0 = w;
	    wptr1 = w + windex*1*j;
	    wptr2 = w + windex*2*j;
	    wptr3 = w + windex*3*j;
	    u0 = *wptr0;
	    u1 = *wptr1;
	    u2 = *wptr2;
	    u3 = *wptr3;
	    for (i = j ; i < n ; i = i + 4*le) {
		xip0 = x + i;              /* find butterfly elements */
		xip1 = xip0 + le;
		xip2 = xip1 + le;
		xip3 = xip2 + le;

		/* multiple w factor */
		tm.real = xip0->real*u0.real - xip0->imag*u0.imag;
		tm.imag = xip0->real*u0.imag + xip0->imag*u0.real;
		*xip0 = tm;
		tm.real = xip1->real*u1.real - xip1->imag*u1.imag;
		tm.imag = xip1->real*u1.imag + xip1->imag*u1.real;
		*xip1 = tm;
		tm.real = xip2->real*u2.real - xip2->imag*u2.imag;
		tm.imag = xip2->real*u2.imag + xip2->imag*u2.real;
		*xip2 = tm;
		tm.real = xip3->real*u3.real - xip3->imag*u3.imag;
		tm.imag = xip3->real*u3.imag + xip3->imag*u3.real;
		*xip3 = tm;

		/* ( 1, 1, 1, 1) */
		temp0.real = xip0->real + xip1->real + xip2->real + xip3->real;
		temp0.imag = xip0->imag + xip1->imag + xip2->imag + xip3->imag;
		/* ( 1, -j, -1, j) */
		temp1.real = xip0->real + xip1->imag - xip2->real - xip3->imag;
		temp1.imag = xip0->imag - xip1->real - xip2->imag + xip3->real;
		/* ( 1, -1, 1, -1) */
		temp2.real = xip0->real - xip1->real + xip2->real - xip3->real;
		temp2.imag = xip0->imag - xip1->imag + xip2->imag - xip3->imag;
		/* ( 1, j, -1, -j) */
		temp3.real = xip0->real - xip1->imag - xip2->real + xip3->imag;
		temp3.imag = xip0->imag + xip1->real - xip2->imag - xip3->real;

		*xip0 = temp0;
		*xip1 = temp1;
		*xip2 = temp2;
		*xip3 = temp3;
	    }
	}
    }
}

/**************************************************************************
rad4ifft - 基4 DIT IFFT 函数

输入参数:
	 COMPLEX *x : IFFT 输入和输出数据区指针;
	      int m : IFFT 长度 ( length = 2^m );
输出参数:
	 输出数据放在 x 所指的输入数据区.
	 无输出参数.

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

void rad4ifft(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 u0,u1,u2,u3,temp,temp0,temp1,temp2,temp3,tm;
    COMPLEX *xi,*xj,*xip0,*xip1,*xip2,*xip3,*wptr0,*wptr1,*wptr2,*wptr3;

    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 = 4**m = inverse fft length */

	n = pow(4,m);
	le = n*3/4;

/* allocate the storage for w */

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

/* calculate the w values recursively */

	arg = 2.0*4.0*atan(1.0)/n;         /* 2*PI/n calculation */
	wrecur_real = w_real = cos(arg);
	wrecur_imag = w_imag = sin(arg);  /* opposite sign from fft */
	xj = w;
	xj->real = 1.0;
	xj->imag = 0.0;
	xj++;
	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;
	}
    }

/* rearrange data by 4 reversing */

    for (i=0; i<n; i++) {
       j = 0;
       l = i;
       for (k = 0; k < m; ++k) {
	  j = (j<<2)|(l%4);
	  l = l/4;
       }
       if (i < j) {
	  xi = x + i;
	  xj = x + j;
	  temp = *xj;
	  *xj = *xi;
	  *xi = temp;
       }
    }

/* start inverse fft */

    le = 1;
    windex = n;
    for (l = 0 ; l < m ; l++) {
	windex = windex>>2;
	for (j = 0 ; j < le ; j++) {
	    wptr0 = w;
	    wptr1 = w + windex*1*j;
	    wptr2 = w + windex*2*j;
	    wptr3 = w + windex*3*j;
	    u0 = *wptr0;
	    u1 = *wptr1;
	    u2 = *wptr2;
	    u3 = *wptr3;
	    for (i = j ; i < n ; i = i + 4*le) {
		xip0 = x + i;              /* find butterfly elements */
		xip1 = xip0 + le;
		xip2 = xip1 + le;
		xip3 = xip2 + le;

		/* multiple w factor */
		tm.real = xip0->real*u0.real - xip0->imag*u0.imag;
		tm.imag = xip0->real*u0.imag + xip0->imag*u0.real;
		*xip0 = tm;
		tm.real = xip1->real*u1.real - xip1->imag*u1.imag;
		tm.imag = xip1->real*u1.imag + xip1->imag*u1.real;
		*xip1 = tm;
		tm.real = xip2->real*u2.real - xip2->imag*u2.imag;
		tm.imag = xip2->real*u2.imag + xip2->imag*u2.real;
		*xip2 = tm;
		tm.real = xip3->real*u3.real - xip3->imag*u3.imag;
		tm.imag = xip3->real*u3.imag + xip3->imag*u3.real;
		*xip3 = tm;

		/* ( 1, 1, 1, 1) */
		temp0.real = xip0->real + xip1->real + xip2->real + xip3->real;
		temp0.imag = xip0->imag + xip1->imag + xip2->imag + xip3->imag;
		/* ( 1, j, -1, -j) */
		temp1.real = xip0->real - xip1->imag - xip2->real + xip3->imag;
		temp1.imag = xip0->imag + xip1->real - xip2->imag - xip3->real;
		/* ( 1, -1, 1, -1) */
		temp2.real = xip0->real - xip1->real + xip2->real - xip3->real;
		temp2.imag = xip0->imag - xip1->imag + xip2->imag - xip3->imag;
		/* ( 1, -j, -1, j) */
		temp3.real = xip0->real + xip1->imag - xip2->real - xip3->imag;
		temp3.imag = xip0->imag - xip1->real - xip2->imag + xip3->real;

		*xip0 = temp0;
		*xip1 = temp1;
		*xip2 = temp2;
		*xip3 = temp3;
	    }
	}
	le = le<<2;
    }

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