王济《Matlab在振动信号处理中的应用》书中的疑问
本帖最后由 牛小贱 于 2014-6-25 11:10 编辑各位大神，，为什么我用了王济《Matlab在振动信号处理中的应用》书中的ITD法中的程序，代码是一样的。包括例题中的输入数据文件都是一样的，无法得出结果。。显示的矩阵超出的问题。跪求大神指点，万分感谢。现在附上程序代码
%ITD法模态参数识别
%%%%%%%%%%%%%%%%%%
clear
clc
close all hidden
format long
%%%%%%%%%%%%%%%%%%
fni=input('ITD法模态参数识别-输入数据文件：','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;
%取时域数据1/2的长度
n=fix(length(b)/2); %向0圆整
%将输入时域数据赋值给列向量h
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,:)'*1i;
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
%计算自由振动响应矩阵
L=length(h);
M=L/2;%不知道为什么之后就没再出现过。。
for k=1:nm
x1(k,:)=h(k:L-(nm-k+1))';
x2(k,:)=h(k+1:L-(nm-k))';
end
%用最小二乘法求解特征发成矩阵
B=x1\x2;%B=x2*x1'*inv(x1*x1');
%计算特征值及特征向量
=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=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.6\n',F(k),D(k)*100,imag(S(k)));
end
status=fclose(fid);
报错是Attempted to access V(5,5); index out of bounds because size(V)=.
Error in ITDmethod (line 64)
U(k)=V(k,k);
数据文件为
5
2
0.1953125
out.txt
0.000000 0.000000
0.000782 -0.0081078
0.003141 0.001224
0.007117 0.027366
0.012780 0.051120
0.020231 0.127391
0.029610 0.090335
……
跪求各位指点。。。。
{:{19}:} 求求各位了，没有大神路过吗。。{:{19}:} 如果是取对角元素，可以用U=diag(V) chybeyond 发表于 2014-6-14 13:06
如果是取对角元素，可以用U=diag(V)
这个语句好像不能用呢。。大神能说详细点吗 滴滴滴滴的 发表于 2014-6-14 13:46
这个语句好像不能用呢。。大神能说详细点吗
a=;u=diag(a)%提取a的对角元素
u =
1
5
9
chybeyond 发表于 2014-6-14 13:54
懂了，感动不已，改过之后没有出问题{:{19}:}。似乎整个程序计算得到的V的大小为4，，而实际要求的要到大nm的值，也就是10 。。所以一直到后面设计nm循环的问题都会超出，，，不知道这是什么问题。。Attempted to access F2(5); index out of bounds because numel(F2)=4.
Error in ITDmethod (line 108)
if F2(k)~=F2(k+1) 滴滴滴滴的 发表于 2014-6-14 14:10
懂了，感动不已，改过之后没有出问题。似乎整个程序计算得到的V的大小为4，，而实际要求的要到大 ...
具体nm值得设置跟你需求有关，这个就不清楚了 chybeyond 发表于 2014-6-14 14:29
具体nm值得设置跟你需求有关，这个就不清楚了
已经很感谢您了，萍水相逢，就肯帮助我，想想比我的导师要更靠谱。给您一万个赞。
滴滴滴滴的 发表于 2014-6-14 14:34
已经很感谢您了，萍水相逢，就肯帮助我，想想比我的导师要更靠谱。给您一万个赞。
欢迎常来论坛学习交流，互帮互助 {:{19}:}有相似经历的大神吗。。。求帮助啊。。 求告知啊，，。，{:{10}:}已经不行了。。 滴滴滴滴的 发表于 2014-6-14 20:13
求告知啊，，。，已经不行了。。
http://forum.vibunion.com/thread-103118-1-2.html chybeyond 发表于 2014-6-14 20:14
http://forum.vibunion.com/thread-103118-1-2.html
为啥我所在的用户组，我法查看或下载附件 滴滴滴滴的 发表于 2014-6-15 08:09
为啥我所在的用户组，我法查看或下载附件
〖新手必读〗之如何获取积分，提高权限（新）http://forum.vibunion.com/thread-132101-1-1.html
页:
[1]
2