levy.m

对于振动信号处理时域信号

%有理分式多项式模态参数识别
clear
clc
close all hidden
format long
%%%%%%%%%%%%%%%%%%
fni=input('有理分式多项式法模态参数识别-输入数据文件名：','s');
fid=fopen(fni,'r');
mn=fscanf(fid,'%d',1); %模态阶数
df=fscanf(fid,'%f',1);%频率间隔
fno=fscanf(fid,'%s',1);%输出数据文件名
b=fscanf(fid,'%f',[2,inf]);%实测频响实部和虚部数据
status=fclose(fid)
%定义幂多项式的阶次
nm=mn*2;
%取频响函数的长度
n=length(b(1,:));
%建立离散频率向量
f=0:df:(n-1)*df;
%建立离散圆频率向量
w=2*pi*f;
%建立归一化离散圆频率向量
wi=w/max(w);
%建立实测频响函数复数向量
H=b(1,:)+i*b(2,:);
%计算拟合频响函数的分子和分母系数向量
[A ,B]=invfreqs(H,wi,nm,nm,[],100);
%幂多项时求根（零点）
P=roots(B);
%计算模态频率向量
F1=abs(P)*max(w)/(2*pi);
%计算阻尼比向量
D1=-real(P)./(abs(P));
%计算振型系数向量

for k=1:nm
  if k==1
     p(1:nm-1)=P(2:nm);
   else
    p(1:k-1)=P(1:k-1);
    p(k:nm-1)=P(k+1:nm);
  end
y=poly(p);
S1 (k)=polyval(A,P(k))/polyval(y,P(k));
end

%计算拟合的频响函数
H1=freqs(A,B,wi);
%绘制频响函数实部拟合曲线图
figure(2)
nn=1:n;
subplot(2,1,1);
plot(f,real(H(nn)),':',f,real(H1(nn)),'r');
xlabel('频率(HZ)');
ylabel('实部');
legend('实部','拟合');
grid on
%绘制频响函数虚部拟合曲线图
subplot(2,1,2);
plot(f,b(2,:),':',f,image(H1(nn)),'r');
xlabel('频率(HZ)');
ylabel('实部');
legend('实部','拟合');
grid on
%将模态频率从小到大排序
[F2,I]=sort(F1);
 %剔除方程解中的非模态项（非共轭根）和共轭项（重复项）
m=0;
for k=1:nm-1
  if F2(k)~=F2(k+1)
  continue;
 end
 m=m+1;
F(m)=F1(I);%模态频率
D(m)=D1(I);%阻尼比
S(m)=S1(I);%振型系数
end
%打开文件输出模态参数数据
fid=fopen(fno,'w');
fprintf(fid,'频率(HZ)  阻尼比（%%） 振型系数\n');
for k=1:m
  fprint(fid,'10.4f %10.4f %10.6f\n',F(k),D(k)*100.0,image(S(k)));
end

status=fclose(fid);