Hcanonical(mass,damp,stiff,Wmax,N_Sample,Noise_Factor,out,inp)

function [H,w,modal_par] = Hcanonical(mass,damp,stiff,Wmax,N_Sample,Noise_Factor,out,inp)

%HCANONICAL   Frequency response function and modal parameter estimation for a N 
% degrees of freedom lineal forced system with viscous damping.
% 
% Syntax: [H,w,par] = Hcanonical(mass,damp,stiff,Wmax,N_Sample,Noise_Factor,out,inp)
%
% Example for system with 3 degrees of freedom:
%       mass         = 0.3*[1,0,0;0,1,0;0,0,1];       %mass matrix
%       damp         = 0.1*[2,-1,0;-1,2,-1;0,-1,2];   %damp matrix
%       stiff        = 10000*[2,-1,0;-1,2,-1;0,-1,2]; %stiff matrix
%       Wmax         = 500;                           %max. frequency (rad/sec)
%       N_Sample     = 512;                           %number of frequency samples
%       Noise_Factor = 0;                             %noise factor
%       out          = 1;                             %response point
%       inp          = 1;                             %input point
%
%       H   = Frequency response function
%       w   = frequency range (rad/sec)
%       modal_par = Modal Parameters [freq,damp,Ci,Oi]: 
%                   freq = Natural frequencies (rad/sec)
%                   damp = Damping ratio
%                   Ci   = Amplitude modal constant
%                   Oi   = Phase modal constant (degrees)
% 
%-----------------------------------------------------------------------------------
% Cristian Gutierrez Acuna, February 2002. [e-mail]
%-----------------------------------------------------------------------------------

dw = Wmax/N_Sample;
[N,c] = size(mass);
null = zeros(N,N);

%state-space formulation
A = [damp,mass;mass,null];
B = [stiff,null;null,-mass];

[V,D] = eig(B,-A); 
	lambda = diag(D);
	[Y,I] = sort(imag(lambda));
	polos = lambda(I);
    psi = V(:,I);
a = psi.'*A*psi;
b = psi.'*B*psi;

%Frequency Response Function
FRF = 0;
for n = 1:N
   w = dw:dw:Wmax;
   residuo = psi(1:N,N+n)*psi(1:N,N+n).'./a(N+n,N+n);
   H = residuo(out,inp)./(j.*w - polos(N+n)) + conj(residuo(out,inp))./(j.*w - conj(polos(N+n)));
   FRF = FRF + H;   
   matrix_residuo(n,1) = residuo(out,inp);
end

%generating IRF. If Noise_Factor = 0, this part can be omitted
IRF = ifft(FRF,N_Sample);
max_irf = max(abs(IRF));
ruido = Noise_Factor*max_irf*randn(size(IRF));
%sprintf('SNR = %0.5g [dB].',20*log10(std(IRF)/std(ruido))),
IRF = IRF + ruido;
H = fft(IRF,N_Sample);

polos = polos((length(polos)/2)+1:end,1); 
	freq = abs(polos); %Natural frequencies (rad/sec)
    damp = -real(polos)./abs(polos); %Damping ratio
   
Ai = -2*(real(matrix_residuo).*real(polos) + imag(matrix_residuo).*imag(polos));
Bi = 2*real(matrix_residuo);
const_modal = complex(Ai,abs(polos).*Bi);
	Ci = abs(const_modal); %Amplitude modal constant
	Oi = angle(const_modal).*180/pi; %Phase modal constant (degrees)

modal_par = [freq, damp, Ci, Oi] %Modal Parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%