| 查看: 2049 | 回复: 4 | ||
[求助]
关于王济老师《Matlab在振动信号处理中的应用》书中第八章8.3节 程序8.2a的问题
|
|
如题用该方法识别结构自振模态，总会提示：“??? Error using ==> qr Complex sparse QR is not yet available.” 程序原代码如下： %随机信号谱分析 %%%%%%%%%%%%%%%%%%%% clear clc close all hidden format global mn %%%%%%%%%%%%%%%%%%% %打开数据文件 [filename,pathname]=uigetfile('*.xls','selectfile','MultiSelect', 'on'); l=length(filename); fun=input('function type'); %(谱分析类型，1=自谱，2=互谱，3=频响，4=相干） windos=input('windos type'); %（窗函数类型，1=矩形，2=汉宁，3=海宁，4=布莱克曼） sf=input('simpling frequency'); %（采样频率） N=input('lenght of fft'); %(FFT长度） f0=[4,4.17,5,12.5,14.29,16.68]; %（模态频率初始值，估计情况确定） d0=[0.05,0.05,0.05,0.05,0.05,0.05]; %(模态阻尼比初始值) mn=6; %建立输出矩阵 outmodex=zeros(l,18); outmodey=zeros(l,18); for ii=1:l; openpath=strcat(pathname,filename{ii}); [data,text]=xlsread(openpath); [m,n]=size(data); tt=0:1/sf  m-1)/sf;tt=tt'; %消除趋势项 for jj=4:n p3=polyfit(tt,data(:,jj),3); data(:,jj)=data(:,jj)-polyval(p3,tt); end x1=data(:,2); %(x为输入激励） x2=data(:,3); y1=data(:,4)-data(:,2); %（y为相对时程） y2=data(:,5)-data(:,3); %频率向量 f=0:sf/N sf/2-sf/N);switch windos case 1 %（矩形窗） w=boxcar(N); case 2 w=hanning(N); %（汉宁窗） case 3 w=hamming(N); %（海宁窗） case 4 w=triang(N); %（三角窗） otherwise w=boxcar(N); end switch fun case 1 %（自谱） z1=spectrum.psd(y1,N,sf,w,N/2); z2=spectrum.psd(y2,N,sf,w,N/2); case 2 %（互谱） z1=csd(x1,y1,N,sf,w,N/2); z2=csd(x2,y2,N,sf,w,N/2); case 3 %（频响） z1=tfestimate(x1,y1,w,N/2,N,256); z2=tfestimate(x2,y2,w,N/2,N,256); case 4 z1=cohere(x1,y1,N,sf,w,N/2); z2=cohere(x2,y2,N,sf,w,N/2); otherwise end figure(ii) %（可附加逻辑判断决定绘制及输出数据） nn=1:N/4; subplot(2,2,1); plot(f(nn),real(z1(nn))); xlabel('Hz'); ylabel('real'); subplot(2,2,2); plot(f(nn),imag(z1(nn))); xlabel('Hz'); ylabel('imag'); subplot(2,2,3); plot(f(nn),real(z2(nn))); xlabel('Hz'); ylabel('real'); subplot(2,2,4); plot(f(nn),imag(z2(nn))); xlabel('Hz'); ylabel('imag'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%55 %最小二乘法计算模态，阻尼比及模态系数 ff=0:0.5*sf/(length(z1)):0.5*(length(z1)-1)*sf/(length(z1)); ww=2*pi*ff; W0=2*pi*f0; for j=1:mn; L=4*(j-1); X0(L+1:L+4)=[-W0(j)*d0(j),W0(j)*sqrt(1-d0(j)^2),1,1]; end z1=z1'; z2=z2'; X1=lsqcurvefit('fun82',X0,ww,z1); X2=lsqcurvefit('fun82',X0,ww,z2); for j=1:mn; L=4*(j-1); %X向模态分析 C1=X1(L+1)+1i*X1(L+2); D1=X1(L+3)+1i*X1(L+4); outmodex(ii,j)=abs(C1)/(2*pi);%模态频率 outmodex(ii,j+6)=-real(C1)/abs(C1);%模态阻尼 outmodex(ii,j+12)=D1; %负振型系数 %Y向模态分析 C2=X2(L+1)+1i*X2(L+2); D2=X2(L+3)+1i*X2(L+4); outmodey(ii,j)=abs(C2)/(2*pi); outmodey(ii,j+6)=-real(C2)/abs(C2); outmodey(ii,j+12)=D2; end outmode=[outmodex;zeros(1,18);outmodey]; modepath=strcat(pathname,'model',filename{ii}); [status,message]=xlswrite(modepath,outmode); end 其中“fun82”代码如下： function M=fun82(X,W) %通过模态参数计算拟合频响函数 %输入参数 %X-复模态参数向量 %W-频率变量向量 %输出参数 %M-拟合频响函数 global mn M=zeros(1,length(W)); for K=0:mn-1; L=4*K; X1=X(L+1);%特征值实部 X2=X(L+2);%特征值虚部 X3=X(L+3);%特征值向量实部 X4=X(L+4);%特征值向量虚部 M=M-W.^2.*((X3+1i*X4)./(W*1i-(X1+1i*X2))+(X3-1i*X4)./(W*1i-(X1-1i*X2))); end 全部为按照书中代码编程，为什么运行不了，望大神帮忙！ 附件为原时程信号，方便大家试运行。 |
回复此楼» 本帖附件资源列表
-
欢迎监督和反馈：小木虫仅提供交流平台，不对该内容负责。
本内容由用户自主发布，如果其内容涉及到知识产权问题，其责任在于用户本人，如对版权有异议，请联系邮箱：xiaomuchong@tal.com - 附件 1 : w1.xlsx
- 附件 2 : w2.xlsx
2014-09-04 15:59:58, 4.08 M
2014-09-04 16:00:37, 4.08 M
» 猜你喜欢
0855求调剂材料
已经有12人回复
-
297分083200求助
已经有7人回复
-
070300化学学硕311分求调剂
已经有10人回复
-
专硕0854初试考材科基，求调剂
已经有8人回复
-
0857大类环境工程B区求调剂
已经有3人回复
-
0817化学工程与技术求调剂，一志愿中海洋319
已经有14人回复
-
0703化学调剂325分
已经有11人回复
-
286求调剂
已经有8人回复
-
求调剂
已经有5人回复
-
生物学308分求调剂（一志愿华东师大）
已经有4人回复
hytao2012
铁杆木虫 (正式写手)
木头虫子
- 应助: 53 (初中生)
- 金币: 6343.8
- 散金: 115
- 红花: 11
- 帖子: 479
- 在线: 206.6小时
- 虫号: 2050091
- 注册: 2012-10-08
- 性别: GG
- 专业: 固体力学
2楼2014-09-04 16:12:56
surrounding8
铜虫 (初入文坛)
- 应助: 1 (幼儿园)
- 金币: 136.6
- 红花: 1
- 帖子: 20
- 在线: 15.6小时
- 虫号: 3434106
- 注册: 2014-09-22
- 专业: 信号理论与信号处理
3楼2014-09-22 20:55:45
pdikang
新虫 (小有名气)
- 应助: 0 (幼儿园)
- 金币: 651.1
- 散金: 317
- 帖子: 203
- 在线: 30.4小时
- 虫号: 3082556
- 注册: 2014-03-23
- 专业: 传动机械学
4楼2014-11-23 21:13:49
wanlingstar
新虫 (正式写手)
- 应助: 0 (幼儿园)
- 金币: 3415.1
- 帖子: 318
- 在线: 49.4小时
- 虫号: 2686341
- 注册: 2013-09-27
- 专业: 岩土与基础工程
5楼2017-01-15 23:51:13
