王济《Matlab在振动信号处理中的应用》书中的疑问

jiangong

在这本王济——《Matlab在振动信号处理中的应用》书中第八章8.3节最小二乘迭代法 程序“8.2a最小二乘法模态参数识别"程序运行时老出错,不知道为什么？因为这种方法对自己很重要，想强高手指点一下，谢谢！

yuebeifan

至少把你的错误发上来啊

jiangong

回复 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',[2,inf]); %实测频响函数实部虚部数据
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
         [R,permR] = feval(pcmtx,H,pcoptions,DM,DG,varargin{:});

Error in ==> snls at 334
      [sx,snod,qp,posdef,pcgit,Z] = trdog(x,g,A,D,delta,dv,...

Error in ==> lsqncommon at 153
    [xC,FVAL,LAMBDA,JACOB,EXITFLAG,OUTPUT,msg]=...

Error in ==> lsqcurvefit at 258
[xCurrent,Resnorm,FVAL,EXITFLAG,OUTPUT,LAMBDA,JACOB] = ...

yuebeifan

从字面上看到qr不存在，但是你的程序没有找到qr，难道是其他的子程序？

jiangong

回复 3 # jiangong 的帖子

应该是吧，这一块我也看不明白

tenglang

help qr

你会看到下面的内容,正交三角分解的东东,
楼主是不是用的绿色版的matlab, 如果是,先help一下,然后就能运行了. 如果还不行,就要看输入该函数参数正不正确.

QR     Orthogonal-triangular decomposition.
    [Q,R] = 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.

    [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 [Q,R] = QR(A).

    If A is full:

    [Q,R,E] = 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.

    [Q,R,E] = 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 [Q,R] = 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).

    [Q,R,E] = 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.

    [Q,R,E] = QR(A,0) produces an "economy size" decomposition in which E
    is a permutation vector, so that A(:,E) = Q*R.

    [C,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.

    [C,R,E] = QR(A,B), also returns a fill-reducing ordering.
    The least-squares solution to A*X = B is X = E*(R\C).

    [C,R] = 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 [C,R] = QR(A,B).

    [C,R,E] = 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;

navy423

我个感觉这本书上很多程序有点问题，即使是书上的源程序都无法运行。不知何解？

wuzhijun117420

这也正常。不少教程都是这样的。这就要我们用心先研究。

jiangong

回复 8 # wuzhijun117420 的帖子

但是很奇怪的就是如果用06版本的matalb就可以运行，我刚开始用09的你就不行

wuzhijun117420

个人感觉是lsqcurvefit调用格式存在问题。help lascurvefit
X=LSQCURVEFIT(FUN,X0,XDATA,YDATA) starts at X0 and finds coefficients X
    to  best fit the nonlinear functions in FUN to the data YDATA (in the
    least-squares sense).
希望楼主参考

wuzhijun117420

不知道将w令为整体变量，而不作为传递的参数是否可以。

jiangong

回复 11 # wuzhijun117420 的帖子

好，我再试试吧啊，谢谢你们:@)

longdesigner

这本书怎么样啊，我做的是地震方面的，求相干函数的相关知识 有帮助吗？

梦泉

正在学习这边书,不知道怎么获取这本书的源程序？
那位大侠能把源程序发给我一下：1634208021@qq.com
小弟不胜感激！

leooja

你好，关于发帖所涉王济老师书中模态识别程序8.2a出现提示“Complex sparse QR is not yet available”该如何解决