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',[1,mn]);
fno=fscanf(fid,'%s',1);
H0=fscanf(fid,'%f',[2,inf]);
status=fclose(fid);
f=0:df:(length (H0(1,:))-1)*df;
w=2*pi*f;
fds=zeros(1,3*mn);
for n=1:4
for j=1:mn
l=3*(j-1);
fds(l+1:l+3)=zeros(1,3);
H1=fun81(fds,w);
H=H0-H1;
k=round(f0(j)/df)+1;
[m,kc]=max(abs(H(2,k-4:k+4)));
kc=kc+k-5;
v=0.15*abs(H(2,kc));
for k=2:30
l=kc-k;
if abs(H(2,l))>v&abs(H(2,l))<=abs(H(2,l+1))
k1=l;
else
break;
end
end
for k=2:30
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=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;
x0=-c(1)/2;
y0=-c(2)/2;
r=sqrt(x0^2+y0^2-c(3));
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(l+1)-f(l-1))/(F(j)*(tan(a1)+tan(a2)));
S(j)=yc*2*D(j);
l=3*(j-1);
fds(l+1:l+3)=[2*pi*F(j) D(j) S(j)];
end
end
H1=fun81(fds,w);
%绘制频响函数实部拟合曲线
figure(1)
subplot(2,1,1);
plot(f,H0(1,:),':',f,H1(1,:));
xlabel('频率（HZ）');
ylabel('实部');
legend('实测','拟合');
grid on;
%绘制频响函数虚部拟合曲线
subplot(2,1,2);
plot(f,H0(2,:),':',f,H1(2,:));
xlabel('频率（HZ）');
ylabel('虚部');
legend('实测','拟合');
grid on;
%绘制频响函数导纳拟合曲线
figure(2)
plot(H0(1,:),H0(2,:),':',H1(1,:),H1(2,:));
axis('equal');
xlabel('实部');
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,S(k));
end
status=fclose(fid);
