关于正交多项式拟合频响函数

fyx103

在《matlab在振动信号处理中的应用》一书中有正交多项式进行模态参数识别的程序：
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]=fun84(H1,w,1,nm);
[q,cq]=fun84(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);
该程序的错误在哪？如何改正呢？