--- a
+++ b/featurebased-approach/subfunctions/residualfeats.m
@@ -0,0 +1,119 @@
+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