%%%%%%%%%%%%%%%%%%%%
%
clear
clc
close all hidden
format short
%%%%%%%%%%%%%%%%%%%%%%
fni=input('ITD法模態參數識別-輸入數據文件名','s');
fid=fopen(fni,'r');
mn=fscanf(fid,'%d',1);%
%
%
%
ig=fscanf(fid,'%d',1);
%
f=fscanf(fid,'%f',1);
fno=fscanf(fid,'%s',1);%
b=fscanf(fid,'%f',[ig,inf]);%
status=fclose(fid);
%
nm=2*mn;
%
if ig==1%
    %
    sf=f;
    %
    n=fix(length(b)/2);
    %
    h=b(1,1:2*n)';
    %
    dt=1/sf;
    %
    t=0:dt:(2*n-1)*dt;
else
    %
    %
    df=f;
    n=length(b(1,:));
    %
    f=0:df:(n-1)*df;
    %
    H=b(1,:)'+b(2,:)'*1i;
    H(n+1)=real(H(n));
    H(n+2:2*n)=conj(H(n:-1:2));
    %
    h=real(ifft(H));
    %
    t=linspace(0,1/df,2*n);
    %
    dt=t(2)-t(1);
end
L=length(h);
M=L/2;
for k=1:nm
    x1(k,:)=h(k:L-(nm-k+1))';
    x2(k,:)=h(k+1:L-(nm-k))';
end
%
B=x1\x2;%
%
[A,V]=eig(B);
for k=1:nm
    U(k)=V(k,k);
end
%f
F1=abs(log(U'))./(2*pi*dt);
%a
D1=sqrt(1./(((imag(log(U'))./real(log(U'))).^2)+1));
%
l=1;
for k=0:(2*n-1)
    Va(k+1,:)=[conj(U).^k];
end
S1=(inv(conj(Va')*Va))*conj(Va')*h;
%
h1=real(Va*S1);
figure(1)
plot(t,h,'r',t,h1)
xlabel('時間（s）');
ylabel('幅值');
legend('實測','擬合');
hold on
grid on
if ig>1
    H1=fft(Va*S1);
    figure(2)
    nn=1:n;
    subplot(2,1,1);
    plot(f,real(H(nn)),':',f,real(H1(nn)),'r');
    hold on;
    xlabel('頻率Hz');
    ylabel('實部');
    legend('實測','擬合');
    grid on
    %
    subplot(2,1,2);
    plot(f,b(2,:),':',f,imag(H1(nn)),'r');
    xlabel('頻率（Hz）');
    ylabel('虛部');
    legend('實測','擬合');
    grid on
end

%
[F2,I]=sort(F1);
%
m=0;
for k=1:nm-1
    if F2(k)~=F2(k+1)
        continue;
    end
    m=m+1;
    l=I(k);
    F(m)=F1(l)%
    D(m)=D1(l)%
    S(m)=S1(l);%
end
%
fid=fopen(fno,'w');
fprintf(fid,'頻率（Hz）  阻尼比（%%）  陣型系數\n');
for k=1:m;
    fprintf(fid,'%10.4f %10.4f  %10.6f\n ',F(k),D(k)*100,imag(S(k)));
end

status=fclose(fid);