function [sqiEcg, error_code] = sqi(ecg, R_loc, fs, beatWidth)
% sqi estimates the quality of an ECG segment
%
% Author: Mohamed Abdelazez
%
% Inputs
%   ecg: column vector containing an ECG segment.
%   ann: column vector containing the location of the R wave in samples
%   fs: scalar containing the sampling frequnecy 
%   beatWidth: scalar containing the duration of a single beat in sec (e.x. 0.7 sec)
%
% Outputs
%   sqiEcg: scalar value representing the quality of an ECG segment
%   error_code: scalar value representing the error encountered (0 - No
%   Error, 1 - less than three beats provided)
%
% The function will return NaN if it encounters an error
%
% Please reference the following papers when using the function:
%   Patrick X. Quesnel, Adrian D. C. Chan, Homer Yang, 'Signal quality and false 
%   myocardial ischemia alarms in ambulatory electrocardiograms.', 2014 IEEE International 
%   Symposium on Medical Measurements and Applications (MeMeA), Lisbon, Portugal, 2014, pp. 1-5
%
%   Abdelazez, Mohamed, et al. "Signal quality analysis of ambulatory electrocardiograms 
%   to gate false myocardial ischemia alarms." IEEE Transactions on Biomedical Engineering 
%   64.6 (2016): 1318-1325.
%
% Modifications
% 14-Apr-2014 PQ created original SQI functionality
% 15-Dec-2019 MA integrated all functionality in a single function

if isrow(ecg)
   ecg = ecg'; 
end

if isrow(R_loc)
    R_loc = R_loc';
end

beatWidthSamp = floor(beatWidth * fs);
beatWidthSampHalf = floor(beatWidthSamp / 2);
nbeats = length(R_loc);

if R_loc(1) - beatWidthSampHalf < 1
    nbeats = nbeats - 1;
    R_loc(1) = [];
end

if R_loc(end) + beatWidthSampHalf - 1 > length(ecg)
    nbeats = nbeats - 1;
    R_loc(end) = [];
end

if nbeats < 3
    sqiEcg = NaN;
    error_code = 1;
    return
end

% Isolating Beats
beats = zeros(beatWidthSampHalf * 2, nbeats);
for k = 1:nbeats
    sIdx = R_loc(k) - beatWidthSampHalf;
    eIdx = R_loc(k) + beatWidthSampHalf - 1;
    beats(:,k) = ecg(sIdx: eIdx);
end

% Calculating average beat
avgBeat = mean(beats, 2);
pwrAvgBeat = sqrt(mean(avgBeat.^2));

% Calculating SNR for each beat
snr = zeros(nbeats, 1);
for m = 1:nbeats 
    beat = beats(:,m);
    noise = beat - avgBeat;
    pwrNoise = sqrt(mean(noise.^2));
    snr(m) = 20*log10(pwrAvgBeat/pwrNoise);
end

% Calculating SQI of the ECG
sqiEcg = prctile(snr, 25);
error_code = 0;

end

