声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 2083|回复: 5

[综合讨论] 模态识别程序问题

[复制链接]
发表于 2016-7-20 16:25 | 显示全部楼层 |阅读模式

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

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

x
我想识别单自由度弹簧振子的频率和阻尼,数据是质量块加速度响应信号,用了《 MATLAB在振动信号处理中的应用》里的STD法程序,运行后输入文件数据名,然后就报错:
STD法模态参数识别-输入数据文件名:xfile.txt
Undefined function or variable 'x1'.

Error in std (line 60)
B(:,nm)=x1\x2(:,nm);%B(:,nm)=inv(x1'*x1)*x1'*x2(:,nm);



下面是原程序代码:
%STD法模态参数识别
%%%%%%%%%%%%%%%%%%
clear
clc
close all hidden
format long
%%%%%%%%%%%%%%%%%%
fni=input('STD法模态参数识别-输入数据文件名:','s');
fid=fopen(fni,'r');
mn = fscanf(fid,'%d',1);      %模态阶数
%定义输入实测数据类型:
%ig=1 时域数据如冲击响应、自由振动、互相关函数、随机减量法处理结果
%ig=2 频域数据如频响函数实部虚部数据
ig = fscanf(fid,'%d',1);      
%ig=1时f为采样频率sf,ig=2时f为频率间隔df
f = fscanf(fid,'%f',1);      
fno = fscanf(fid,'%s',1);     %输出数据文件名
b = fscanf(fid,'%f',[ig,inf]);%实测时域或频域数据
status=fclose(fid);

%建立特征方程矩阵的阶数(为模态阶数的2倍)
nm=2*mn;
if ig==1 %实测时域数据
  sf=f;
  n=fix(length(b)/2);
  h=b(1,1:2*n)';
%计算时间间隔
  dt=1/sf;
%建立离散时间向量
  t=0:dt:(2*n-1)*dt;
else     %实测频域数据
%取实测频响函数的长度
  df=f;
  n=length(b(1,:));
%建立离散频率向量
  f=0:df:(n-1)*df;
%建立对应正负频率的实测频响函数向量
  H=b(1,:)'+b(2,:)'*i;
  H(n+1) = real(H(n));
  H(n+2:2*n) = conj(H(n:-1:2));
%频响函数经IFFT并取实部变换成脉冲响应函数,
  h = real(ifft(H));
%建立离散时间向量
  t = linspace(0,1/df,2*n);
%计算时间间隔
  dt = t(2)-t(1);
end

%计算自由振动响应矩阵
M=length(h)-nm;
for k = 1:M
  x1(k,:) = h(k:k+nm-1)';
  x2(k,:) = h(k+1:k+nm)';
end

%计算Hessenberg矩阵
B=zeros(nm,nm);
B(2:nm,1:nm-1)=eye(nm-1,nm-1);
%用最小二乘法求解待定系数列向量
B(:,nm)=x1\x2(:,nm);%B(:,nm)=inv(x1'*x1)*x1'*x2(:,nm);
%计算特征值及特征向量
[A,V] = eig(B);
%变换特征值对角阵为一向量   
for k = 1:nm  
  U(k)=V(k,k);  
end
%计算模态频率
F1 = abs(log(U'))./(2*pi*dt);
%计算阻尼比
D1 = sqrt(1./(((imag(log(U'))./real(log(U'))).^2)+1));
%计算振型系数向量(特征向量)
l=1;
for k = 1:nm
  if abs(real(U(k)))<= 1 & abs(imag(U(k)))<= 1
    V0(l) = U(k);
    l = l + 1;
  end
end
for k = 0:(2*n - 1)  
  Va(k+1,:) = [conj(V0).^k ];  
end
S1 = (inv(conj(Va')*Va)*conj(Va')*h);
%计算生成的脉冲响应函数
h1=real(Va*S1);
%绘制脉冲响应函数拟合曲线图
figure(1)
plot(t,h,':',t,h1);
xlabel('时间 (s)');  
ylabel('幅值');
legend('实测','拟合');
grid on;                        
if ig>1
%计算生成的频响函数
  H1 = fft(Va*S1);
%绘制频响函数实部拟合曲线图
  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;                        
end
%将模态频从小到大排序
[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:m
  fprintf(fid,'%10.4f  %10.4f  %10.6f\n',F(k),D(k)*100.0,imag(S(k)));
end
status=fclose(fid);

错误在哪里呀,该怎么修改,大神指点一下。
回复
分享到:

使用道具 举报

发表于 2016-9-5 10:04 | 显示全部楼层
大兄弟,你的单自由度弹簧质量块的仿真是怎么做的

点评

http://forum.vibunion.com/thread-10513-1-1.html  详情 回复 发表于 2016-9-5 13:26
发表于 2016-9-5 13:24 | 显示全部楼层
个人猜测你这个估计是输入文件有问题
发表于 2016-9-5 13:26 | 显示全部楼层
a397874795 发表于 2016-9-5 10:04
大兄弟,你的单自由度弹簧质量块的仿真是怎么做的

http://forum.vibunion.com/thread-10513-1-1.html
发表于 2016-9-6 13:13 | 显示全部楼层
不会是文件的问题吧
发表于 2016-9-7 08:54 | 显示全部楼层
找到问题了吗?
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

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

GMT+8, 2024-11-9 10:13 , Processed in 0.092233 second(s), 26 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

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