iird.c

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

/*********************************************************************
* 双线性变换法设计 IIR 数字滤波器程序
/*********************************************************************/
#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, imag;
} COMPLEX;

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

#define    PI	      (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

/***********************************************************************/
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;

 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("\nNot 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[i]);

 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],i,ptr->b[i]);

 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.imag=0;
    hwdb2.real=0; hwdb2.imag=0;
    for(i=0;i<=nz;i++)  {
       hwdb1.real += ptr->a[i]*cos(i*jw);
       hwdb2.real += ptr->b[i]*cos(i*jw);
       hwdb1.imag += ptr->a[i]*sin(i*jw);
       hwdb2.imag += ptr->b[i]*sin(i*jw);
    }
    amp1 = (pow(hwdb1.real,2)+pow(hwdb1.imag,2));
    amp2 = (pow(hwdb2.real,2)+pow(hwdb2.imag,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);
}

/***********************************************************************/
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,Ap: Passband 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,Ap: Passband 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<fr)||(*ap>*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<fr)||(*ap>*ar)||(fs<2*fp));
 }
 *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;
}
/***********************************************************************/
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,Ap: Passband 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 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)||(fp1<fr1)||(fp2>fr2)||(*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)||(fp1<fr1)||(fp2>fr2)||(*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) < fabs(pws2-pwp2)) {
    *wp=pwp1;
    *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; *(f2+2)=1.0;
 *nf = 2;
}

/************************************************************************
bcg - Chebyshev 和 Buttterworth 型模拟原型传输函数生成子程序
      即程序得到系统函数H(s).
      输出格式为: Ha(s)=1/(b0+b1s+b2s^2...+bns^n). 分母为升幂排列.
输入参数:
      double ap,as : 通带和阻带衰耗值,单位为 dB;
      double wp,ws : 通带和阻带特征频率,单位为 rad/s;
输出参数:
      int *n       : 传输函数 H(s)  的阶数;
      double *h    : 传输函数 Ha(s) 的系数指针,系数为升幂排列;
      int *type    : 滤波器类型  1.BW 型  2.CB 型
*************************************************************************/
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 *b1,*b2,*b;

 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)->imag = -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)->imag=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->imag = -(b->imag);