function [list_clean_r_peaks, no_r_peak_areas_start, no_r_peak_areas_end, removed_beats] = clean_r_peaks(list_r_peaks, nobaseline_ecg, ecg_corrected_snr, snr_stride_length, samplerate)


removed_beats = zeros(0,1);
list_clean_r_peaks = zeros(0,1);
no_r_peak_areas_start = zeros(0,1);
no_r_peak_areas_end = zeros(0,1);

%Remove the last found peak, as it quite offten selekts the T- or P-Wave
removed_beats(end+1) = list_r_peaks(end);
list_r_peaks(end) = [];


min_peak_distance = floor(0.5*samplerate);

peak_list = list_r_peaks;
peak_distance = diff(list_r_peaks);
peak_distance = [peak_distance(1), peak_distance];

short_beats = (peak_distance < min_peak_distance);

if any(short_beats)
    [segments_start, segments_end, segments_value] = enumerate_continuous_segments(short_beats);

    %Find all short segments of "real" beats, and ignore them...
    short_segments = ~segments_value & ((segments_end - segments_start + 1) <= 2);
    if any(short_segments)
        %these segments are max two long...change manually....
        short_beats(segments_start(short_segments)) = true;
        short_beats(segments_end(short_segments)) = true;

        [segments_start, segments_end, segments_value] = enumerate_continuous_segments(short_beats);
        assert(~any(~segments_value & ((segments_end - segments_start + 1) <= 2)));
    end


    %Remove all short_beats
    removed_beats = [removed_beats, peak_list(short_beats)];
    list_clean_r_peaks = peak_list(~short_beats);


    %Area around peak, that should not be considered needing a peak...
    min_peak_radius = floor(0.3*samplerate);

    for ii=1:length(segments_start)
        %Is it a short beat segment?
        if (segments_value)
            if ii==1
                no_r_peak_areas_start(end+1)=0;
            else
                no_r_peak_areas_start(end+1)=peak_list(segments_end(ii-1)) + min_peak_radius;
            end
            if ii==length(segments_start)
                no_r_peak_areas_end(end+1)=length(nobaseline_ecg);
            elsesnr_segments_start
                no_r_peak_areas_end(end+1)=peak_list(segments_start(ii+1)) - min_peak_radius;
            end
        else
            %The last segment might end before the end of the ecg....
            if ii==length(segments_end)
                if (peak_list(segments_end) < (length(nobaseline_ecg) - min_peak_distance - min_peak_radius))
                    no_r_peak_areas_start(end+1)=peak_list(segments_end(ii)) + min_peak_radius;
                    no_r_peak_areas_end(end+1)=length(nobaseline_ecg);
                end
            end
        end
    end

else
%    removed_beats = removed_beats;
    list_clean_r_peaks = peak_list;
end

ecg_corrected_snr_interp1 = interp1(1:snr_stride_length:min((length(ecg_corrected_snr))*snr_stride_length, length(nobaseline_ecg)), ecg_corrected_snr, 1:length(nobaseline_ecg));
low_snr = (ecg_corrected_snr_interp1 < 10);

if any(low_snr)
    no_r_peak_areas = false(size(nobaseline_ecg));
    for ii=1:length(no_r_peak_areas_start)
        no_r_peak_areas(no_r_peak_areas_start(ii):no_r_peak_areas_end(ii)) = true;
    end

     [snr_segments_start, snr_segments_end, low_snr_value] = enumerate_continuous_segments(low_snr);
    snr_segments_start(~low_snr_value) = [];
    snr_segments_end(~low_snr_value) = [];

    for ii=1:length(snr_segments_start)
        if ((snr_segments_end(ii) - snr_segments_start(ii) + 1) > (2 * samplerate))
            no_r_peak_areas(snr_segments_start(ii):snr_segments_end(ii)) = true;
            high_noise_peaks = (list_clean_r_peaks > snr_segments_start(ii)) & (list_clean_r_peaks < snr_segments_end(ii));
	    removed_beats = [removed_beats, list_clean_r_peaks(high_noise_peaks)];
	    list_clean_r_peaks(high_noise_peaks) = [];
        end
    end

    [no_r_peak_areas_start, no_r_peak_areas_end, no_r_peak_areas_value] = enumerate_continuous_segments(no_r_peak_areas);
    no_r_peak_areas_start(~no_r_peak_areas_value) = [];
    no_r_peak_areas_end(~no_r_peak_areas_value) = [];
end






