IIR数字滤波器设计C程序 - 包括图形显示[转帖] IIR, 图形, 数字滤波器, 转帖, 程序

yuepei

/**********************************************************
IIRD_DF.c-双线性变换法设计IIR数字滤波器
Infinite Impulse Response Digital Filter Design By Double Converting
====================================================
包括1.Lowpass       2.Highpass      3.bandpass三部分，如果有特殊需要
把其他部分删除掉即可。
***********************************************************/
/*#include <windows.h>*/
#include <math.h>
#include <dos.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <conio.h>
#include <graphics.h>

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

struct rptr{
double *a;
double *b;
};

#define  I     (double)(4.0*atan(1.0))
#define  FNSSH(x)    log(x+sqrt(x*x+1))
#define  FNCCH(x)    log(x+sqrt(x*x-1))
#define  FNSH1(x)    (exp(x)-exp(-x))/2
#define  FNCH1(x)    (exp(x)+exp(-x))/2

yuepei

/***********************************************************************************************/
double *bcg(double ap,double as,double wp,double ws,int *n,double *h,int *type);
struct rptr *bsf(double *c,int ni,double *f1,double *f2,int nf,struct rptr *ptr,int *no);
double *pnpe(double *a,int m,int n,double *b,int *mn);
double *ypmp(double *a,int m,double *b,int n,double *c,int *mn);
void lowpass_input(double *wp,double *ws,double *ap,double *ar,double *f1,double *f2,int *nf);
void highpass_input(double *wp,double *ws,double *ap,double *ar,double *f1,double *f2,int *nf);
void bandpass_input(double *wp,double *ws,double *ap,double *ar,double *f1,double *f2,int *nf);
void draw_image(double *x,int m,char *title1,char *title2,char *xdis1,char *xdis2,int dis_type);
/***********************************************************************************************/

void main(void)
{
double *f1,*f2,*hwdb;
double *h=NULL;
int N,nf,ns,nz,i,j,k,ftype,type;
COMPLEX hwdb1,hwdb2;
double wp,ws,ap,as,jw,amp1,amp2;
char title[80],tmp[20];
struct rptr *ptr=NULL;

/* printf("%lf",PI);
getch();*/
N=500;
f1=(double *)calloc(4,sizeof(double));
f2=(double *)calloc(4,sizeof(double));
hwdb=(double *)calloc(N+2,sizeof(double));
if (hwdb==NULL){
printf("\n Not enough memory to allocate!");
exit(0);
}
printf("\n1.Lowpass 2.Highpass 3.bandpass");
printf("\nPlease select the filter type:");
scanf("%d",&ftype);
switch(ftype){
case 1:lowpass_input(&wp,&ws,&ap,&as,f1,f2,&nf);
break;
case 2:highpass_input(&wp,&ws,&ap,&as,f1,f2,&nf);
break;
case 3:bandpass_input(&wp,&ws,&ap,&as,f1,f2,&nf);
break;
default:lowpass_input(&wp,&ws,&ap,&as,f1,f2,&nf);
}

h=bcg(ap,as,wp,ws,&ns,h,&type);
printf("\nThe analog filter denominator coefficients of Ha(s):");
for(i=0;i<=ns;i++)
printf("\nb[%2d]=%16f",i,h);

ptr=bsf(h,ns,f1,f2,nf,ptr,&nz);
printf("\nThe digital filter coefficients of H(z):");
printf("\n(a is numerator coefficient. b is denominator coefficient.)");
for(i=0;i<=nz;i++)
printf("\na[%2d]=%10f,b[%2d]=%16f",i,ptr->a,i,ptr->b);

printf("\n\nPress any key to calculate the filter response of H(z)...");
getch();
printf("\nWaitting for calculating...");


/* Calculate the magnitude-frequency response*/
for(k=0;k<=N;k++){
jw=k*PI/N;
hwdb1.real=0;hwdb1.image=0;
hwdb2.real=0;hwdb2.image=0;
for(i=0;i<=nz;i++){
hwdb1.real+=ptr->a*cos(i*jw);
hwdb2.real+=ptr->b*cos(i*jw);
hwdb1.image+=ptr->a*sin(i*jw);
hwdb2.image+=ptr->b*sin(i*jw);
}
amp1=(pow(hwdb1.real,2)+pow(hwdb1.image,2));
amp2=(pow(hwdb2.real,2)+pow(hwdb2.image,2));
if(amp1==0) amp1=1e-90;
if(amp2==0) amp2=1e-90;
hwdb[k]=10*log10(amp1/amp2);
if(hwdb[k]<-200) hwdb[k]=-200;
if(k%10==9) printf("*");
}

strcpy(title,"Transfer Property");
if(type==1) strcat(title,"(BW Type) N=");
if(type==2) strcat(title,"(CB Type) N=");
itoa(ns,tmp,10);
strcat(title,tmp);
strcpy(tmp,"PI(rad)");
draw_image(hwdb,N,title,"The Attenuation(dB)","0",tmp,0);

free(ptr->b);
free(ptr->a);
free(h);
free(hwdb);
free(f2);
free(f1);
}

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

