chengxu.rar_模态分析

%利用拟合圆法，求三阶固有频率、阻尼比。 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear clc close all hidden %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% global mn; %声明为全局变量 z=2; %采用z号激振点 fs=4096; %采样频率 N=4096; %采样点数 mn=1; %模态阶数 df=fs/N; %频率间隔 f0=[72 150 258]; %模态频率初值数组 fid=fopen(['pinxianghanshu',num2str(z),'.txt'],'r'); H=fscanf(fid,'%f',[2,inf]); %读入频响函数实部与虚部数据 fclose(fid); %定义离散频率向量 f=0:df:(length(H(1,:))-1)*df; %f=0:（频响函数数据点个数-1)(4095?) %定义离散圆频率向量 w=2*pi*f; %建立模态参数向量 fds=zeros(1,3*mn); %zeros(m,n)产生m×n的double类零矩阵 for j=1:mn %确定拟合导纳圆的数据 k=round(f0(j)/df)+1; % round,四舍五入函数 [m,kc]=max(abs(H(2,k-4:k+4))); % abs= Absolute value % kc为最大值所在的位置 [Y,I] = max([1 2 3 3 22 5]) Y = 22 I = 5 kc=kc+k-5; %设频响函数虚部模态峰值的0.15作为取数据的条件 v=0.05*abs(H(2,kc)); %依照条件进行取数据 for k=1:10 l=kc-k; if abs(H(2,l))>v&&abs(H(2,l))<=abs(H(2,l+1)) k1=l; else break; end end k2=kc+1; %这条语句是后来加的 for k=1:10 l=kc+k; if abs(H(2,l))>v&&abs(H(2,l))<=abs(H(2,l-1)) k2=l; else break; end end %将确定的频响函数的实部和虚部分别存入x和y x=H(1,k1:k2); y=H(2,k1:k2); %用最小二乘法进行导纳圆的拟合 A=[sum(x.^2),sum(x.*y),sum(x);sum(x.*y),sum(y.^2),sum(y);sum(x),sum(y),length(x)]; B=[-sum(x.^3+x.*y.^2);-sum(x.^2.*y+y.^3);-sum(x.^2+y.^2)]; %解线性方程 c=inv(A)*B; %求逆矩阵 %计算圆心X坐标 x0=-c(1)/2; %计算圆心Y坐标 y0=-c(2)/2; %计算圆半径 r=sqrt(x0^2+y0^2-c(3)); %sqrt开根号 alpha=0:pi/100:2*pi; %角度[0,2*pi] xx=r*cos(alpha)+x0; yy=r*sin(alpha)+y0; switch j case 1 figure(1) case 2 figure(2) case 3 figure(3) end plot(xx,yy); axis equal hold on plot(x0,y0,'r*'); plot(x,y,'r*'); switch j case 1 title('一阶拟合圆') case 2 title('二阶拟合圆') case 3 title('三阶拟合圆') end grid on %用抛物线插值法计算模态频率值 %y1、y2、y3分别为频响函数虚部模态峰值及左右的3个点 %x1、x2、x3是对应于y1、y2、y3的频率值 x1=f(kc-1); x2=f(kc); x3=f(kc+1); y1=H(2,kc-1); y2=H(2,kc); y3=H(2,kc+1); fc=((y2-y1)*(x3^2-x2^2)-(y3-y2)*(x2^2-x1^2))/(2*(2*y2-y1-y3)*(x2-x1)); yc=(y1*(x2-fc)^2-y2*(x1-fc)^2)/((x2-fc)^2-(x1-fc)^2); if abs(y1)>abs(y3) xc=H(1,kc-1)+(H(1,kc)-H(1,kc-1))*(fc-x1)/(x2-x1); else xc=H(1,kc)+(H(1,kc+1)-H(1,kc))*(fc-x2)/(x3-x2); end F(j)=fc; %固有频率 %计算对应模态频率的阻尼比值 [m,k]=max(abs(y)); %根据三角形的余弦定理计算峰值数据点与前后数据点的夹角 a=(xc-x(k-1))^2+(yc-y(k-1)^2); b=(x0-xc)^2+(y0-yc)^2; c=(x0-x(k-1))^2+(y0-y(k-1))^2; a1=acos((b+c-a)/(2*sqrt(b)*sqrt(c)))/2; a=(x(k+1)-xc)^2+(y(k+1)-yc)^2; b=(x0-x(k+1))^2+(y0-y(k+1))^2; c=(x0-xc)^2+(y0-yc)^2; a2=acos((b+c-a)/(2*sqrt(b)*sqrt(c)))/2; D(j)=(f(k+1)-f(k-1))/(F(j)*(tan(a1)+tan(a2)));%阻尼比 end %打开文件输出模态参数数据 fid=fopen(['motaicanshu',num2str(z),'.txt'],'w'); fprintf(fid,' 阶次 频率（Hz） 阻尼比（%%）\n'); for k=1:mn fprintf(fid,'%d %10.4f %10.4f\n',k,F(k),D(k)*100.0); end fclose(fid);