模态识别程序问题

维格纳

我想识别单自由度弹簧振子的频率和阻尼，数据是质量块加速度响应信号，用了《 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);

错误在哪里呀，该怎么修改，大神指点一下。

a397874795

大兄弟，你的单自由度弹簧质量块的仿真是怎么做的

suffer

个人猜测你这个估计是输入文件有问题

suffer

a397874795 发表于 2016-9-5 10:04
大兄弟，你的单自由度弹簧质量块的仿真是怎么做的

http://forum.vibunion.com/thread-10513-1-1.html

陌影

不会是文件的问题吧

think2015

找到问题了吗？