firfd.c

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

/**********************************************************************
*  FIR 滤波器频率取样法设计程序
*  程序运行中要求输入FIR滤波器的阶数N,及其频率域取样点的幅度值H(n).
/*********************************************************************/

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

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

#define	PI	(4.0*atan(1.0))

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

void idft(COMPLEX *Datain, COMPLEX *Dataout, int N);
void add_ph1(double *Datain, COMPLEX *Dataout, int N,int fif);
void add_ph2(double *Datain, COMPLEX *Dataout, int N,int fif);
void draw_image(double *x,int m,char *title1,char *title2,
		char *xdis1,char *xdis2,int dis_type);

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

void main(void)
{
 double *hwdb,*hm;
 COMPLEX *hk,*hw,*hn;
 int L,N,i,j,k,type;
 COMPLEX iw;
 double jw,p2n,cc,cc1,wp2n;
 char title1[80],title2[80],tmp[20];

 L=500;
 hwdb=(double*)calloc(L,sizeof(double));
 if(hwdb==NULL) {
       printf("\nNot enough memory to allocate!");
       exit(0);
 }

 do {
    printf("Input pointer Number N:");
    scanf("%d",&N);
 }while(N==0);
 hm=(double*)calloc(N,sizeof(double));
 hk=(COMPLEX*)calloc(N,sizeof(COMPLEX));
 hn=(COMPLEX*)calloc(N,sizeof(COMPLEX));
 hw=(COMPLEX*)calloc(L,sizeof(COMPLEX));
 if(hw==NULL) {
       printf("\nNot enough memory to allocate!");
       exit(0);
 }

 k=0;
 do {
    printf("Input H(%3d)=",k);
    scanf("%lf",hm+k);
    k++;
 }while((2*k/N)<1);

 if(((float)N/2.0-N/2!=0)&&(*hm!=0)) type=1;
 if(((float)N/2.0-N/2==0)&&(*hm!=0)) type=2;
 if(((float)N/2.0-N/2!=0)&&(*hm==0)) type=3;
 if(((float)N/2.0-N/2==0)&&(*hm==0)) type=4;
 if ((type==1)||(type==2)) add_ph1(hm,hk,N,type);
 else add_ph2(hm,hk,N,type);

 printf("\nPlease waitting for calculation...");
 idft(hk,hn,N);
 printf("\nThe result h(n) is:\n");
 j=0;
 for(i=0;i<N;i++) {
    j=j+1;
    if(j>22) {
       printf("Press any key to continue...\n");
       j=0;
       getch();
    }
    printf(" h(%3d)=%10f + j%10f\n",i,(hn+i)->real,(hn+i)->imag);
 }

 printf("\nPress any key to calculate the Frequency Response...");
 getch();
 printf("\nPlease waitting for calculation...");

 /* Calculate the magnitude-frequency response */
 for(i=0;i<L;i++)  {
    (hw+i)->real=0;
    (hw+i)->imag=0;
    jw=(float)i*PI/L;
    for(k=0;k<N;k++) {
      p2n=2*PI*(float)k/N;
      cc=sin((jw-p2n)/2.0);
      cc1=sin((double)N*(jw-p2n)/2.0);
      wp2n=-(jw-p2n)*(float)(N-1.0)/2.0;
      if(fabs(jw-p2n)>1.0e-80) {
	iw.real=cc1/cc/N*cos(wp2n);
	iw.imag=cc1/cc/N*sin(wp2n);
      }
      else {
	iw.real=cos(wp2n);
	iw.imag=sin(wp2n);
      }
      (hw+i)->real += (hk+k)->real*iw.real-(hk+k)->imag*iw.imag;
      (hw+i)->imag += (hk+k)->real*iw.imag+(hk+k)->imag*iw.real;
    }
    if (i % 21 == 20) printf("*");
 }

 for(i=0;i<L;i++)  {
    hwdb[i]=10*log10(pow((hw+i)->real,2)+pow((hw+i)->imag,2));
    if (hwdb[i]<-100) hwdb[i]=-100;
    if (i % 21 == 20) printf("*");
 }
 strcpy(title1,"FIR Transfer Property N=");
 strcpy(title2,"The Filter Attenuation (dB)");
 itoa(N,tmp,10);
 strcat(title1,tmp);
 strcpy(tmp,"PI(rad)");
 draw_image(hwdb,L,title1,title2,"0",tmp,0);

 free(hwdb);
 free(hw);
 free(hn);
 free(hk);
 free(hm);
}

