#include <math.h>
#include <stdio.h>
void fircad(n,type,nbands,lgrid,edge,fx,wtx,h)
int n,type,lgrid,nbands;
double h[],fx[],edge[],wtx[];
{int j,k,l,ix,nz,neg,nml,kup,nfcns,ngrid;
int iext[67],band,nodd,nzmj,nf2j,nf3j;
double pi,xt,dev,pi2,delf,fup,temp,change;
double ad[67],x[67],y[67],alpha[67];
double grid[1046],des[1046],wt[1046],deviat[11];
double eff(),wate();
void remez();
pi=4.0*atan(1.0);
pi2=2.0*pi;
if(nbands<=0) nbands=1;
if(lgrid<=0) lgrid=16;
neg=1;
if(type==1) neg=0;
nodd=n/2;
nodd=n-2*nodd;
nfcns=n/2;
if((nodd==1)&&(neg==0)) nfcns++;
grid[1]=edge[1];
delf=lgrid*nfcns;
delf=0.5/delf;
if(neg!=0)
{if(edge[1]<delf) grid[1]=delf;}
j=1;l=1;band=1;
do
{fup=edge[l+1];
do
{temp=grid[j];
des[j]=eff(temp,fx,band,type);
wt[j]=wate(temp,fx,wtx,band,type);
j++;
grid[j]=fup;
des[j-1]=eff(fup,fx,band,type);
wt[j-1]=wate(fup,fx,wtx,band,type);
band++;
l+=2;
if(band<=nbands) grid[j]=edge[l];
}while(band<=nbands);
ngrid=j-1;
if((neg==nodd)&&(grid[ngrid]>(0.5-delf))) ngrid--;
if(neg<=0)
{if(nodd!=1)
{for(j=1;j<=ngrid;j++)
{change=cos(pi*grid[j]);
des[j]/=change;
wt[j]*=change;}}}
else
{if(nodd!=1)
{for(j=1;j<ngrid;j++)
{change=sin(pi*grid[j]);
des[j]/=change;
wt[j]*=change;}}
else
{for(j=1;j<ngrid;j++)
{change=sin(pi2*grid[j]);
des[j]/=change;
wt[j]=*change;}}}
temp=(ngrid-1)/((double)nfcns);
for(j=1;j<=nfcns;j++)
{xt=j-1;
iext[j]=xt*temp+1.0;}
iext[nfcns+1]=ngrid;
nm1=nfcns-1;
nz=nfcns+1;
remez(nfcns,ngrid,&dev,ad,x,y,alpha,grid,des,wt,iext);
if(neg<=0)
{if(nodd==0)
{h[1]=0.25*alpha[nfcns];
for(j=2;j<=nm1;j++)
{nzmj=nz-j;
nf2j=nfcns+2-j;
h[j]=0.25*(alpha[nzmj]+alpha[nf2j]);}
h[nfcns]=0.5*alpha[1]+0.25*alpha[2];}
else
{for(j=1;j<=nm1;j++)
{nzmj=nz-j;
h[nfcns]=alpha[1];}}
else
{if(nodd!=0)
{h[1]=0.25*alpha[nfcns];
h[2]=0.25*alpha[nm1];
for(j=3;j<=nm1;j++)
{nzmj=nz-j;
nf3j=nfcns+3-j;
h[j]=0.25*(alpha[nzmj]-alpha[nf3j]);}
h[nfcns]=0.5*alpha[1]-0.25*alpha[3];
h[nz]=0.0;}
else
{h[1]=0.25*alpha[nfcns];
for(j=2;j<=nm1;j++)
{nzmj=nz-j;
nf2j=nfcns+2-j;
h[j]=0.25*(alpha[nzmj]-alpha[nf2j]);}
h[nfcns]=0.5*alpha[1]-0.25*alpha[2];}}
printf(“\n\n”);
for(j=1;j<=70;j++)
{printf(“*”);}
printf(“\n Finite Impulse Response (FIR)”);
printf(“\n Linear Phase Digital Filter Design”);
printf(“\n Remez Exchange Algorithm\n”);
if(type==1) printf(“\n bandpass filter\n”);
if(type==2) printf(“\n differentiator\n”);
if(type==3) printf(“\n hilbert transformer\n”);
printf(“\n * * * * * impulse response * * * * * \n”);
for(j=1;j<=nfcns;j++)
{k=n+1-j;
if(neg==0)
{printf(“ h(%2d)=%14.8lf=h(%3d)\n”,j-1,h[j],k-1);
h[k]=h[j];}
if(neg==1)
{printf(“ h(%2d)=%14.8lf=h(%3d)\n”,j-1,h[j],k-1);
h[k]=-h[j];}}
if((neg==1)&&(nodd==1))
{printf(“ h(%2d)=%14.8lf\n”,nz-1,0.0);}
for(j=0;j<n;j++)
{h[j]=h[j+1];}
h[n]=0.0;
for(k=1;k<=nbands;k+=4)
{kup=k+3;
if(kup>nbands) kup=nbands;
printf(“\n\n ”);
for(j=k;j<=kup;j++)
{printf(“%s%3d ”,”Band”,j);}
printf(“\n lower band edge”);
for(j=k;j<=kup;j++)
{printf(“%14.7lg”,edge[2*j-1]);}
printf(“\n upper band edge”);
for(j=k;j<=kup;j++)
{printf(“%14.7lf”,edge[2*j]);}
if(type!=2)
{printf(“\n desired value”);
for(j=k;j<=kup;j++)
{printf(“%14.7lf”,fx[j]);}}
if(type==2)
{printf(“\n desired slope”);
for(j=k;j<=kup;j++)
{printf(“%14.7lf”,fx[j]);}}
printf(“\n weight ”);
<< 上一页 [21] [22] [23] [24] [25] [26] [27] [28] 下一页