声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 1881|回复: 2

[编程技巧] prony 法对一信号进行模态参数辨识

[复制链接]
发表于 2013-4-16 09:36 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
王济著的《MATLAB在振动信号处理中的应用》这本书,最近看了看。第九章例9.5用prony法对系统进行模态辨识,想问一下,如果把数据文件改成一个正弦信号,比如y=sin5t+sin6t,程序应该怎么改呢,不解中。
原程序如下:
clear                           %消除内存中所有变量和函数
format long
%%%%%%%%%%%%%%%%%%%%%%%%%
%声明为全局变量
global mn;
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:mn
    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,imag(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;
    l=I(k);
    F(m)=F1(l);  %模态频率
    D(m)=D1(l);  %阻尼比
    S(m)=S1(l);  %振型系数
end
%打开文件输出模态参数数据
fid=fopen(fno,'w');
fprintf(fid,'频率(Hz) 阻尼比(%%) 振型系数\n');
for k=1:mn
   fprintf(fid,'%10.4f %10.4f %10.6f\n',F(k),D(k)*100.0,S(k));
end
status=fclose(fid);

回复
分享到:

使用道具 举报

发表于 2014-9-9 16:20 | 显示全部楼层
这个是有理分式多项式法(Levy法)的程序啊
发表于 2015-4-8 20:50 | 显示全部楼层
学习怎样弄到积分
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-25 12:32 , Processed in 0.079296 second(s), 17 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表