yuepei

void lowpass_input(double *wp,double *ws,double *ap,double *ar,double *f1,double *f2,int *nf)
{
double fp,fr,fs;
printf("\nPlease input the Fp,Ap,Fr,Ar,Fs value");
printf("\nFp,Apassband frequency(Hz) and MAXAttenuation(dB)");
printf("\nFr,Ar:Stopband frequency(Hz) and MINAttenuation(dB)");
printf("\nFs is the sample frequency(Hz) Lowpass filter");
printf("\nInput parameters Fp,Ap,Fr,Ar,Fs:");
scanf("%lf,%lf,%lf,%lf,%lf",&fp,ap,&fr,ar,&fs);
if((fp>fr)||(*ap>*ar)||(fs<2*fr)){
do{
sound(1000);delay(200);nosound();
printf("Invalid input! Please Reinput:");
scanf("%lf,%lf,%lf,%lf,%lf",&fp,ap,&fr,ar,&fs);
}while((fp>fr)||(*ap>*ar)||(fs<2*fr));
}
*wp=tan(PI*fp/fs);
*ws=tan(PI*fr/fs);
*f1=-1.0;*(f1+1)=1.0;
*f2=1.0;*(f2+1)=1.0;
*nf=1;
}


/**********************************************************************************************/
void highpass_input(double *wp,double *ws,double *ap,double *ar,double *f1,double *f2,int *nf)
{
double fp,fr,fs;
printf("\nPlease input the Fp,Ap,Fr,Ar,Fs value");
printf("\nFp,Apassband frequency(Hz) and MAXAttenuation(dB)");
printf("\nFr,Ar:Stopband frequency(Hz) and MINAttenuation(dB)");
printf("\nFs is the sample frequency(Hz) highpass filter");
printf("\nInput parameters Fp,Ap,Fr,Ar,Fs:");
scanf("%lf,%lf,%lf,%lf,%lf",&fp,ap,&fr,ar,&fs);
if((fp*ar)||(fs<2*fp)){
do{
sound(1000);delay(200);nosound();
printf("Invalid input! Please Reinput:");
scanf("%lf,%lf,%lf,%lf,%lf",&fp,ap,&fr,ar,&fs);
}while((fp*ar)||(fs<2*fr));
}
*wp=fabs(1/tan(PI*fp/fs));
*ws=fabs(1/tan(PI*fr/fs));
*f1=1.0;*(f1+1)=1.0;
*f2=-1.0;*(f2+1)=1.0;
*nf=1;
}


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

yuepei

