mremez1.c

数字信号处理_胡广书(笫二版)光盘 包含该书的一些C

extern float ad[66],x[66],y[66],grid[1045],des[1045],wt[1045],alpha[66];
extern int iext[66];
extern float pi2,dev;
extern int nfcns,ngrid,niter;
/*-------------------------------------------------------------------*/
        float d(int k,int n,int m)
        {float z,q;
         int l,j;
         z=1.;
         q=x[k];
         for(l=1;l<=m;l++)
            {for(j=l;j<=n;j+=m)
                 if((j-k)!=0)
                     z*=2.*(q-x[j]);
             }
         z=1./z;
         return(z);
         }
/*------------------------------------------------------------------*/
        float gee(int k,int n)
        {float p,xf,d1,c;
         int j;
         p=0.0;
         xf=grid[k];
         xf=cos(pi2*xf);
         d1=0.;
         for(j=1;j<=n;j++)
            {c=xf-x[j];
             c=ad[j]/c;
             d1+=c;
             p+=c*y[j];
             }
         p/=d1;
         return(p);
         }
/*---------------------------------------------------------------------
        Subroutine mremez1.c
                                       in chapter 8
---------------------------------------------------------------------*/
void mremez1()
{
        float a[66],p[66],q[66];
        float devl,dnum,dden,dtemp,fsh,gtemp,aa,bb,ft,xt,xe;
        float ynz,comp,y1,err,delf,xt1;
        int itrmax,nz,nzz,jxt,jet,k,l,nu,jchnge;
        int nut,klow,k1,knz,kup,luck,nut1;
        int j,nzzmj,nzmj,kn,nm1,kkk,cn,jm1,nf1j;
        int flag;
        comp=0.0;
        itrmax=25;
        devl=-1.;
        nz=nfcns+1;
        nzz=nfcns+2;
        niter=0;
        while(1)
          {
            iext[nzz]=ngrid+1;
            niter++;
            if(niter>itrmax)break;
            for(j=1;j<=nz;j++)
               {jxt=iext[j];
                dtemp=grid[jxt];
                dtemp=cos(dtemp*pi2);
                x[j]=dtemp;
                }
            jet=(nfcns-1)/15+1;
            for (j=1;j<=nz;j++)
                 ad[j]=d(j,nz,jet);
            dnum=0.0;
            dden=0.0;
            k=1;
            for(j=1;j<=nz;j++)
               {
                    l=iext[j];
                    dtemp=ad[j]*des[l];
                    dnum+=dtemp;
                    dtemp=k*ad[j]/wt[l];
                    dden+=dtemp;
                    k=-k;
                }
            dev=dnum/dden;
            nu=1;
            if(dev>0.)
               nu=-1;
            dev*=-nu;
            k=nu;
            for(j=1;j<=nz;j++)
               {
                    l=iext[j];
                    dtemp=k*dev/wt[l];
                    y[j]=des[l]+dtemp;
                    k=-k;
                }
            if(dev<=devl)break;
            devl=dev;
            jchnge=0;
            k1=iext[1];
            knz=iext[nz];
            klow=0;
            nut=-nu;
            j=1;
            while(1)
            {
                while(1)
                    {
                        if(j==nzz)
                           ynz=comp;
                        if(j>=nzz)break;
                        kup=iext[j+1];
                        l=iext[j]+1;
                        nut=-nut;
                        if(j==2)
                           y1=comp;
                        comp=dev;
                        if(l>=kup)break;
                        err=gee(l,nz);
                        err=(err-des[l])*wt[l];
                        dtemp=nut*err-comp;
                        if(dtemp<=0.)break;
                        comp=nut*err;
                        while(1)
                            {
                                l++;
                                if(l>=kup)break;
                                err=gee(l,nz);
                                err=(err-des[l])*wt[l];
                                dtemp=nut*err-comp;
                                if(dtemp<=0.)break;
                                comp=nut*err;
                            }
                        iext[j]=l-1;
                        j++;
                        klow=l-1;
                        jchnge++;
                    }
                if(j<nzz)
                {
                    l--;
                    flag=0;
                    while(1)
                        {
                            l--;
                            if(l<=klow)break;
                            err=gee(l,nz);
                            err=(err-des[l])*wt[l];
                            dtemp=nut*err-comp;
                            if(dtemp>0.)break;
                            if(jchnge>0.)
                                {
                                    klow=iext[j];
                                    j++;
                                    flag=1;
                                    break;
                                }
                        }
                    if(flag==1)continue;
                    flag=0;
                    if(l>klow)
                        {
                           comp=nut*err;                 
                            while(1)
                                {
                                    l--;
                                    if(l<=klow)break;
                                    err=gee(l,nz);
                                    err=(err-des[l])*wt[l];
                                    dtemp=nut*err-comp;
                                    if(dtemp<=0.)break;
                                    comp=nut*err;
                                }
                            klow=iext[j];
                            iext[j]=l+1;
                            j++;
                            jchnge++;
                            flag=1;
                        }
                    if(flag==1)continue;
                    l=iext[j]+1;
                    flag=0;
                    if(jchnge>0.)
                       {
                        iext[j]=l-1;
                        j++;
                        klow=l-1;
                        jchnge++;
                        flag=1;
                       }
                    if(flag!=1)
                        {
                         while(1)
                            {
                                l++;
                                if(l>=kup)break;
                                err=gee(l,nz);
                                err=(err-des[l])*wt[l];
                                dtemp=nut*err-comp;
                                if(dtemp>0.)
                                    {  
                                        comp=nut*err;
                                            while(1)
                                            {
                                                l++;
                                                if(l>=kup)break;
                                                err=gee(l,nz);
                                                err=(err-des[l])*wt[l];
                                                dtemp=nut*err-comp;
                                                if(dtemp<=0.)break;
                                                comp=nut*err;
                                            }
                                        iext[j]=l-1;
                                        jchnge++;                     
                                        break;
                                    }    
                            }
                            klow=iext[j];
                            j++;
                        }
                    continue;
                }
              if(j>nzz)               
                     {
                        flag=0;
                        if(luck>9)
                            {
                                kn=iext[nzz];
                                for(j=1;j<=nfcns;j++)
                                    iext[j]=iext[j+1];
                                iext[nz]=kn;
                                flag=2;
                                break;
                            }
                        if(comp>y1)
                           y1=comp;
                        k1=iext[nzz];
                      }
                else
                    {   
                        flag=0;
                if(k1>iext[1])
                   k1=iext[1];
                if(knz<iext[nz])
                   knz=iext[nz];
                nut1=nut;
                nut=-nu;
                l=0;
                kup=k1;
                comp=ynz*1.00001;
                luck=1;
                while(1)
                {
                    l++;
                    if(l>=kup)break;
                    err=gee(l,nz);
                    err=(err-des[l])*wt[l];
                    dtemp=nut*err-comp;
                    flag=0;
                    if(dtemp>0.)
                        {
                            comp=nut*err;
                            j=nzz;
                            while(1)           
                                {
                                    l++;
                                    if(l>=kup)break;
                                    err=gee(l,nz);
                                    err=(err-des[l])*wt[l];
                                    dtemp=nut*err-comp;
                                    if(dtemp<=0.)break;
                                    comp=nut*err;
                                }
                            iext[j]=l-1;
                            j++;
                            klow=l-1;
                            jchnge++;
                            flag=1;
                            break;
                        }
                }
                        if(flag==1)continue;
                        luck=6;
                    }
            l=ngrid+1;
            klow=knz;
            nut=-nut1;
            comp=y1*1.00001;
            while(1)
            {
                l--;
                if(l<=klow)break;
                err=gee(l,nz);
                err=(err-des[l])*wt[l];
                dtemp=nut*err-comp;
                flag=0;
                if(dtemp>0.)
                    {
                        j=nzz;
                        comp=nut*err;
                        luck+=10;
                        while(1)
                            {
                                l--;
                                if(l<=klow)break;
                                err=gee(l,nz);
                                err=(err-des[l])*wt[l];
                                dtemp=nut*err-comp;
                                if(dtemp<=0.)break;
                                comp=nut*err;
                            }
                        klow=iext[j];
                        iext[j]=l+1;
                        j++;
                        jchnge++;
                        flag=1;
                        break;
                    }
            }
                if(flag!=1)break;
            }
            if(flag==2)continue;
            if(luck!=6)
                {
                    for(j=1;j<=nfcns;j++)
                       {
                        nzzmj=nzz-j;
                        nzmj=nz-j;
                        iext[nzzmj]=iext[nzmj];
                        }
                    iext[1]=k1;
                    continue;
                }
         if(jchnge>0)continue;   
        }
        nm1=nfcns-1;
        fsh=1.e-06;
        gtemp=grid[1];
        x[nzz]=-2.;
        cn=2*nfcns-1;
        delf=1./cn;
        l=1;
        kkk=0;
        if(grid[1]<=.01&&grid[ngrid]>0.49)
           kkk=1;
        if(nfcns<=3)
           kkk=1;
        if(kkk!=1)
            {
                dtemp=cos(pi2*grid[1]);
                dnum=cos(pi2*grid[ngrid]);
                aa=2./(dtemp-dnum);
                bb=-(dtemp+dnum)/(dtemp-dnum);
            }
        for(j=1;j<=nfcns;j++)
           {
                ft=j-1;
                ft*=delf;
                xt=cos(pi2*ft);
                if(kkk!=1)
                    {
                        xt=(xt-bb)/aa;
                        xt1=sqrt(1.-xt*xt);
                        ft=atan2(xt1,xt)/pi2;
                    }
               while(1)
                    {
                        xe=x[l];
                        if(xt>xe)break;
                        if((xe-xt)<fsh)break;
                        l++;
                    }
               if((xe-xt)<fsh||(xt-xe)<fsh)
                   a[j]=y[l]; 
               if(xt>xe)
                   {
                       grid[1]=ft;
                       a[j]=gee(1,nz);
                   }
               if(l>1)l--;
           }
        grid[1]=gtemp;
        dden=pi2/cn;
        for(j=1;j<=nfcns;j++)
           {
                dtemp=0.;
                dnum=j-1;
                dnum*=dden;
                if(nm1>=1)
                    for(k=1;k<=nm1;k++)
                        dtemp+=a[k+1]*cos(dnum*k);
                dtemp=dtemp*2.+a[1];
                alpha[j]=dtemp;
            }
        for(j=2;j<=nfcns;j++)
            alpha[j]*=2./cn;
        alpha[1]/=cn;
        if(kkk!=1)
        {
            p[1]=2.*alpha[nfcns]*bb+alpha[nm1];
            p[2]=2.*alpha[nfcns]*aa;
            q[1]=alpha[nfcns-2]-alpha[nfcns];
            for(j=2;j<=nm1;j++)
               {
                    if(j>=nm1)
                      {
                       aa*=.5;
                       bb*=.5;
                       }
                    p[j+1]=0.;
                    for(k=1;k<=j;k++)
                       {
                        a[k]=p[k];
                        p[k]=2.*bb*a[k];
                        }
                    p[2]+=a[1]*2.*aa;
                    jm1=j-1;
                    for(k=1;k<=jm1;k++)
                        p[k]+=q[k]+aa*a[k+1];
                    for(k=3;k<=j+1;k++)
                        p[k]+=aa*a[k-1];
                    if(j==nm1) continue;
                    for(k=1;k<=j;k++)
                        q[k]=-a[k];
                    nf1j=nfcns-j-1;
                    q[1]+=alpha[nf1j];
                }
            for(j=1;j<=nfcns;j++)
                alpha[j]=p[j];
         }
       if(nfcns>3)
           return;
        alpha[nfcns+1]=0.0;
        alpha[nfcns+2]=0.0;
        return;
}