clc %emol
clearvars %emol
close all %emol

load impact_data_raw
Navg=10; % for unprocessed data set (emol)

fs=256;
N=length(in1); 
f=fs*(0:N-1)/N;
T=N/fs; 
t=0:1/fs:T-1/fs; 

for n=1:Navg
  In = ['In', int2str(n),' = fft(in', int2str(n), ');'];
  eval(In);
  Out = ['Out', int2str(n),' = fft(out', int2str(n), ');'];
  eval(Out);
  Sxx = ['Sxx', int2str(n),' = conj(In',int2str(n),')' '.*In',int2str(n), ';'];
  eval(Sxx);
  Sxy = ['Sxy', int2str(n),' = conj(In',int2str(n),')' '.*Out',int2str(n), ';'];
  eval(Sxy);
end

Sxx=[]; Sxy=[];  
for n=1:Navg
  tmp1= ['Sxx', int2str(n), ';'];
  Sxx=[Sxx; eval(tmp1)];
  tmp2= ['Sxy', int2str(n), ';'];
  Sxy=[Sxy; eval(tmp2)];
end
Sxx=mean(Sxx); Sxy=mean(Sxy); 
H1=Sxy./Sxx; 

figure(1)
subplot(2,2,1) %emol
plot(f(41:761),20*log10(abs(H1(41:761))))
axis([5 95 -30 50])
xlabel('Frequency (Hz)'); ylabel('FRF (Modulus, dB)')
grid on %emol
subplot(2,2,2) %emol
plot(f(41:761),unwrap(angle(H1(41:761))))
axis([5 95 -3.5 3.5])
xlabel('Frequency (Hz)'); ylabel('FRF (Phase, rad)')
grid on %emol

% % repetion of same procedure used in case A (emol)
% load impact_data_pre_processed
% Navg=3; % for pre-processed data set
% 
% for n=1:Navg
%   In = ['In', int2str(n),' = fft(in', int2str(n), ');'];
%   eval(In);
%   Out = ['Out', int2str(n),' = fft(out', int2str(n), ');'];
%   eval(Out);
%   Sxx = ['Sxx', int2str(n),' = conj(In',int2str(n),')' '.*In',int2str(n), ';'];
%   eval(Sxx);
%   Sxy = ['Sxy', int2str(n),' = conj(In',int2str(n),')' '.*Out',int2str(n), ';'];
%   eval(Sxy);
% end
% 
% Sxx=[]; Sxy=[];  
% for n=1:Navg
%   tmp1= ['Sxx', int2str(n), ';'];
%   Sxx=[Sxx; eval(tmp1)];
%   tmp2= ['Sxy', int2str(n), ';'];
%   Sxy=[Sxy; eval(tmp2)];
% end
% Sxx=mean(Sxx); Sxy=mean(Sxy); 
% H1=Sxy./Sxx; 
% 
% figure(1)
% subplot(2,2,3) %emol
% plot(f(41:761),20*log10(abs(H1(41:761))))
% axis([5 95 -30 50])
% xlabel('Frequency (Hz)'); ylabel('FRF (Modulus, dB)')
% grid on %emol
% subplot(2,2,4) %emol
% plot(f(41:761),unwrap(angle(H1(41:761))))
% axis([5 95 -3.5 3.5])
% xlabel('Frequency (Hz)'); ylabel('FRF (Phase, rad)')
% grid on %emol
% 
% % --> Averaging improves the FRF estimate (its effect is to remove the
% %     noises that are 'uncorrelated' with the input and output signals);
% % --> Preprocessing results in a much better FRF estimate using far 
% %     fewer averages