求教:哪位高手帮我看看修改一下程序?

ddmm

下面是在频域下利用幅度平方误差最小法设计的IIR数字滤波器的程序,但在运行错误,我修改了两天,还是没改对,麻烦看看,指导一下.
function [A,para]=IIRdirFreq(Hd,k)
u=0.001;N=5000;E=zeros(1,N+1);
M=length(Hd);
w=0:pi/(M+1):pi;
z1=exp(-j.*w);
z2=z1.*z1;
fai=zeros(1,4*k);
difffai=zeros(1,4*k);
E(2)=-0.5;n=3;
while(abs(E(n-1)-E(n-2))>0.001)&(n<=N)
P=zeros(1,M);
a1=zeros(1,M);
b1=zeros(1,M);
     for n1=1:k
        m1=4*(n1-1);
        a1=1+fai(m1+1)*z1+fai(m1+2);
        b1=1+fai(m1+3)*z1+fai(m1+4);
    end
    P=P*(a1/b1);
    P=abs( P);
     A0=abs(P*Hd')/(P*P');
    E(n)=sum(A0*P-abs(Hd).^2);
    for n1=1:k
        m1=4*(n1-1);
        c=2*A0*(A0.*P-abs(Hd));
        diff_fai(m1+1)=c*(P.*real(z1./(1+fai(m1+1)*z1+fai(m1+2)*z2)))';
        diff_fai(m1+2)=c*(P.*real(z2./(1+fai(m1+1)*z1+fai(m1+2)*z2)))';
        diff_fai(m1+3)=c*(-P.*real(z1./(1+fai(m1+3)*z1+fai(m1+4)*z2)))';
        diff_fai(m1+4)=c*(-P.*real(z1./(1+fai(m1+3)*z1+fai(m1+4)*z2)))';
    end
    fai=fai-u.*diff_fai;
    n=n+1;
end
para=fai;
A=A0;

ericlin

低手帮你顶一下.