function [ Fw_mag, Fw_phase, F0_mag, F0_phase, f] = IdentSpectral( MeasurementData, ModeValues, ModeAnalysis, length_window_barlet )
%IdentSpectral - function for non-parametric spectral identification method
%   [INPUT]
%   MeasurementData - class of type MotionIdentTwinCAT with input/output
%                     system data
%   ModeValues - mode of values for computation:
%                       - 'InterestingValues' - compute transfer from values that was measured at the time of sampling PRBS
%                       - 'AllValues' - compute transfer from values that was measured at the time of sampling with cycle time
%   ModeAnalysis - mode of spectral analys:
%                       - 'Basic' - basic
%                       - 'Barlett' - Barlett method with Hann window
%   [OUTPUT]
%   Fw_mag - magnitude of closed loop
%   Fw_phase - phase of closed loop [degrees - °]
%   F0_mag - magnitude of open loop
%   F0_phase - phase of open loop [degrees - °]
%   f - frequency for x label [Hz]

NumberRepeat = size(MeasurementData.CommandValue,1);

% % 1 step delay
% MeasurementData.CommandValueInterest = MeasurementData.CommandValueInterest(1:end-1);
% MeasurementData.FeedbackValueInterest = MeasurementData.FeedbackValueInterest(2:end);

% % 2 steps delay
% MeasurementData.CommandValueInterest = MeasurementData.CommandValueInterest(1:end-2);
% MeasurementData.FeedbackValueInterest = MeasurementData.FeedbackValueInterest(3:end);

% Init variables
if strcmp(ModeValues,'InterestingValues')  % Interesting data - at the time of sampling PRBS 
    if strcmp(ModeAnalysis,'Basic')
         Fw = complex(zeros(NumberRepeat,floor(size(MeasurementData.CommandValueInterest,2)/2)));
     elseif strcmp(ModeAnalysis,'Barlett')
         Fw = complex(zeros(NumberRepeat,floor(length_window_barlet)));
     else
         error('Wrong method analysis');
    end        
elseif strcmp(ModeValues,'AllValues')  % ALL data - at the time of sampling with cycle time
    if strcmp(ModeAnalysis,'Basic')
         Fw = complex(zeros(NumberRepeat,floor(size(MeasurementData.CommandValue,2)/2)));
     elseif strcmp(ModeAnalysis,'Barlett')
         Fw = complex(zeros(NumberRepeat,floor(length_window_barlet)));
     else
         error('Wrong method analysis');
    end        
else
        error('Wrong mode of values for computation');
end
Fw_avg_mag = zeros(size(Fw));
Fw_avg_phase = zeros(size(Fw));
F0_avg_mag = zeros(size(Fw));
F0_avg_phase = zeros(size(Fw));
    
for i = 1:NumberRepeat
    
    if strcmp(ModeValues,'InterestingValues')  % Interesting data - at the time of sampling PRBS          
        CommandValue = MeasurementData.CommandValueInterest(i,:);
        FeedbackValue = MeasurementData.FeedbackValueInterest(i,:);
        
    elseif strcmp(ModeValues,'AllValues')  % ALL data - at the time of sampling with cycle time
        CommandValue = MeasurementData.CommandValue(i,:);
        FeedbackValue = MeasurementData.FeedbackValue(i,:);
        
    else
        error('Wrong mode of values for computation');
    end
    
    if MeasurementData.AmplAuxMove ~= 0 
        MeanCommand = mean(CommandValue);
        CommandValue = CommandValue - MeanCommand;
               
        FeedbackValue = FeedbackValue - MeanCommand;
%         MeanFeedBack = mean(FeedbackValue);
%         FeedbackValue = FeedbackValue - MeanFeedBack;      
    end

     % Computed transfer
     if strcmp(ModeAnalysis,'Basic')
         [ Fw_i] = spa_0( CommandValue, FeedbackValue);
         
          % throw away:
          %               -    half of spectrum Fw_i -> discreete system, repeated
          %                    spectrum
          Fw(i,:) = Fw_i(1:floor(length(Fw_i)/2));
          
         
         % Magnitude and phase for closed loop Fw
         Fw_avg_mag(i,:) = abs(Fw(i,:));       % Magnitude
         Fw_avg_phase(i,:) = phaseE(Fw(i,:));    % Phase
    
     elseif strcmp(ModeAnalysis,'Barlett')
         [ Fw_i] = spa_b( CommandValue, FeedbackValue, length_window_barlet);
         
          % throw away:
          %               -    half of spectrum Fw_i -> discreete system, repeated
          %                    spectrum
          Fw(i,:) = Fw_i(1:ceil(length(Fw_i)/2));
         
         % Magnitude and phase for closed loop Fw
         Fw_avg_mag(i,:) = abs(Fw(i,:));       % Magnitude
         Fw_avg_phase(i,:) = phaseM(Fw(i,:));    % Phase
        
     else
         error('Wrong method analysis');
     end
    
    % compute F0(open loop) from Fw(closed loop)
    if strcmp(ModeAnalysis,'Basic')
        [ F0_avg_mag(i,:), F0_avg_phase(i,:)] = FwtoF0( Fw_avg_mag(i,:), Fw_avg_phase(i,:), 'Base' );
    else
        [ F0_avg_mag(i,:), F0_avg_phase(i,:)] = FwtoF0( Fw_avg_mag(i,:), Fw_avg_phase(i,:));
    end
    
% %     Filtering
% %     non-linear
%     Fw_avg_mag(i,:) = medfilt1(Fw_avg_mag(i,:),3);
%     Fw_avg_phase(i,:) = medfilt1(Fw_avg_phase(i,:),3);
%     F0_avg_mag(i,:) = medfilt1(F0_avg_mag(i,:),3);
%     F0_avg_phase(i,:) = medfilt1(F0_avg_phase(i,:),3);
end

% Init variables
Fw_mag = zeros(1,size(Fw_avg_mag,2));
Fw_phase = zeros(1,size(Fw_avg_mag,2));
F0_mag = zeros(1,size(Fw_avg_mag,2));
F0_phase = zeros(1,size(Fw_avg_mag,2));

% % Averaging with median
% for i = 1:size(Fw_avg_mag,2)
%     Fw_mag(1,i) = median(Fw_avg_mag(:,i));
%     Fw_phase(1,i) = median(Fw_avg_phase(:,i));
%     
%     F0_mag(1,i) = median(F0_avg_mag(:,i));
%     F0_phase(1,i) = median(F0_avg_phase(:,i));
% end

% Averaging with mean
for i = 1:size(Fw_avg_mag,2)
    Fw_mag(1,i) = mean(Fw_avg_mag(:,i));
    Fw_phase(1,i) = mean(Fw_avg_phase(:,i));
    
    F0_mag(1,i) = mean(F0_avg_mag(:,i));
    F0_phase(1,i) = mean(F0_avg_phase(:,i));
end

% Convert phase from radians to degrees
Fw_phase = Fw_phase.*(180/pi);
F0_phase = F0_phase.*(180/pi);

% Frequency for x label in Hz - x axis
if strcmp(ModeValues,'InterestingValues')    
    f = (1:length(Fw))/MeasurementData.PeriodInputSignal/size(Fw_i,2);
elseif strcmp(ModeValues,'AllValues')
    f = (1:length(Fw))/MeasurementData.BaseTime/size(Fw_i,2);
else
    error('Wrong mode of values for computation');
end

end

