function [ ecg ] = preprocess_wavelet(ecg)
%preprocess_wavelet   Preprocess the ECG using wavelets

%Flip ECG if necessary
flipBack = 0;
if size(ecg,1) > size(ecg,2)
    ecg = ecg';
    flipBack = 1;
end

for i = 1:min(size(ecg))
    ecg(i,:) = ecg(i,:) - bwqvr(10000, ecg(i,:)); %Remove baseline
    [C, L] = wavedec(ecg(i,:), 4, 'db4');
    
    a1 = appcoef(C, L, 'db4', 1);
    d1 = detcoef(C, L, 1);
    a2 = appcoef(C, L, 'db4', 2);
    d2 = detcoef(C, L, 2);
    a3 = appcoef(C, L, 'db4', 3);
    d3 = detcoef(C, L, 3);
    a4 = appcoef(C, L, 'db4', 4);
    d4 = detcoef(C, L, 4);
    
    % set detail coefficients below threshold equal to zero
    
    sigma1 = median(abs(d1))/0.6745;
    threshd1 = sigma1*sqrt(2*log(length(ecg(i,:))));
    % thresholding on d1
    for j=1:length(d1)
        
        if abs(d1(j))>threshd1
            d1_kept(j) = d1(j);
        elseif abs(d1(j))<threshd1
            d1_kept(j) = 0;
        end
        
    end
    
    sigma2 = median(abs(d2))/0.6745;
    threshd2 = sigma2*sqrt(2*log(length(ecg(i,:))));
    
    % thresholding on d2
    
    for j=1:length(d2)
        if abs(d2(j))>threshd2
            d2_kept(j) = d2(j);
            
        elseif abs(d2(j))<threshd2
            d2_kept(j) = 0;
            
        end
        
    end
    coeff = [a4, d4, d3, d2_kept, d1_kept];
    % reconstruction
    ecg(i,:) = waverec(coeff, L, 'db4');
end

if flipBack == 1
    ecg = ecg';
end
end