void bandpass_input(double *wp,double *ws,double *ap,double *ar,double *f1,double *f2,int *nf)
{
double fp1,fp2,fr1,fr2,fs,wp1,wp2,wr1,wr2,cw0,pwp1,pwp2,pws1,pws2;
printf("\nPlease input the Fp1,Fp2,Ap,Fr1,Fr2,,Ar,Fs value");
printf("\nFp1,Fp2,Apassband frequency(Hz) and MAXAttenuation(dB)");
printf("\nFr1,Fr2,Ar:Stopband frequency(Hz) and MINAttenuation(dB)");
printf("\nFs is the sample frequency(Hz) bandpass filter");
printf("\nInput parameters Fp1,Fp2,Ap,Fr1,Fr2,Ar,Fs:");
scanf("%lf,%lf,%lf,%lf,%lf,%lf,%lf",&fp1,&fp2,ap,&fr1,&fr2,ar,&fs);
if((fp1>fp2)||(fp1fr2)||(*ap>*ar)||(fs<2*fr2)){
do{
sound(1000);delay(200);nosound();
printf("Invalid input! Please Reinput:");
scanf("%lf,%lf,%lf,%lf,%lf,%lf,%lf",&fp1,&fp2,ap,&fr1,&fr2,ar,&fs);
}while((fp1>fp2)||(fp1fr2)||(*ap>*ar)||(fs<2*fr2));
}
wp1=2*PI*fp1/fs;wr1=2*PI*fr1/fs;
wp2=2*PI*fp2/fs;wr2=2*PI*fr2/fs;
cw0=sin(wp1+wp2)/(sin(wp1)+sin(wp2));
pwp1=fabs((cw0-cos(wp1))/sin(wp1));
pws1=fabs((cw0-cos(wr1))/sin(wr1));
pwp2=fabs((cw0-cos(wp2))/sin(wp2));
pws2=fabs((cw0-cos(wr2))/sin(wr2));
if(fabs(pws1-pwp1) *wp=pws1;
*ws=pws1;

}
else{
*wp=pwp2;
*ws=pws2;
}
*f1=1.0;*(f1+1)=-2.0*cw0;*(f1+2)=1.0;
*f2=-1.0;*(f2+1)=0.0*cw0;*(f2+2)=1.0;
*nf=2;
}

yuepei

/**************************************************************************
bcg-Chebyshev 和 Butterworth 型模拟原型传输函数生成子程序
即程序得到系统函数 H(s).
输出格式为：

****************************************************************************/
double *bcg(double ap,double as,double wp,double ws,int *n,double *h,int *type)
{
int i,j,k;
double a,c,e,p,q,x,y,wc,cs1,cs2;
COMPLEX *b;
double pp[20];
double xs[8][8]={{1.0},
{1.0,1.41421356},
{1.0,2.0,2.0},
{1.0,2.61312593,3.41421356,2,61312593},
{1.0,3.23606789,5.23606789,5.23606789,3.23606789},
{1.0,3.86370331,7.46410162,9.14162017,7.46410162,3.86370331},
{1.0,4.49395921,10.09783468,14.59179389,14.59579389,10.09783468,4.49395921},
{1.0,5.12583090,13.13707118,21.84615097,25.68835593,21.84615097,13.13707118,5.12583090}};

printf("\nTYPE 1.Butterworth 2.Chebyshev?(1/2):");
scanf("%d",type);
if(*type==2){
c=sqrt(pow(10,as/10.0)-1.0);
e=sqrt(pow(10,ap/10.0)-1.0);
*n=(int)(fabs(FNCCH(c/e)/FNCCH(ws/wp))+0.99999);
b=(COMPLEX*)calloc(*n+2,sizeof(COMPLEX));
if(b==NULL){
printf("\nNot enough memory to allocate!");
exit(0);
}
wc=wp;
a=pow(wc,(*n))/(e*pow(2.0,(*n)-1));
q=1/e;
x=FNSSH(q)/(*n);
for(i=0;i<*n;i++){
y=(2.0*i+1.0)*PI/(2.0*(*n));
(b+i)->real=-wc*FNSH1(x)*sin(y);
(b+i)->image=-wc*FNCH1(x)*cos(y);
}
}
else{
c=(pow(10.0,(0.1*ap))-1.0)/(pow(10.0,(0.1*as))-1.0);
*n=(int)(fabs(log(c)/log(wp/ws)/2.0)+0.99999);
b=(COMPLEX *)calloc(*n+2,sizeof(COMPLEX));
if(b==NULL){
printf("\nNot enough memory to allocate!");
exit(0);
}
wc=wp/pow(pow(10.0,0.1*ap)-1.0,1.0/(2.0*(*n)));
a=pow(wc,(double)(*n));
for(i=0;i<*n;i++){
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 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=(b1+i)->real;
for(i=0;i<=*n;i++)
printf("\nz[%2d]=%16f",i,h);
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=xs[*n-1];
for(i=0;i<=*n-1;i++)
h=pp/pow(wc,i);

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

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

yuepei

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

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*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=a;
if(n==1) {free(c);return b;}
else {
nk=m;
for(i=1;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;
}


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

yuepei

/******************************************************************************
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*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 if((*(x+i)>0)&&(*(x+i)>ym)) ym=*(x+i);
if((*(x+i)<0)&&(-*(x+i)>ym)) ym=-*(x+i);
}
for(i=0;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) 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) 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 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();
}