ex09051.m

复指数法是根据结构的自由振动响应或脉冲响应函数可以表示为复指数函数和的形式

%
% ex09051
clear all; clc;
mn=7;
ig=2;
f=0.1953125;
b=load('x09051.txt');       %读入数据
b=b';
fno='out09051.txt';      
%b=fscanf (fid,'%f', Eig, inf]); 
%status=fclose (fid);
nm=2 * mn;
rad=pi/180;
if ig==1
   sf=f;
%1~ 1/2 I~:f~
   n=fix(length(b)/2);
   h=b(1, 1 : 2*n)';
   dt=1/sf;
   t=0 : dt : (2*n-1)*dt;
else      %~:~t~t~
   df=f;
   n=length(b(1,:));
   f=0 : df : (n-1)*df;
   H=b(1,:)'+b(2,:)'*i;
   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: M-1+k);
 end
for k=1: M
  x2(k,:) =-h(nm+k);
end
B=x1\x2;   %B=inv (x1'*x1) * (x1'*x2);
%~ B (nm+1) ~ 1
B(nm+1) =1;
B1=B(nm+1 : -1: 1);
V=roots(B1);
F1=abs(log(V)) / (2*pi*dt);
D1=sqrt(1./(((imag(log(V))./real(log(V))).^2)+1));
for k=0: (2*n-1)
   Va(k+1,:) = conj(V').^k;
end
S1=(inv(conj(Va') * Va) *conj(Va') * h);

h1=real(Va * S1);

figure (1)
plot(t, h,'r','linewidth',2); hold on;
plot(t,h1,'k'); hold off;
xlabel ('Time (s)');
ylabel ('Amplitude');
legend ('Measured','Fitness');
grid on;

if ig>1
   H1=fft(Va * S1);
   figure(2)
   nn=1:  n;
   subplot (2, 1, 1);
   plot(f, real(H(nn)),'r', 'linewidth',2); hold on; 
   plot(f, real(H1(nn)),'k'); hold off;
   xlabel ('Frequency (Hz)');
   ylabel ('Real Part');
   legend ('Measured','Fitness');
   grid on;

  subplot (2, 1, 2);
  plot(f, b(2,:),'r', 'linewidth',2); hold on;
  plot(f, imag(H1(nn)),'k'); hold off;
  xlabel ('Frequency (Hz)');
  ylabel ('Imag Part');
  legend ('Measured','Fitness');
  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);  %[~g
  D(m) =D1(l);   %l~l]~l~
  S(m) =S1(l);  %M~tR
end
fid=fopen(fno,'w');
fprintf (fid,'Freq(Hz)   Rate (MM)   Coef\n');
for k=1: m
  fprintf (fid,'%10.4f  %10.4f  %10.6f  %10.4f  %10.4f\n', F(k), D(k)*100.0, imag(S(k)), 2*abs(S(k)), angle(S(k))/rad);
end
status=fclose (fid);