/***********************************************************************
add_ph1 - 为频域取样数据增加相位因子(相位类型 1 或 3)

输入参数: double *Datain  : 输入数据区指针;
	  COMPLEX *Dataout: 输出数据区指针;
		 int   N  : 数据长度;
		 int  fif : 相位类型.
输出参数:
	  输出数据存放在 Dataout 所指的数据区.
	  无输出参数.

void add_ph1(double *Datain, COMPLEX *Dataout, int N,int fif)
************************************************************************/
void add_ph1(Datain,Dataout,N,fif)
 double *Datain;
 COMPLEX *Dataout;
 int N,fif;
{
 int i,nn;
 if((float)N/2.0-N/2==0) nn=N/2;
 else nn=(N-1)/2;
 for(i=nn+1;i<N;i++)  {
    if(fif==1) Datain[i]=Datain[N-i];
    else Datain[i]=-Datain[N-i];
 }
 for(i=0;i<N;i++) {
    (Dataout+i)->real=Datain[i]*cos(-(N-1)*PI*i/N);
    (Dataout+i)->imag=Datain[i]*sin(-(N-1)*PI*i/N);
 }
}

/***********************************************************************
add_ph2 - 为频域取样数据增加相位因子(相位类型 2 或 4)

输入参数: double *Datain  : 输入数据区指针;
	  COMPLEX *Dataout: 输出数据区指针;
		 int   N  : 数据长度;
		 int  fif : 相位类型.
输出参数:
	  输出数据存放在 Dataout 所指的数据区.
	  无输出参数.

void add_ph2(double *Datain, COMPLEX *Dataout, int N,int fif)
************************************************************************/

void add_ph2(Datain,Dataout,N,fif)
 double *Datain;
 COMPLEX *Dataout;
 int N,fif;
{
 int i,nn;
 if((float)N/2.0-N/2==0) nn=N/2;
 else nn=(N-1)/2;
 for(i=nn+1;i<N;i++)  {
    if(fif==4) Datain[i]=Datain[N-i];
    else Datain[i]=-Datain[N-i];
 }
 for(i=0;i<N;i++) {
    (Dataout+i)->real=Datain[i]*cos(-(N-1)*PI*i/N-PI/2);
    (Dataout+i)->imag=Datain[i]*sin(-(N-1)*PI*i/N-PI/2);
 }
}
/***********************************************************************
idft - 离散傅里叶反变换子程序
输入参数:
	  COMPLEX *Datain : 输入数据区指针;
	  COMPLEX *Dataout: 输出数据区指针;
		   int  N : 数据长度;
输出参数:
	  输出数据存放在 Dataout 所指的数据区;
	  无输出参数.

void idft(COMPLEX *Datain, COMPLEX *Dataout, int N)
************************************************************************/

void idft(Datain, Dataout, N)
    COMPLEX *Datain, *Dataout;
    int N;
{
    int i,k,n,p;
    static int nstore = 0;      /* store N for future use */
    static COMPLEX *cf;         /* coefficient storage */
    COMPLEX *cfptr,*Dinptr;
    double arg;

/* Create the coefficients if N has changed */

    if(N != nstore) {
        if(nstore != 0) free((char *) cf);    /* free previous */

        cf = (COMPLEX  *) calloc(N, sizeof(COMPLEX));
	if (cf == 0) {
            printf("\nUnable to allocate memory for coefficients.\n");
            exit(1);
        }

/* scale stored values by 1/N */
	arg = 8.0*atan(1.0)/N;      /* 2*PI/N */
        for (i=0 ; i<N ; i++) {
            cf[i].real = (float)(cos(arg*i)/(double)N);
            cf[i].imag = (float)(sin(arg*i)/(double)N);
        }
    }

/* Perform the DFT calculation */

    printf("\n");
    for (k=0 ; k<N ; k++) {

	Dinptr = Datain;
	Dataout->real = Dinptr->real * cf[0].real;
	Dataout->imag = Dinptr->imag * cf[0].real;
        Dinptr++;
        for (n=1; n<N; n++) {

        p = (int)((long)n*k % N);
            cfptr = cf + p;         /* pointer to cf modulo N */

            Dataout->real += Dinptr->real * cfptr->real
                             - Dinptr->imag * cfptr->imag;

	    Dataout->imag += Dinptr->real * cfptr->imag
			     + Dinptr->imag * cfptr->real;
	    Dinptr++;
	}
	if (k % 32 == 31) printf("*");
	Dataout++;          /* next output */
    }
}

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