24小时热门版块排行榜

返回列表

leooja

铁虫 (初入文坛)

应助: 0 (幼儿园)
金币: 13.6
帖子: 12
在线: 4.9小时
虫号: 3365819
注册: 2014-08-14
专业: 结构工程

[求助] 关于王济老师《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
全部为按照书中代码编程，为什么运行不了，望大神帮忙！
附件为原时程信号，方便大家试运行。

回复此楼