模态识别程序问题
我想识别单自由度弹簧振子的频率和阻尼，数据是质量块加速度响应信号，用了《 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',);%实测时域或频域数据
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);
%计算特征值及特征向量
= 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,:) = ;
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
%将模态频从小到大排序
=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);
错误在哪里呀，该怎么修改，大神指点一下。
大兄弟，你的单自由度弹簧质量块的仿真是怎么做的{:{13}:} 个人猜测你这个估计是输入文件有问题 a397874795 发表于 2016-9-5 10:04
大兄弟，你的单自由度弹簧质量块的仿真是怎么做的
http://forum.vibunion.com/thread-10513-1-1.html 不会是文件的问题吧 找到问题了吗？
页:
[1]