简支梁模态分析程序-motaifenxi_2.m

%利用拟合圆法，求三阶固有频率、阻尼比。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
clc
close all hidden
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global mn;                           %声明为全局变量
z=3;                                 %采用z号激振点
fs=4096;                             %采样频率
N=4096;                              %采样点数
mn=1;                                %模态阶数
df=fs/N;                             %频率间隔
f0=[26 150 258];                     %模态频率初值数组
fid=fopen（[‘pinxianghanshu‘num2str（z）‘.txt‘]‘r‘）;
H=fscanf（fid‘%f‘[2inf]）;         %读入频响函数实部与虚部数据
fclose（fid）;
%定义离散频率向量
f=0:df:（length（H（1:））-1）*df;
%定义离散圆频率向量
w=2*pi*f;
%建立模态参数向量
fds=zeros（13*mn）;
    for j=1:mn
        %确定拟合导纳圆的数据
        k=round（f0（j）/df）+1;
        [mkc]=max（abs（H（2k-4:k+4）））;
        kc=kc+k-5;
        %设频响函数虚部模态峰值的0.15作为取数据的条件
        v=0.05*abs（H（2kc））;
        %依照条件进行取数据
        for k=1:10
            l=kc-k;
            if abs（H（2l））>v&abs（H（2l））<=abs（H（2l+1））
                k1=l;
            else
                break;
            end
        end
        for k=1:10
            l=kc+k;
            if abs（H（2l））>v&abs（H（2l））<=abs（H（2l-1））
                k2=l;
            else
                break;
            end
        end
        %将确定的频响函数的实部和虚部分别存入x和y
        x=H（1k1:k2）;
        y=H（2k1: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））;
        alpha=0:pi/100:2*pi;            %角度[02*pi] 
        xx=r*cos（alpha）+x0