splitfft.c

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

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

#define  M       9
#define  N       ( 1<<M )

void draw_image(double *x,int m,char *title1,char *title2,
		char *xdis1,char *xdis2,int dis_type);
void splitfft(double  *x,double *y,int length,int m);
/*********************************************************
*                                                        *
*  double *x  : 输入数据的实部缓冲区。当函数计算完后，将 *
*               FFT 的实部结果也存放在以x 开始的缓冲区中 *
*  double *y  : 输入数据的虚部缓冲区。当函数计算完后，将 *
*               FFT 的虚部结果也存放在以y 开始的缓冲区中 *
*  int length : FFT 分析的点数                           *
*                                                        *
*  int  m     : FFT 分析的阶数m和length要满足 length=2^m *
*                                                        *
**********************************************************/

void main()
{
  int     i;
  char    title[80],tmp[20];
  double  x[N+1],y[N+1];
  double  buffer[N+1];

	for(i=0;i<N;i++) {
	  x[i]=0; y[i]=0;
	}

	for(i=0;i<N/32;i++) {
	  x[i]=10; y[i]=0;
	}

	for(i=0;i<N;i++)  buffer[i]=x[i]*x[i]+y[i]*y[i];
	strcpy(title,"The Input Sampling Signal Data");
	draw_image(buffer,N,title,"The Magnitude","0",itoa(N,tmp,10),0);

	splitfft(x,y,N,M);

	for(i=0;i<N;i++)  buffer[i]=x[i]*x[i]+y[i]*y[i];

	strcpy(title,"The SplitFFT Result Data");
	draw_image(buffer,N,title,"The Magnitude","0",itoa(N,tmp,10),0);
}

void    splitfft(double  *x,double *y,int length,int m)
{
  int     i,j,k,step,start,index0,index1,index2,index3;
  int     n,sub_length;
  double   w,cc1,cc3,ss1,ss3,r1,r2,s1,s2,s3,temp;

  n=length;
  for(i=0;i<m-1;i++)  {
	sub_length=n/4;
	w=4.0*asin(1.0)/n;
	for(j=0;j<sub_length;j++) {
	    cc1=cos( j*w );         ss1=sin( j*w );
	    cc3=cos( j*w*3 );       ss3=sin( j*w*3 );
	    step=n*2;
	    start = j;
	    while( start < length ) {
		for( index0=start;index0<length;index0+=step) {
		   index1=index0+sub_length;
		   index2=index1+sub_length;
		   index3=index2+sub_length;
		   r1=x[index0]-x[index2];
		   x[index0]+=x[index2];
		   r2=x[index1]-x[index3];
		   x[index1]+=x[index3];

		   s1=y[index0]-y[index2];
		   y[index0]+=y[index2];
		   s2=y[index1]-y[index3];
		   y[index1]+=y[index3];

		   s3=r1-s2;
		   r1=r1+s2;
		   s2=r2-s1;

		   r2=r2+s1;
		   x[index2]=r1*cc1-s2*ss1;
		   y[index2]=-s2*cc1-r1*ss1;

		   x[index3]=s3*cc3+r2*ss3;
		   y[index3]=r2*cc3-s3*ss3;
		}
		start=2*step-n+j;
		step=4*step;
	    }
	}
	n=n/2;
  }

  start=0;step=4;
  while( start <length )  {
     for( index0=start; index0 <length ;index0+=step )	{
	 index1=index0+1;
	 r1=x[index0];
	 x[index0]=r1+x[index1];
	 x[index1]=r1-x[index1];

	 r1=y[index0];
	 y[index0]=r1+y[index1];
	 y[index1]=r1-y[index1];
     }
     start=step*2-2;
     step=step*4;
  }

  j=0;
  for(i=1;i<length-1;i++) {
     k=length/2;
     while( j >= k ) {
	j=j-k;
	k=k/2;
     }
     j=j+k;

     if( i < j ) {
	temp=x[i];
	x[i]=x[j];
	x[j]=temp;
	temp=y[i];
	y[i]=y[j];
	y[j]=temp;
     }
  }
}

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