%利用擬合圓法，求三階固有頻率、阻尼比。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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