% Script to reproduce the results shown in Table IV (Case 1) and Table VII
% for hm=0.05 (minimum),  
% of the manuscript "Classification and Characterization of a Step in 
% the Magnitude or Phase of an AC Signal", submitted for publication as 
% research article to Journal.....

% Data has been produce by running the
% script SCCM_evaluation_step_characterization_magnitude_step_hm_min.m
% to evaluate the SCCM performance 
% with regards location charactezization of magnitude steps


% Tested for Matlab 2012b and R2023a.

% Script written by Paulo Esquef: pesquef@lncc.br
% Last modified: May 3rd, 2024

clear all; close all; clc

load sccm_eval_case1_hm_min_55dB
% loads vectors hm_ref_v hm_est_v error_hm tau_v tau_est_v epsilon_tau

MC=length(tau_est_v);
Nc=sum(isnan(tau_est_v));


aux00=isnan(tau_est_v);
epsilon_tau(aux00)=[];  % removes undetected (misclassified) mag steps, i.e., tau_est=NaN.
% These errors (step type misclassification) will be taken into account in count_Nc
%Nc=count_Nc;  % Nc takes into account false negatives (estimated step type =0) and
% wrong classification (estimated step type = 2, i.e., phase step)

range_mtau=[0:1:80];  % range of values for m_tau to evaluate
% metric function \varepsilon_f(m_tau)
Q=length(range_mtau);

varepsilon_f_v=zeros(Q,1);  % vector to store \varepislon_f(m_tau)

for k=1:Q
    m_tau=range_mtau(k); 
    Nu=sum(abs(epsilon_tau)> m_tau); % number of times the estimated location
    % of a detected step is above the error tolerance m_\tau
    varepsilon_f = 100 * (Nc+Nu)./MC; % calculates \varepislon_f(m_tau)
    varepsilon_f_v(k)=varepsilon_f; % stores calculated value
end

% Plot the results
stem(range_mtau,varepsilon_f_v,'linewidth',2)
xlabel('$m_\tau$ (in samples)','interpreter','latex');
ylabel('$\varepsilon_{\textrm{f}}(m_\tau)$ (in \%)','interpreter','latex')
hold on
plot(range_mtau,4.55*ones(size(varepsilon_f_v)),'r','linewidth',2) % plots threshold of 4.55%

aux=find(varepsilon_f_v<4.55); ind=aux(1); 
m_tau_min=range_mtau(ind);
varepsilon_f=varepsilon_f_v(ind);
disp('===================================')
disp('Report Results of Table IV (Case 1), for 55 dB:')

m_tau_min

disp(['varepsilon_f(', num2str(m_tau_min),') = ',num2str(varepsilon_f),' %'])

disp('===================================')
disp('Report results for 55 dB and hm=0.05 (not reported in the paper)')

% Report results of Table V
% Performance metric for step height estimation

% This is only for valid steps (correctly classified) 
% finds height invalid estimates (NaN) and remove them 
aux11=isnan(hm_est_v);
hm_est_v(aux11)=[];
hm_ref_v(aux11)=[];
error_hm(aux11)=[];


disp('Mean relative error (MRE) of the step height estimates (in %):')
MRE=mean(100*(hm_est_v-hm_ref_v)./hm_ref_v)

disp('Mean \mu of the error distribution (height estimate):')
mu=mean(error_hm)

disp('Standard deviation \sigma of the error distribution (height estimate):')
sigma=std(error_hm)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all; 

load sccm_eval_case1_hm_min_40dB
% loads vectors hm_ref_v hm_est_v error_hm tau_v tau_est_v epsilon_tau

MC=length(tau_est_v);
Nc=sum(isnan(tau_est_v));


aux00=isnan(tau_est_v);
epsilon_tau(aux00)=[];  % removes undetected (misclassified) mag steps, i.e., tau_est=NaN.
% These errors (step type misclassification) will be taken into account in count_Nc
%Nc=count_Nc;  % Nc takes into account false negatives (estimated step type =0) and
% wrong classification (estimated step type = 2, i.e., phase step)

range_mtau=[0:1:80];  % range of values for m_tau to evaluate
% metric function \varepsilon_f(m_tau)
Q=length(range_mtau);

varepsilon_f_v=zeros(Q,1);  % vector to store \varepislon_f(m_tau)

for k=1:Q
    m_tau=range_mtau(k); 
    Nu=sum(abs(epsilon_tau)> m_tau); % number of times the estimated location
    % of a detected step is above the error tolerance m_\tau
    varepsilon_f = 100 * (Nc+Nu)./MC; % calculates \varepislon_f(m_tau)
    varepsilon_f_v(k)=varepsilon_f; % stores calculated value
