% Script to run the Monte Carlo simulation* associated with
% the performance results of the SCCM, reported in Section VIII.B
% (Table IV - Case 1 and Table VII, for hm=0.05) 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 1 (magnitude steps), with minimum step change hm=0.05.
% Since SCCM does an almost perfect
% step type classification, here all AC signals have a magnitude
% step of hm=0.05. 

% 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_case1_hm_min_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=1; % true step type: always a magnitude step (Case 1)
% 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
hm_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

% value (positive and negative as well). Here, hm = \pm 0.05.
hm1=0.05; hm2=0.05;  % range of relative values of the magnitude step
hm_ref_v=(hm1+(hm2-hm1).*rand(MC,1)); % .05 <=  hm <= .05 
sign_switch=2*(randi([0 1],MC,1)-.5); % uniform switch process with values -1 and 1
hm_ref_v=hm_ref_v.*sign_switch;  % true values of hm (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= [hm_ref_v(ll); 0; 0]; % relative magnitde step change
    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 ~= 1), % 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)
        hm_est=NaN;  % invalid step height estimate
        tau_est=NaN; % invalid step location estimate
    end
        
    hm_est_v(ll)=hm_est;  % stores step height estimates
    % For Xm not equal to 1, replace line above by: hm_est_v(ll)=hm_est/Xm;
    epsilon_tau(ll)=tau_est-tau; % stores step location errors  
    error_hm(ll)=hm_est-hm_ref_v(ll);% stores step height errors
    tau_est_v(ll)=tau_est; % stores estimated step locations
end

% Uncomment the lines below to save the results 
%save sccm_eval_case1_hm_min_55dB hm_ref_v hm_est_v error_hm tau_v tau_est_v epsilon_tau
%save sccm_eval_case1_hm_min_40dB hm_ref_v hm_est_v error_hm tau_v tau_est_v epsilon_tau
%save sccm_eval_case1_hm_min_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(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 %):')
mean(100*(hm_est_v-hm_ref_v)./hm_ref_v)

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

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

