王济《Matlab在振动信号处理中的应用》书中的疑问
在这本王济——《Matlab在振动信号处理中的应用》书中第八章8.3节最小二乘迭代法 程序“8.2a最小二乘法模态参数识别"程序运行时老出错,不知道为什么？因为这种方法对自己很重要，想强高手指点一下，谢谢！ 至少把你的错误发上来啊 回复 2 # yuebeifan 的帖子%程序8-2a
%最小二乘法模态参数识别（复模态-频响函数实虚部）
%%%%%%%%%%%%%%%%%%%%%%%
clear
clc
close all hidden
format long
%%%%%%%%%%%%%%%%%%%%%%%
%声明全局变量
global mn
fni=input('最小二乘法模态参数识别-输入数据文件名:','s');
fid=fopen(fni,'r');
mn = fscanf(fid,'%d',1); %模态阶数
df = fscanf(fid,'%f',1); %频率间隔
f0 = fscanf(fid,'%f',mn); %输入模态频率初值数组
d0 = fscanf(fid,'%f',mn); %输入模态阻尼比初值数组
fno = fscanf(fid,'%s',1); %输出数据文件名
b = fscanf(fid,'%f',); %实测频响函数实部虚部数据
status=fclose(fid);
%建立离散频率向量
f=0:df:(length(b(1,:))-1)*df;
%建立离散圆频率向量
w=2*pi*f;
%建立实测频响函数复数向量
H =b(1,:)+b(2,:)*i;
%计算模态圆频率初值向量
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
%用最小二乘非线性数据拟合法估计复模态参数
x=lsqcurvefit('fun82a',x0,w,H);
%计算模态频率、阻尼比和振型系数
for j=1:mn
l=4*(j-1);
c=x(l+1)+i*x(l+2);
d=x(l+3)+i*x(l+4);
F(j)=abs(c)/(2*pi);%模态频率
D(j)=-real(c)/abs(c);%阻尼比
S(j)=d; %复振型系数
end
%计算拟合的频响函数向量
H1=fun82a(x,w);
%绘制频响函数实部拟合曲线图
subplot(2,1,1);
plot(f,real(H),':',f,real(H1));
xlabel('频率 (Hz)');
ylabel('实部');
legend('实测','拟合');
grid on;
%绘制频响函数虚部拟合曲线图
subplot(2,1,2);
plot(f,imag(H),':',f,imag(H1),'r');
xlabel('频率 (Hz)');
ylabel('虚部');
legend('实测','拟合');
grid on;
%打开文件输出模态参数数据
fid=fopen(fno,'w');
fprintf(fid,' 频率(Hz) 阻尼比(%%) 振型系数\n');
for k=1:mn
fprintf(fid,'%10.4f%10.4f%10.6f\n',F(k),D(k)*100.0,imag(S(k)));
end
status=fclose(fid);
%函数8-2a(用于程序8-2a)
function m=fun82a(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+i*x4)./(w*i-(x1+i*x2))+(x3-i*x4)./(w*i-(x1-i*x2)));
end
就是这个程序，导入数据文件后运行的结果为
??? Error using ==> qr
Complex sparse QR is not yet available.
Error in ==> aprecon at 57
RPCMTX = qr(TM(:,p));
Error in ==> trdog at 47
= feval(pcmtx,H,pcoptions,DM,DG,varargin{:});
Error in ==> snls at 334
= trdog(x,g,A,D,delta,dv,...
Error in ==> lsqncommon at 153
=...
Error in ==> lsqcurvefit at 258
= ...
从字面上看到qr不存在，但是你的程序没有找到qr，难道是其他的子程序？ 回复 3 # jiangong 的帖子
应该是吧，这一块我也看不明白
本帖最后由 tenglang 于 2010-12-22 16:39 编辑
help qr
你会看到下面的内容,正交三角分解的东东,
楼主是不是用的绿色版的matlab, 如果是,先help一下,然后就能运行了. 如果还不行,就要看输入该函数参数正不正确.
QR Orthogonal-triangular decomposition.
= QR(A), where A is m-by-n, produces an m-by-n upper triangular
matrix R and an m-by-m unitary matrix Q so that A = Q*R.
= QR(A,0) produces the "economy size" decomposition.
If m>n, only the first n columns of Q and the first n rows of R are
computed. If m<=n, this is the same as = QR(A).
If A is full:
= QR(A) produces unitary Q, upper triangular R and a
permutation matrix E so that A*E = Q*R. The column permutation E is
chosen so that ABS(DIAG(R)) is decreasing.
= QR(A,0) produces an "economy size" decomposition in which E
is a permutation vector, so that A(:,E) = Q*R.
X = QR(A) and X = QR(A,0) return the output of LAPACK's *GEQRF routine.
TRIU(X) is the upper triangular factor R.
If A is sparse:
R = QR(A) computes a "Q-less QR decomposition" and returns the upper
triangular factor R. Note that R = CHOL(A'*A). Since Q is often nearly
full, this is preferred to = QR(A).
R = QR(A,0) produces "economy size" R. If m>n, R has only n rows. If
m<=n, this is the same as R = QR(A).
= QR(A) produces unitary Q, upper triangular R and a
permutation matrix E so that A*E = Q*R. The column permutation E is
chosen to reduce fill-in in R.
= QR(A,0) produces an "economy size" decomposition in which E
is a permutation vector, so that A(:,E) = Q*R.
= QR(A,B), where B has as many rows as A, returns C = Q'*B.
The least-squares solution to A*X = B is X = R\C.
= QR(A,B), also returns a fill-reducing ordering.
The least-squares solution to A*X = B is X = E*(R\C).
= QR(A,B,0) produces "economy size" results. If m>n, C and R have
only n rows. If m<=n, this is the same as = QR(A,B).
= QR(A,B,0) additionally produces a fill-reducing permutation
vector E.In this case, the least-squares solution to A*X = B is
X(E,:) = R\C.
Example: The least squares approximate solution to A*x = b can be found
with the Q-less QR decomposition and one step of iterative refinement:
if issparse(A), R = qr(A); else R = triu(qr(A)); end
x = R\(R'\(A'*b));
r = b - A*x;
e = R\(R'\(A'*r));
x = x + e;
我个感觉这本书上很多程序有点问题，即使是书上的源程序都无法运行。不知何解？ 这也正常。不少教程都是这样的。这就要我们用心先研究。 回复 8 # wuzhijun117420 的帖子
但是很奇怪的就是如果用06版本的matalb就可以运行，我刚开始用09的你就不行 个人感觉是lsqcurvefit调用格式存在问题。help lascurvefit
X=LSQCURVEFIT(FUN,X0,XDATA,YDATA) starts at X0 and finds coefficients X
tobest fit the nonlinear functions in FUN to the data YDATA (in the
least-squares sense).
希望楼主参考 不知道将w令为整体变量，而不作为传递的参数是否可以。 回复 11 # wuzhijun117420 的帖子
好，我再试试吧啊，谢谢你们:@) 这本书怎么样啊，我做的是地震方面的，求相干函数的相关知识 有帮助吗？ 正在学习这边书,不知道怎么获取这本书的源程序？
那位大侠能把源程序发给我一下：1634208021@qq.com
小弟不胜感激！ 你好，关于发帖所涉王济老师书中模态识别程序8.2a出现提示“Complex sparse QR is not yet available”该如何解决
页:
[1]