matlab在振动信号处理中的应用,关于《matlab在振动信号处理中的应用》

本帖最后由 fyx10319 于 2019-4-11 20:23 编辑

本人正在学习使用正交多项式法拟合频响函数的matlab实现，发现王济老师的《matlab在振动信号处理中的应用》里面的程序运行得不到正确结果，请问该程序是否有错误？能否指导一下如何改正？谢谢！

clear

format long

global mn;

fni=input('输入数据文件名称：','s');

fi=fopen(fni,'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);

以上是前部分程序，F1得不到正确解。

以下是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