end

% Plot the results
figure
stem(range_mtau,varepsilon_f_v,'linewidth',2)
xlabel('$m_\tau$ (in samples)','interpreter','latex');
ylabel('$\varepsilon_{\textrm{f}}(m_\tau)$ (in \%)','interpreter','latex')
hold on
plot(range_mtau,4.55*ones(size(varepsilon_f_v)),'r','linewidth',2) % plots threshold of 4.55%

aux=find(varepsilon_f_v<4.55); ind=aux(1); 
m_tau_min=range_mtau(ind);
varepsilon_f=varepsilon_f_v(ind);
disp('===================================')
disp('Report Results of Table IV, for 40 dB and hm=0.05:')

m_tau_min

disp(['varepsilon_f(', num2str(m_tau_min),') = ',num2str(varepsilon_f),' %'])

disp('===================================')
disp('Report results of for 40 dB and hm=0.05 (not reported in the paper) ')

% Report results of Table V
% Performance metric for step height estimation

% This is only for valid steps (correctly classified) 
% finds height invalid estimates (NaN) and remove them 
aux11=isnan(hm_est_v);
hm_est_v(aux11)=[];
hm_ref_v(aux11)=[];
error_hm(aux11)=[];


disp('Mean relative error (MRE) of the step height estimates (in %):')
MRE=mean(100*(hm_est_v-hm_ref_v)./hm_ref_v)

disp('Mean \mu of the error distribution (height estimate):')
mu=mean(error_hm)

disp('Standard deviation \sigma of the error distribution (height estimate):')
sigma=std(error_hm)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all; 
load sccm_eval_case1_hm_min_30dB
% loads vectors hm_ref_v hm_est_v error_hm tau_v tau_est_v epsilon_tau

MC=length(tau_est_v);
Nc=sum(isnan(tau_est_v));


aux00=isnan(tau_est_v);
epsilon_tau(aux00)=[];  % removes undetected (misclassified) mag steps, i.e., tau_est=NaN.
% These errors (step type misclassification) will be taken into account in count_Nc
%Nc=count_Nc;  % Nc takes into account false negatives (estimated step type =0) and
% wrong classification (estimated step type = 2, i.e., phase step)

range_mtau=[0:1:80];  % range of values for m_tau to evaluate
% metric function \varepsilon_f(m_tau)
Q=length(range_mtau);

varepsilon_f_v=zeros(Q,1);  % vector to store \varepislon_f(m_tau)

for k=1:Q
    m_tau=range_mtau(k); 
    Nu=sum(abs(epsilon_tau)> m_tau); % number of times the estimated location
    % of a detected step is above the error tolerance m_\tau
    varepsilon_f = 100 * (Nc+Nu)./MC; % calculates \varepislon_f(m_tau)
    varepsilon_f_v(k)=varepsilon_f; % stores calculated value
end

% Plot the results
figure
stem(range_mtau,varepsilon_f_v,'linewidth',2)
xlabel('$m_\tau$ (in samples)','interpreter','latex');
ylabel('$\varepsilon_{\textrm{f}}(m_\tau)$ (in \%)','interpreter','latex')
hold on
plot(range_mtau,4.55*ones(size(varepsilon_f_v)),'r','linewidth',2) % plots threshold of 4.55%

aux=find(varepsilon_f_v<4.55); ind=aux(1); 
m_tau_min=range_mtau(ind);
varepsilon_f=varepsilon_f_v(ind);
disp('===================================')
disp('Report Results of Table IV, for 30 dB and hm=0.05:')

m_tau_min

disp(['varepsilon_f(', num2str(m_tau_min),') = ',num2str(varepsilon_f),' %'])

disp('===================================')
disp('Report results of Table VII, for 30 dB and hm=0.05: ')

% Report results of Table V
% Performance metric for step height estimation

% This is only for valid steps (correctly classified) 
% finds height invalid estimates (NaN) and remove them 
aux11=isnan(hm_est_v);
hm_est_v(aux11)=[];
hm_ref_v(aux11)=[];
error_hm(aux11)=[];


disp('Mean relative error (MRE) of the step height estimates (in %):')
MRE=mean(100*(hm_est_v-hm_ref_v)./hm_ref_v)

disp('Mean \mu of the error distribution (height estimate):')
mu=mean(error_hm)

disp('Standard deviation \sigma of the error distribution (height estimate):')
sigma=std(error_hm)



