|
楼主 |
发表于 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] = ...
|
|