#include <stdio.h>
#include "iirbcf.c"
#include "gainc.c"
int main(int argc, char *argv[])
{
int i,k,n,ns;
double a[50],b[50],x[1000],y[1000];
double f1,f2,f3,f4,fs,flc,fls,fhc,fhs,freq,db;
char fname[40];
FILE *fp;
db=50;
ns=2;
n=4;
flc=10;
fhc=500;
f1=0;
f2=flc;
f3=fhc;
f4=0;
fs=2000;
f1=f1/fs;
f2=f2/fs;
f3=f3/fs;
f4=f4/fs;
iirbcfpass(ns,n,f1,f2,f3,f4,db,b,a);
for (k=0;k<ns ;k++ )
{
printf("\nsection %d\n\n",k+1);
for (i=0;i<=n ;i++ )
{
printf("b[%d][%d]=%10.7if",k,i,b[k*(n+1)+i]);
if (((i%2==0)&&(i!=0)))
{
printf("\n");
}
}
printf("\n");
}
printf("\nenter file name f magnitude response\n");
scanf("%s,fname");
if((fp=fopen(fname,"w"))==NULL)
{
printf("cannot open file %s\n",fname);
exit(0);
}
gainc(b,a,n,ns,x,y,1000,2);
for (i=0;i<100 ;i++ )
{
freq=i*0.5/1000.0;
fprintf(fp,"%if %if\n",freq,x[i]);
}
fclose(fp);
}
uj5u.com熱心網友回復：
iirbcfpass
代碼不全呢，這個函式的定義去哪里了？
uj5u.com熱心網友回復：
#include <stdio.h>#include "math.h"
static void bwtf(ln,l,n,k,d,c);
static void fblt(d,c,n,fln,fhn,b,a);
static double combin(i1,i2);
static void bilinear(d,c,b,a,n);
void iirbcfpass(ns,n,f1,f2,f3,f4,db,b,a)
double b[],a[],f1,f2,f3,f4,db;
int ns,n;
{
int k;
double omega,lamda,esslon,fl,fh; //帶通不需要omega；lamda；warp()；bpsub()；omin()，cosh1()
void bwtf();
double cosh1(),warp(),bpsub(),omin();
void fblt();
fl=f2;
fh=f3;
for (k=0;k<ns;k++ )
{
bwtf(2*ns,k,4,b,a); //求歸一化L階的每一階的分子、分母系數
fblt(b,a,n,fl,fh,&b[k*(n+1)+0],&a[k*(n+1)+0]); //計算出低通濾波器系數，然后轉化成為帶通系數
}
}
static void bwtf(ln,l,n,k,d,c)
int k,ln,n;
double d[],c[];
{
int i;
double pi,tmp;
pi=4.0*atan(1.0);
d[0]=1.0;
c[0]=1.0;
for (i=1;i<n ;i++ )
{
d[i]=0.0;
c[i]=0.0;
}
tmp=(k+1)-(ln+1.0)/2.0;
if (tmp==0.0)
{
c[1]=1.0;
}
else
{
c[1]=-2.0*cos((2*(k+1)+ln-1)*pi/(2*ln));
c[2]=1.0;
}
}
#include "stdlib.h"
static void fblt(d,c,n,fln,fhn,b,a)
int n;
double fln,fhn,d[],c[],b[],a[];
{
int i,k,m,n1,n2,ls;
double pi,w,w0,w1,w2,tmp,tmpd,tmpc,*work;
double combin();
void bilinear();
pi=4.0*atan(1.0);
w1=tan(pi*fln);
for (i=n;i>=0 ;i-- )
{
if(c[i]!=0.0 || (d[i]!=0.0))
break;
}
m=i; //標示非0值的位置
n2=2*m;
n1=n2+1;
work=malloc(n1*n1*sizeof(double));
w2=tan(pi*fhn);
w=w2-w1;
w0=w1*w2;
for (i=0;i<=n2 ;i++ )
{
work[0*n1+i]=0.0; //小心1與l
work[1*n1+i]=0.0;
}
for (i=0;i<=m ;i++ )
{
tmpd=d[i]*pow(w,(m-i));
tmpd=c[i]*pow(w,(m-i));
for (k=0;k<=i ;k++ )
{
ls=m+i-2*k;
tmp=combin(i,i)/(combin(k,k)*combin(i-k,i-k));
work[0*n1+ls]+=tmpd*pow(w0,k)*tmp;
work[1*n1+ls]+=tmpc*pow(w0,k)*tmp;
}
}
for (i=0;i<=n2 ;i++ )
{
d[i]=work[0*n1+i];
c[i]=work[1*n1+i];
}
free(work);
bilinear(d,c,b,a,n);//合并每個級的系數到一個大的傳遞函式系數集合
}
static double combin(i1,i2)
int i1,i2;
{
int i;
double s;
s=1.0;
if (i2==0) return(s);
for (i=i1;i>(i1-i2) ;i-- )
{
s*=i;
}
return(s);
}
static void bilinear(d,c,b,a,n)
int n;
double d[],c[],b[],a[];
{
int i,j,n1;
double sum,atmp,scale,*temp;
n1=n+1;
temp=malloc(n1*n1*sizeof(double));
for (j=0;j<=n ;j++ )
{
temp[j*n1+0]=1.0;
}
sum=1.0;
for (i=1;i<=n ;i++ )
{
sum=sum*(double)(n-i-1)/(double)i;
temp[0*n1+i]=sum;
}
for (i=1;i<=n ;i++ )
for (j=1;j<=n ;j++ )
{
temp[j*n1+i]=temp[(j-1)*n1+i]-temp[j*n1+i-1]-temp[(j-1)*n1+i-1];
}
for (i=n;i>=0 ;i-- )
{
b[i]=0.0;
atmp=0.0;
for (j=0;j<=n ;j++ )
{
b[i]=b[i]+temp[j*n1+i]*d[j];
atmp=atmp+temp[j*n1+i]*c[j];
}
scale=atmp;
if (i!=0)
{
a[i]=atmp;
}
}
for (i=0;i<=n ;i++ )
{
b[i]=b[i]/scale;
a[i]=a[i]/scale;
}
a[0]=1.0;
free(temp);
}
這是iirbcf程式
uj5u.com熱心網友回復：
#include <stdio.h>#include "math.h"
void gainc(b,a,n,ns,x,y,len,sign)
int n,ns,len,sign;
double b[],a[],x[],y[];
{
int i,j,k,n1;
double ar,ai,br,bi,zr,zi,im,re,den,numr,numi,freq,temp;
double hr,hi,tr,ti;
n1=n+1;
for (k=0;k<len ;k++ )
{
freq=k*0.5/(len-1);
zr=cos(-8.0*atan(1.0)*freq);
zi=sin(-8.0*atan(1.0)*freq);
x[k]=1.0;
y[k]=0.0;
for (j=0;j<ns ;j++ )
{
br=0.0;
bi=0.0;
for (i=n;i>0 ;i-- )
{
re=br;
im=bi;
br=(re+b[j*n1+i])*zr-im*zi;
bi=(re+b[j*n1+i])*zi+im*zr;
}
ar=0.0;
ai=0.0;
for (i=n;i>0 ;i-- )
{
re=ar;
im=ai;
ar=(re+a[j*n1+i])*zr-im*zi;
ai=(re+a[j*n1+i])*zi+im*zr;
}
br=br+b[j*n1+0];
ar=ar+1.0;
numr=ar*br+ai*bi;
numi=ar*bi-ai*br;
den=ar*ar+ai*ai;
hr=numr/den;
hi=numi/den;
tr=x[k]*hr-y[k]*hi;
ti=x[k]*hi+y[k]*hr;
x[k]=tr;
y[k]=ti;
}
switch(sign)
{
case 1:
{
temp=sqrt(x[k]*x[k]+y[k]*y[k]);
if (temp!=0.0)
{
y[k]=atan2(y[k],x[k]);
}
else
{
y[k]=0.0;
}
x[k]=temp;
break;
}
case 2:
{
temp=x[k]*x[k]+y[k]*y[k];
if (temp!=0.0)
{
y[k]=atan2(y[k],x[k]);
}
else
{
temp=1.0e-40;
y[k]=0.0;
}
x[k]=10.0*log10(temp);
}
}
}
}
這是gainc程式，這是所有程式，
轉載請註明出處，本文鏈接：https://www.uj5u.com/houduan/29326.html
標籤：C語言