function [FRF,fi] = ComputeFRFs(time,U,Y,nens,plotting);
%
% [FRF,fi] = ComputeFRFs(time,U,Y,nens,plotting);
%
% This function calculates RAOs based on averaging multiple output files
% generated with white noise excitation.  Averaging four files appears to give decent
% results.  Six files gives good results.
%
% Inputs:   time = time signal
%           U = input signal
%           Y = output signal
%           nens = number of ensembles to average (i.e., 4 or 6)
%           plotting = flag for whether you want the program to generate
%           plots for you.
%
% Outputs:  FRF = frequency response function for each of the outputs
%           desired (same as the RAO)
%           fi = associated frequency vector
%
% -------------------------------------------------------------------------
%% Get the time information about the simulation
N = length(time);
dt = time(2)-time(1);
time_span = time(end)-time(1);
% -------------------------------------------------------------------------
%% Compute RAOs- using ensemble averaging
% Compute ensemble length and new time/frequency vectors for each ensemable
ni = floor(N/nens);         % ensemble length 
time_spani = time(ni)-time(1);
dfi = 1/time_spani;
fi = [0:ni-1]/(ni-1)*(1/dt); fi = fi(1:floor(end/2));
% Initialize auto and cross-spectral densities
UU_psd = zeros(ni,1);
YU_psd = zeros(ni,1);
% Auto spectral density - U*U
for i=1:nens
    U_fft = fft(U(ni*(i-1)+1:ni*i));
    UU_psd =  UU_psd + 2*abs( U_fft.*conj(U_fft) )./(ni*ni*dfi);
end;
% UU_psd = UU_psd./nens;
UU_psd = UU_psd(floor(1:ni/2));  %Take only the 1-sided spectrum
    
% Cross spectral density - Y*U
for i=1:nens
    U_fft = fft(U(ni*(i-1)+1:ni*i));
    Y_fft = fft(Y(ni*(i-1)+1:ni*i));
    YU_psd = YU_psd + 2*abs( Y_fft.*conj(U_fft) )./(ni*ni*dfi);
end;
    
% YU_psd = YU_psd./nens;
YU_psd = YU_psd(floor(1:ni/2));  %Take only the 1-sided spectrum    
 
% FRF - YU/UU = Y/U
FRF = YU_psd./UU_psd; % - baseline(j);   
 
 
% -------------------------------------------------------------------------
% Plot the results
if plotting
    figure;
    semilogy( fi,FRF);
    xlabel('Frequency (Hz)'); 
    ylabel(outputs); 
    grid on;
    axis('tight');
end;
% if plotting
%     figure;
%     for j=1:no
%         subplot(3,2,j);
%         semilogy( fi,FRF{j});
%         xlabel('Frequency (Hz)'); 
%         ylabel(outputs{j}); 
%         grid on;
%         axis('tight');
%     end;
% end;