关于《matlab在振动信号处理中的应用》中程序错误

fyx103

书中的正交多项式法模态参数识别部分程序如下：
clear
format long
global mn;
%fni=input('输入数据文件名称：','s');
fi=fopen('t1.txt','r');       %打开输入数据文件
mn=fscanf(fi,'%d',1);    %输入模态阶数
df=fscanf(fi,'%f',1);    %输入频率间隔
fo=fscanf(fi,'%s',1);    %输出文件名称
h=fscanf(fi,'%f',[2,inf]);   %输入实测频响函数实部和虚部
status=fclose(fi);       %关闭文件
nm=mn*2;                %定义频响函数有理分式的阶数
n=length(h(1,:));       %频响函数长度
f=0:df:(n-1)*df;        %定义频率向量
w=[-f(n:-1:1) f(1:n)]/max(f);
H1=[h(1,n:-1:1)-i*h(2,n:-1:1) h(1,1:n)+i*h(2,1:n)];%建立扩展的实测频响函数向量
w=w.';
H1=H1.';
[p,cp]=fun85(H1,w,1,nm);
[q,cq]=fun85(H1,w,2,nm);
for k=1:nm
    Q(:,k)=(H1.*q(:,k));
end
P=p;
T=-real(P'*Q);
g=real(P'*(H1.*q(:,nm+1)));
D=-inv(eye(nm)-T'*T)*T'*g;
C=g-T*D;
D(nm+1)=1;
A=cp*C;
B=cq*D;
A=A(nm+1:-1:1).';
B=B(nm+1:-1:1).';
[R,U,K]=residue(A,B);
F1=abs(U)*max(f);
运行后模态频率得不到正确值，请问有研究过这本书的有正确的程序吗？
以下是fun85：
if ig==1
   W=ones(size(w));
else
   W=(abs(H)).^2;
end
R0=zeros(size(w));
R1=ones(size(w))./sqrt(sum(W));
R=[R0,R1];
cf=zeros(m+1,m+2);
cf(1,2)=1/sqrt(sum(W));
for k=1:m
    V=sum(w.*R(:,k+1).*R(:,k).*W);
    S=w.*R(:,k+1)-V*R(:,k);
    D=sqrt(sum((S.^2).*W));
    R=[R,S/D];
    cf(:,k+2)=-V*cf(:,k);
    cf(2:k+1,k+2)=cf(2:k+1,k+2)+cf(1:k,k+1);
    cf(:,k+2)=cf(:,k+2)/D;
end
R=R(:,2:m+2);
cf=cf(:,2:m+2);
for k=1:m+1
    P(:,k)=R(:,k)*i^(k-1);
end