function feats = residualfeats(signal,fs,qrs)
% This function derives features out of QRS cancelled ECG signals, aiming at the atrial information.
%
% --
% ECG classification from single-lead segments using Deep Convolutional Neural 
% Networks and Feature-Based Approaches - December 2017
% 
% Released under the GNU General Public License
%
% Copyright (C) 2017  Fernando Andreotti, Oliver Carr
% University of Oxford, Insitute of Biomedical Engineering, CIBIM Lab - Oxford 2017
% fernando.andreotti@eng.ox.ac.uk
%
% 
% For more information visit: https://github.com/fernandoandreotti/cinc-challenge2017
% 
% Referencing this work
%
% Andreotti, F., Carr, O., Pimentel, M.A.F., Mahdi, A., & De Vos, M. (2017). 
% Comparing Feature Based Classifiers and Convolutional Neural Networks to Detect 
% Arrhythmia from Short Segments of ECG. In Computing in Cardiology. Rennes (France).
%
% Last updated : December 2017
% 
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
% 
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
% 
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

feats = zeros(1,5);
signal = signal';

if length(qrs) > 5
    residual = FECGSYN_ts_extraction(qrs,signal,'TS-CERUTTI',0,5,2,fs);
else
    return
end

try
    % Cross-correlation
    R=corrcoef(signal,residual);
    feats(1) = 1 - R(1,2);
    
    % Spectral coherence
    [Pmf,f]=mscohere(signal,residual,[],[],1024,fs);
    feats(2)=1-mean(Pmf(f>3&f<30));
    
    feats(3) = psqi(residual',fs);
    
    % Fundamental frequency on spectrum against rest
    feats(4)=mpsqi1(residual',fs);
    
    % Test for randomness
    [h,p,R] = white_test(residual', 1, rand(1e3,1));
    feats(5) = 1-p;
catch
    disp('Residuals failed!')
end

end

function [h, p_value, R] = white_test(x, maxlag, alpha)
% tests against the null hyptothesis of whiteness
% see http://dsp.stackexchange.com/questions/7678/determining-the-whiteness-of-noise for details
% demo:
% % white:
% >> [h,p,R]=white_test(((filter([1], [1], rand(1e3,1)))))
% >> h = 0, p =1, R = 455
% % non white:
% >> [h,p,R]=white_test(((filter([1 .3], [.4 0.3], rand(1e3,1)))))
% >> h = 1, p = 0, R = 2e3
%
% Copyright (c) 2013, Hanan Shteingart
% All rights reserved.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
%
% * Redistributions of source code must retain the above copyright
% notice, this list of conditions and the following disclaimer.
% * Redistributions in binary form must reproduce the above copyright
% notice, this list of conditions and the following disclaimer in
% the documentation and/or other materials provided with the distribution
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.

N = length(x);
x = x-mean(x);
if ~exist('m','var')
    maxlag = N-1;
end
if ~exist('alpha','var')
    alpha = 0.05;
end
[r, lag] = xcorr(x, maxlag, 'biased');
R = N/r(lag==0)^2*sum(r(lag > 0).^2);
p_value = 1-chi2cdf(R, maxlag);
T = chi2inv(1-alpha, maxlag);
h =  R > T;
end