马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。
您需要 登录 才可以下载或查看，没有账号？我要加入
x
本帖最后由 牛小贱 于 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',[ig,inf]); %实测时域或频域数据
- 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');
- %计算特征值及特征向量
- [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=0:(2*n-1)
- Va(k+1,:)=[conj(U).^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.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)=[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:49
提升卡
变色卡
显身卡
楼主
发表于 2014-6-14 13:06
已经不行了。。