声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 3990|回复: 17

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

  [复制链接]
发表于 2010-12-21 11:00 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

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

使用道具 举报

发表于 2010-12-21 16:51 | 显示全部楼层
至少把你的错误发上来啊
 楼主| 发表于 2010-12-21 19:26 | 显示全部楼层
回复 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] = ...
发表于 2010-12-21 21:29 | 显示全部楼层
从字面上看到qr不存在,但是你的程序没有找到qr,难道是其他的子程序?
 楼主| 发表于 2010-12-22 10:25 | 显示全部楼层
回复 3 # jiangong 的帖子

应该是吧,这一块我也看不明白
发表于 2010-12-22 16:21 | 显示全部楼层
本帖最后由 tenglang 于 2010-12-22 16:39 编辑

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;
发表于 2010-12-22 21:05 | 显示全部楼层
我个感觉这本书上很多程序有点问题,即使是书上的源程序都无法运行。不知何解?
发表于 2010-12-23 10:35 | 显示全部楼层
这也正常。不少教程都是这样的。这就要我们用心先研究。
 楼主| 发表于 2010-12-23 14:56 | 显示全部楼层
回复 8 # wuzhijun117420 的帖子

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

点评

有时新旧版用法有差, 书上原则上应会说明使用版别!  发表于 2010-12-28 00:10
发表于 2010-12-25 09:40 | 显示全部楼层
个人感觉是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).
希望楼主参考
发表于 2010-12-25 09:42 | 显示全部楼层
不知道将w令为整体变量,而不作为传递的参数是否可以。
 楼主| 发表于 2010-12-25 10:17 | 显示全部楼层
回复 11 # wuzhijun117420 的帖子

好,我再试试吧啊,谢谢你们:@)
发表于 2011-3-30 11:13 | 显示全部楼层
这本书怎么样啊,我做的是地震方面的,求相干函数的相关知识 有帮助吗?
发表于 2011-7-23 13:03 | 显示全部楼层
正在学习这边书,不知道怎么获取这本书的源程序?
那位大侠能把源程序发给我一下:1634208021@qq.com
小弟不胜感激!
发表于 2014-9-4 16:20 | 显示全部楼层
你好,关于发帖所涉王济老师书中模态识别程序8.2a出现提示“Complex sparse QR is not yet available”该如何解决
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-8 23:00 , Processed in 0.078121 second(s), 23 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表