有理分式多项式法（Levy）拟合问题

simplebinbin

用有理分式多项式法进行模态参数分析发现拟合曲线漏掉了一些模态这是什么原因造成的呢，无论怎么改变分析模态的阶数，都是这样。
请高手指点

simplebinbin

从图中可以看出，1、3、4阶模态都没有拟合上

simplebinbin

出现了下面的警告：
Warning: Matrix is close to singular or badly scaled.
         Results may be inaccurate. RCOND = 3.526498e-024.
矩阵要怎么修改呢

simplebinbin

1. clear;format long
   
2. %%%%%%%%%%%%%%%%%%%%%%%%%
   
3. %声明为全局变量
   
4. global mn;
   
5. fni=input('有理分式多项式法-输入数据批文件名：','s');  
   
6. fid=fopen(fni,'r');
   
7. mn=fscanf(fid,'%d',1);        %模态阶数
   
8. df=fscanf(fid,'%f',1);        %频率间隔
   
9. fno=fscanf(fid,'%s',1);       %输出数据文件名   
   
10. b=fscanf(fid,'%f',[2,inf]);   %输出实测频响函数实虚部数值数组
    
11. status=fclose(fid);
    
12. %定义幂多项式的阶次
    
13. nm=mn*2;
    
14. %取频响函数长度
    
15. n=length(b(1,:));
    
16. %建立离散频率向量
    
17. f=0:df:(n-1)*df;
    
18. %建立离散圆频率向量
    
19. w=2*pi*f;
    
20. %建立归一化离散频率向量
    
21. wi=w/max(w);
    
22. %建立实测频响函数复数向量
    
23. H=b(1,:)+i*b(2,:);
    
24. %计算拟合频响函数的分子和分母系数向量
    
25. [A,B]=invfreqs(H,wi,nm,nm,[],100);
    
26. %幂多项式方程求根
    
27. P=roots(B);
    
28. %计算模态频率向量
    
29. F1=abs(P)*max(w)/(2*pi);
    
30. %计算阻尼比向量
    
31. D1=-real(P)./(abs(P));
    
32. %计算振型系数向量
    
33. for k=1:mn
    
34.     if k==1
    
35.         p(1:nm-1)=P(2:nm);
    
36.     else
    
37.         p(1:k-1)=P(1:k-1);
    
38.         p(k:nm-1)=P(k+1:nm);
    
39.     end
    
40.     y=poly(p);
    
41.     S1(k)=polyval(A,P(k))/polyval(y,P(k));
    
42. end
    
43. %计算拟合的频响函数
    
44. H1=freqs(A,B,wi);
    
45. %绘制频响函数实部拟合曲线图
    
46. figure(2)
    
47. nn=1:n;
    
48. subplot(2,1,1);
    
49. plot(f,real(H(nn)),':',f,real(H1(nn)),'r');
    
50. xlabel('频率（Hz）');
    
51. ylabel('实部');
    
52. legend('实测','拟合');
    
53. grid on;
    
54. %绘制频响函数虚部拟合曲线图
    
55. subplot(2,1,2);
    
56. plot(f,b(2,:),':',f,imag(H1(nn)),'r');
    
57. xlabel('频率(Hz)');
    
58. ylabel('虚部');
    
59. legend('实测','拟合');
    
60. grid on;
    
61. %将模态频率从小到大排序
    
62. [F2,I]=sort(F1);
    
63. %剔除方程解中的非模态项和共轭项
    
64. m=0;
    
65. for k=1:nm-1
    
66.     if F2(k)~=F2(k+1)
    
67.         continue;
    
68.     end
    
69.     m=m+1;
    
70.     l=I(k);
    
71.     F(m)=F1(l);  %模态频率
    
72.     D(m)=D1(l);  %阻尼比
    
73.     S(m)=S1(l);  %振型系数
    
74. end
    
75. %打开文件输出模态参数数据
    
76. fid=fopen(fno,'w');
    
77. fprintf(fid,'频率（Hz） 阻尼比（%%） 振型系数\n');
    
78. for k=1:mn
    
79.    fprintf(fid,'%10.4f %10.4f %10.6f\n',F(k),D(k)*100.0,S(k));
    
80. end
    
81. status=fclose(fid);
    
82. 
    
83. 这是程序代码

复制代码

simplebinbin