% Script to run the Monte Carlo simulation* associated with
% the performance results of the SCCM, reported in Section VIII.B
% (Table III - Case 2 and Table VI) 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.....

% *Matlab script to run the SCCM performance evaluation
% with regards step Characterization of correctly classified
% Type 2 (phase steps). Since SCCM does an almost perfect
% step type classification, here all AC signals have a phase
% step. 

% Tested in Matlab 2012b and R2023a.

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

% @ Paulo Esquef - March/24 


clear all; close all; clc;

SNR = 30;  % signal-to-noise ratio
% change the SNR above to 55 or 40 to generate the other 2 conditions
% replace also the SNR value in the name of the .mat to be saved (below)
% save sccm_eval_case2_30dB .... 
F0 = 60.0; %nominal frequency
F1 = 60; %fundamental frequency of x[n]
Xm=1; % nominal magnitude of the AC signal 

Fs=4800;  % sampling frequency 
NCycles = 6;  % number of 60 Hz cycles
T = NCycles/F0;  % window duration in seconds
N = floor(NCycles*Fs/F0); % window duration in samples

n=0:N-1; n=n(:); % time basis in samples

% SCCM setup
% PATV parameters for the SLHEM
dm=0;  % polynomial order
lambda_m=2;  % lambda_m
Nit=100;  % number of iterations in PATV filtering
% detection threshold
Lm=1e-5; % detection threshold


% PATV parameters for the SLHEP 
dp=1;  % polynomial order
lambda_p=2;   % lambda_p
lambda_v=20e-3; % lambda_v (for the PATV-based estimation of the linear phase)
Lp=1e-5; % detection threshold

MC=20000;  % number of Monte Carlo runs

tau_v=randi([60 417],MC,1);
% true time instants in samples of the magnitude step, from tau = 60 samples
% to tau = 417 samples, drawn from an integer Uniform distribution

type_step_true=2; % true step type: always a phase step (Case 2)
% type_step_true = 0 implies no step is present in the AC signal
% type_step_true = 1 implies a mag step is present in the AC signal
% type_step_true = 2 implies a phase step is present in the AC signal


epsilon_tau=zeros(MC,1);   % vector to store estimation errors for the step location
error_hm=zeros(MC,1); % vector to store estimation errors for the step height
hp_est_v=zeros(MC,1); % vector to store step height height estimates

type_step_est_v=zeros(MC,1); % vector to store the estimated step type

tau_est_v=zeros(MC,1);  % vector to store the estimated values of the step location

% range of values (positive and negative as well) for hp
hp1=pi/36; hp2=5*pi/18;  % range of values of the phase step depth (in rad)
hp_ref_v=(hp1+(hp2-hp1).*rand(MC,1)); %   pi/36  <= hp <= 5*pi/18 rad
sign_switch=2*(randi([0 1],MC,1)-.5); % uniform switch process with values -1 and 1
hp_ref_v=hp_ref_v.*sign_switch;  % true values of hp (positive and negative) 

Nc=0;  % counter for step type misclassification. This is used
% in the performance metric \varepsilon_f(m_tau)

for ll=1:MC,
    hv= [0; hp_ref_v(ll); 0]; % phase step from pi/36 rad to 5*pi/18 rad (and negative values)
    h=hv; % uses prescribed step height
    tau=tau_v(ll);  % selects a given step location in the range
    phi_0_rad=randi([0,180]).*pi./180;  % selects initial phase randomly 
    % from Uniform distribution from 0 to 180 degrees and converts to rad
    % generates a noisy AC signal with a step
    w1 = 2*pi*F1/Fs; w2 = 2*pi*(h(3))/Fs; % instantaneous frequency
    Xmn = Xm*ones(size(n))+Xm.*(h(1)).*(n>=tau); % instantaneous magnitude
    PHI=w1*n+w2.*(n-tau).*(n>=tau)+(h(2)).*(n>=tau)+ phi_0_rad; % instantaneous phase
    x=Xmn.*cos(PHI);  % deterministic AC signal with a step
    vx=var(x); vruido=vx./(10^(SNR/10));  % variance of the noise component
    xn=x+sqrt(vruido)*randn(size(n));  % noisy AC signal
    % runs SCCM to obtain: step type and related estimates of the step location and height
    [tau_est,type_step_est,hm_est,hp_est]=sccm_ref(xn,dm,dp,lambda_m,lambda_p,lambda_v,Nit,Lm,Lp);
    clear xn w1 w2 Xmn PHI
    if (type_step_est ~= 2), % if a step type misclassification happens
        Nc=Nc+1;  % increments Nc (counter for misclassified steps)
        % Nc takes into account false negatives (estimated step type =0) and
        % wrong classification (estimated step type = 2, i.e., phase step)
        hp_est=NaN;  % invalid step height estimate
        tau_est=NaN; % invalid step location estimate
    end
        
    hp_est_v(ll)=hp_est;  % stores step height estimates
    epsilon_tau(ll)=tau_est-tau; % stores step location errors  
    error_hp(ll)=hp_est-hp_ref_v(ll);% stores step height errors
    tau_est_v(ll)=tau_est; % stores estimated step locations
end

% Uncomment the lines (one at a time) below to save the results 
%save sccm_eval_case2_55dB hm_ref_v hm_est_v error_hm tau_v tau_est_v epsilon_tau
%save sccm_eval_case2_40dB hm_ref_v hm_est_v error_hm tau_v tau_est_v epsilon_tau
%save sccm_eval_case2_30dB hm_ref_v hm_est_v error_hm tau_v tau_est_v epsilon_tau

% Performance metric for step location \varepsilon_f(m_\tau)

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%



% Performance metric for step height estimation

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


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

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

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

