iir.c

irr滤波器的实现与仿真供大家参考望有帮助

   p=PI*(0.5+(2.0*i+1.0)/2.0/(*n));
   (b+i)->real=wc*cos(p);
   (b+i)->image=wc*sin(p);
  }
}
printf("\\nThe order of prototype filter is:%d",*n);
/* b1=(COMPLEX*)calloc(*n+2,sizeof(COMPLEX));
b2=(COMPLEX*)calloc(*n+2,sizeof(COMPLEX));
h=(double *)calloc((*n+2),sizeof(double));

if(h==NULL){
  printf("\\nNot enough memory to allocate!");
  exit(0);
}

b1->real=-(b->real);b1->image=-(b->image);
(b1+1)->real=1.0;(b1+1)->image=0.0;
if(*n!=1){
  for(i=1;i<*n;i++){
   for(k=0;k<i;k++){
    cs1=(b1+k)->real-(b1+k+1)->real*(b+i)->real;
    cs2=(b1+k)->image-(b1+k+1)->real*(b+i)->image;
    (b2+k+1)->real=cs1+(b1+k+1)->image*(b+i)->image;
    (b2+k+1)->image=cs2-(b1+k+1)->image*(b+i)->real;
   }
   b2->real=-(b1->real*(b+i)->real-b1->image*(b+i)->image);
   b2->image=-(b1->real*(b+i)->image+b1->image*(b+i)->real);
   (b2+i+1)->real=((b1+i)->real);
   (b2+i+1)->image=((b1+i)->image);
   for(k=0;k<=i+1;k++){
    (b1+k)->real=(b2+k)->real;
    (b1+k)->image=(b2+k)->image;
    (b2+k)->real=0.0;

    (b2+k)->image=0.0;
   }
  }
}

for(i=0;i<=*n;i++)
  h[i]=(b1+i)->real;
for(i=0;i<=*n;i++)
  printf("\\nz[%2d]=%16f",i,h[i]);
printf("\\nz[0]=%16f,\\nz[1]=%16f,\\nz[2]=%16f,\\nz[3]=%16f,\\nz[4]=%16f",1.0,2.6131/wc,3.4142/pow(wc,2),2.6131/pow(wc,3),1/a);*/


for(i=0;i<=*n-1;i++)
  pp[i]=xs[*n-1][i];
for(i=0;i<=*n-1;i++)
  h[i]=pp[i]/pow(wc,i);

free(b);
/* free(b1);
free(b2);*/
return h;
}

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

  bsf-有理分式变换子程序

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

struct rptr * bsf(double *c,int ni,double *f1,double *f2,int nf,struct rptr *ptr,int *no)
{
int i,j,ii,nd,ne,ng;
double *d=NULL,*e=NULL,*g=NULL;
ptr->a=pnpe(f2,nf,ni,ptr->a,no);  /*calcute the numberator coefficients */
ptr->b=(double *)calloc(*no+2,sizeof(double));
if(ptr->b==NULL){
  printf("\\nNot enough memory to allocate!");
  exit(0);
}
for(i=0;i<=ni;i++){      /* calculate the denominator coefficients */
  d=pnpe(f1,nf,i,d,&nd);
  e=pnpe(f2,nf,ni-i,e,&ne);
  g=ypmp(d,nd,e,ne,g,&ng);
  for(j=0;j<=ng;j++)
   ptr->b[j]+=c[i]*g[j];
  free(d);
  free(e);
  free(g);
  }
  return ptr;
}
/*********************************************************************************/
/*pnpe-*/


/**********************************************************************************/
double *pnpe(double *a,int m,int n,double *b,int *mn)
{
  int i,j,k,nk;
  double *c;
  *mn=m*n;
  c=(double *)calloc(*mn+3,sizeof(double));
  b=(double *)calloc(*mn+3,sizeof(double));
  if(b==NULL){
   printf("\\nNot enough memory to allocate!");
   exit(0);
  }
  if(n==0) {*b=1.00;free(c);return b;}
  else{
   for(i=0;i<=m;i++) b[i]=a[i];
   if(n==1) {free(c);return b;}
   else {
    nk=m;
    for(i=1;i<n;i++){
     for(j=0;j<=m;j++)
      for(k=0;k<=nk;k++) c[k+j]+=a[j]*b[k];
     nk+=m;
     for(k=0;k<=nk;k++) {b[k]=c[k];c[k]=0.0;}
    }
   }
  }

  free(c);
  return b;
}


/******************************************************************************
ypmp-


********************************************************************************/
double *ypmp(double *a,int m,double *b,int n,double *c,int *mn)
{
int i,j;
*mn=m+n;
c=(double *)calloc(*mn+3,sizeof(double));
if(c==NULL){
  printf("\\nNot enough memory to allocate!");
  exit(0);
}


for(i=0;i<=m;i++)
  for(j=0;j<=n;j++)
   c[i+j]+=a[i]*b[j];
  return c;
}

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


draw_image
************************************************************/
void draw_image(double *x,int m,char *title1,char *title2,char *xdis1,char *xdis2,int dis_type)
{
int gdriver=DETECT,gmode,errorcode;
/* registerbgidriver(EGAVGA_driver);
initgraph(&gdriver,&gmode,"c:\\\\tc20");*/
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(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)){
  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,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 amplititude 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();
}