%----- example matlab m-file for Chapter 5 -----

load nsimR05  % This simulation data is for the 'input noise' case
              %   with RMS 5%  noise.
              % Another data 'nsimR20' is for the 'input noise' case
              %   with RMS 20% noise.
M=4100/980.;  % total weight of the model.  see Fig. 5.1

dir='NS';   % Choose 'NS' or 'EW'
if dir=='NS',  chZ=1; chX=2;  else  chZ=7; chX=8;  end

n1=4096;    % n1 : segment length for FFT
            % You may change this value.
n2=length(eq);
ncase=1;    % ncase = 1 : input noise case
            %       = 2 : output noise case

freq=(fs/2)*(0:n1/2)/(n1/2);
omg=2*pi*freq;

fmax=500;  % fmax : the index for the highest frequency to calculate values
kCn = zeros(1,fmax);

if ncase==1    % input-noise case
   Pxz=csd(eq(1:n2,chX),eq(1:n2,chZ),n1,fs,hanning(n1),n1/2);
   Pxx=psd(eq(1:n2,chX),n1,fs,hanning(n1),n1/2);
   for L=1:fmax
      kCn(L)=M*(-omg(L)^2)*(Pxx(L)/(Pxz(L)-Pxx(L)));  % see eq. (5.16)
   end
else           % output-noise case
   Pzz=psd(eq(1:n2,chZ),n1,fs,hanning(n1),n1/2);
   Pzx=csd(eq(1:n2,chZ),eq(1:n2,chX),n1,fs,hanning(n1),n1/2);
   for L=1:fmax,
      kCn(L)=M*(-omg(L)^2)*(Pzx(L)/(Pzz(L)-Pzx(L)));  % see eq. (5.15)
   end
end

kn=real(kCn);

%----- least squares estimation using a linear function for Kn
fm1=0.3; fm2=0.6;  % fm1 : the lower frequency for the calculation
                   % fm2 : the upper frequency for the calculation
[tmp,ifm1]=min(abs(freq-fm1));
[tmp,ifm2]=min(abs(freq-fm2));
Zk=kn(ifm1:ifm2)';
Hf=[freq(ifm1:ifm2)',ones(ifm2-ifm1+1,1)];   %  y=ax+b
% rem : As shown in the text book, you may use a parabolic function
%       for least squares estimation with the following line. 
% Hf=[(freq(ifm1:ifm2).^2)',freq(ifm1:ifm2)',ones(ifm2-ifm1+1,1)];
                                             %  y=ax^2+bx+c      
Pk=(Hf'*Hf)\(Hf'*Zk);
kne=Pk(2);    % identified Kn

%----- plotting -- stiffness
figure(1)
subplot(3,1,1)
plot(freq(1:fmax),kn,'r-')   % The upper figure of Fig. 5.3
xmax=2.0;
ymax=40;   % Max value for the Y-axis.
axis([0 xmax 0 ymax])
xlabel('frequency (Hz)')
ylabel('stiffness (tonf/cm)')
text(xmax*0.1,ymax*0.85,['identified Kn =',num2str(kne,3)])
textt=[date,' : ',type,' : ',chname(chX,:),' / ',chname(chZ,:)];
title(textt)
grid

%----- calculating h
hn=0.5*imag(kCn)/kne;

fm1=0.3; fm2=0.6;
[tmp,ifm1]=min(abs(freq-fm1));
[tmp,ifm2]=min(abs(freq-fm2));
Zh=hn(ifm1:ifm2)';
Hf=[freq(ifm1:ifm2)',ones(ifm2-ifm1+1,1)];   %  y=ax+b
Ph=(Hf'*Hf)\Hf'*Zh;
beta=Ph(2);    % identified beta

cn=imag(kCn);
Zc=cn(ifm1:ifm2)';
Pc=(Hf'*Hf)\Hf'*Zc;
Cne=Pc(1)/(2*pi);     % identified c

%----- plotting -- equivalent-h & viscous damping
subplot(3,1,2)
plot(freq(1:fmax),hn,'r-')   % The upper figure of Fig. 5.5
xmax=2.0;
ymax=1.0;  % Max value for the Y-axis.
axis([0 xmax 0 ymax])
xlabel('frequency (Hz)')
ylabel('hysteretic damping ratio')
text(xmax*0.11,ymax*0.85,['identified beta = ',num2str(beta,3)])
grid

subplot(3,1,3)
plot(freq(1:fmax),cn,'r-')   % The upper figure of Fig. 5.5
ymax=40;   % Max value for the Y-axis.
axis([0 xmax 0 ymax])
xlabel('frequency (Hz)')
ylabel('imaginary part (tonf/cm)')
text(xmax*0.11,ymax*0.85,['identified c =',num2str(Cne,3)])
grid
