声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 3963|回复: 20

[编程技巧] 王济《Matlab在振动信号处理中的应用》书中的疑问

[复制链接]
发表于 2014-6-14 12:49 | 显示全部楼层 |阅读模式

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

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

x
本帖最后由 牛小贱 于 2014-6-25 11:10 编辑

各位大神,,为什么我用了王济《Matlab在振动信号处理中的应用》书中的ITD法中的程序,代码是一样的。包括例题中的输入数据文件都是一样的,无法得出结果。。显示的矩阵超出的问题。跪求大神指点,万分感谢。现在附上程序代码
  1. %ITD法模态参数识别
  2. %%%%%%%%%%%%%%%%%%
  3. clear
  4. clc
  5. close all hidden
  6. format long
  7. %%%%%%%%%%%%%%%%%%
  8. fni=input('ITD法模态参数识别-输入数据文件:','s');
  9. fid=fopen(fni,'r');
  10. mn=fscanf(fid,'%d',1); %模态阶数
  11. %定义输入实测数据类型
  12. %ig=1时域数据如冲击响应、自由振动、互相关函数、随机减量法处理结果
  13. %ig=2频域数据如频响函数实部虚部数据
  14. ig=fscanf(fid,'%d',1);
  15. %ig=1时f为采样频率(每秒采样次数)sf,ig=2时f为频率间隔df
  16. f=fscanf(fid,'%f',1);
  17. fno=fscanf(fid,'%s',1);  %输出数据文件名
  18. b=fscanf(fid,'%f',[ig,inf]); %实测时域或频域数据
  19. status=fclose(fid);
  20. %建立特征发成矩阵的阶数(为模态阶数的2倍)
  21. nm=2*mn;
  22. %组织识别计算所用的时域数据参数
  23. if ig==1 %若数据类型为实测时域数据时
  24. %取采样频率
  25. sf=f;
  26. %取时域数据1/2的长度
  27. n=fix(length(b)/2); %向0圆整
  28. %将输入时域数据赋值给列向量h
  29. h=b(1,1:2*n)';
  30. %计算时间间隔
  31. dt=1/sf;
  32. %建立离散时间向量
  33. t=0:dt:(2*n-1)*dt;
  34. else %若为实测频域数据
  35. %取频率间隔
  36. df=f;
  37. %取实测频响函数的长度
  38. n=length(b(1,:));
  39. f=0:df:(n-1)*df;
  40. %建立对应正负频率的实测频响函数向量
  41. H=b(1,:)'+b(2,:)'*1i;
  42. H(n+1)=real(H(n));
  43. H(n+2:2*n)=conj(H(n:-1:2));
  44. %频响函数经IFFT并取实部变换成脉冲响应函数
  45. h=real(ifft(H));
  46. %建立离散时间向量
  47. t=linspace(0,1/df,2*n);
  48. %计算时间间隔
  49. dt=t(2)-t(1);
  50. end
  51. %计算自由振动响应矩阵
  52. L=length(h);
  53. M=L/2;%不知道为什么之后就没再出现过。。
  54. for k=1:nm
  55.     x1(k,:)=h(k:L-(nm-k+1))';
  56.     x2(k,:)=h(k+1:L-(nm-k))';
  57. end
  58. %用最小二乘法求解特征发成矩阵
  59. B=x1\x2;%B=x2*x1'*inv(x1*x1');
  60. %计算特征值及特征向量
  61. [A,V]=eig(B);
  62. %变换特征值对角阵为一向量
  63. for k=1:nm
  64.    U(k)=V(k,k);(这里会报错,说是矩阵没有这么大)
  65. end
  66. %计算模态频率向量
  67. F1=abs(log(U'))./(2*pi*dt);
  68. %计算阻尼比向量
  69. D1=sqrt(1./(((imag(log(U'))./real(log(U'))).^2)+1));
  70. %计算阵型系数向量
  71. l=1;
  72. for k=0:(2*n-1)
  73.     Va(k+1,:)=[conj(U).^k];
  74. end,
  75. S1=(inv(conj(Va')*Va)*conj(Va')*h);
  76. %计算生成的脉冲响应函数
  77. h1=real(Va*S1);
  78. %绘制脉冲响应函数拟合曲线图
  79. figure(1)
  80. plot(t,h,':',t,h1);
  81. xlabel('时间(s)');
  82. ylabel('幅值');
  83. legend('实测','拟合');
  84. grid on;
  85. if ig>1
  86. %计算生成的频响函数
  87. H1=fft(Va*S1);
  88. %绘制频响函数实部拟合曲线图
  89. figure(2)
  90. nn=1:n;
  91. subplot(2,1,1);%将实、虚部图线画在同一张图中
  92. plot(f,real(H(nn)),':',f,real(H1(nn)),'r');
  93. xlabel('频率(Hz)');
  94. ylabel('实部');
  95. legend('实测','拟合');
  96. grid on;
  97. %绘制频响函数虚部拟合曲线图
  98. subplot(2,1,2);
  99. plot(f,b(2,:),':',f,imag(H1(nn)),'r');
  100. xlabel('频率(Hz)');
  101. ylabel('虚部');
  102. legend('实测','拟合');
  103. grid on;
  104. end
  105. %将自振频率从小到大排序
  106. [F2,I]=sort(F1);%带有标记信息
  107. %剔除方法解中的非模态项(非共轭根)和共轭项(重复项)
  108. m=0;
  109. for k=1:nm-1
  110.     if F2(k)~=F2(k+1)
  111.         continue;
  112.     end
  113.     m=m+1;
  114.     l=I(k);
  115.     F(m)=F1(l);%自振频率
  116.     D(m)=D1(l);%阻尼比
  117.     S(m)=S1(l);%阵型系数
  118. end
  119. %打开文件输出识别的模态参数数据
  120. fid=fopen(fno,'w');
  121. fprintf(fid,'  频率(HZ)   阻尼比(%%)   阵型系数\n')
  122. for k=1:m
  123.     fprintf(fid,'%10.4f %10.4f %10.6\n',F(k),D(k)*100,imag(S(k)));
  124. end
  125. status=fclose(fid);
复制代码

报错是Attempted to access V(5,5); index out of bounds because size(V)=[4,4].

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
……

跪求各位指点。。。。


回复
分享到:

使用道具 举报

 楼主| 发表于 2014-6-14 12:50 | 显示全部楼层

点评

反对: 4.0
反对: 4
  发表于 2014-6-29 13:23
 楼主| 发表于 2014-6-14 12:52 | 显示全部楼层
求求各位了,没有大神路过吗。。
发表于 2014-6-14 13:06 | 显示全部楼层
如果是取对角元素,可以用U=diag(V)

评分

1

查看全部评分

 楼主| 发表于 2014-6-14 13:46 | 显示全部楼层
chybeyond 发表于 2014-6-14 13:06
如果是取对角元素,可以用U=diag(V)

这个语句好像不能用呢。。大神能说详细点吗
发表于 2014-6-14 13:54 | 显示全部楼层
滴滴滴滴的 发表于 2014-6-14 13:46
这个语句好像不能用呢。。大神能说详细点吗
  1. a=[1 2 3;4 5 6;7 8 9];u=diag(a)%提取a的对角元素

  2. u =

  3.      1
  4.      5
  5.      9
复制代码


评分

1

查看全部评分

 楼主| 发表于 2014-6-14 14:10 | 显示全部楼层

懂了,感动不已,改过之后没有出问题。似乎整个程序计算得到的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:29 | 显示全部楼层
滴滴滴滴的 发表于 2014-6-14 14:10
懂了,感动不已,改过之后没有出问题。似乎整个程序计算得到的V的大小为4,,而实际要求的要到大 ...

具体nm值得设置跟你需求有关,这个就不清楚了
 楼主| 发表于 2014-6-14 14:34 | 显示全部楼层
chybeyond 发表于 2014-6-14 14:29
具体nm值得设置跟你需求有关,这个就不清楚了

已经很感谢您了,萍水相逢,就肯帮助我,想想比我的导师要更靠谱。给您一万个赞。
发表于 2014-6-14 14:36 | 显示全部楼层
滴滴滴滴的 发表于 2014-6-14 14:34
已经很感谢您了,萍水相逢,就肯帮助我,想想比我的导师要更靠谱。给您一万个赞。

欢迎常来论坛学习交流,互帮互助
 楼主| 发表于 2014-6-14 18:09 | 显示全部楼层
有相似经历的大神吗。。。求帮助啊。。
 楼主| 发表于 2014-6-14 20:13 | 显示全部楼层
求告知啊,,。,已经不行了。。
发表于 2014-6-14 20:14 | 显示全部楼层
 楼主| 发表于 2014-6-15 08:09 | 显示全部楼层
chybeyond 发表于 2014-6-14 20:14
http://forum.vibunion.com/thread-103118-1-2.html

为啥我所在的用户组,我法查看或下载附件
发表于 2014-6-15 08:14 | 显示全部楼层
滴滴滴滴的 发表于 2014-6-15 08:09
为啥我所在的用户组,我法查看或下载附件

〖新手必读〗之如何获取积分,提高权限(新)http://forum.vibunion.com/thread-132101-1-1.html
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

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

GMT+8, 2024-3-28 18:04 , Processed in 0.076015 second(s), 18 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

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