王济 matlab,关于王济老师《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

全部为按照书中代码编程，为什么运行不了，望大神帮忙！

附件为原时程信号，方便大家试运行。