function noise_estimate = ecg_estimate_noise( ecg_vector, r_peak_list, samplerate )

%% Part 0; Initializations and Error Avoidance


    MATRIX_ROW_COUNT = 100;
    SECONDS_BEFORE_R = 0.5;
    SECONDS_AFTER_R = 0.75;

%    if (nargout > 1)
%       t_onset_list, t_length_list
%        error('cvp:not_implemented', 'signal estimate not implemented yet');
%    end

    if ~isvector(ecg_vector) || ~isnumeric(ecg_vector)
        error('cvp:ecg_estimate_noise:ecg_not_vector', 'ecg_vector is not a vector');
    end

    if ~isvector(r_peak_list) || ~isnumeric(r_peak_list)
        error('cvp:ecg_estimate_noise:r_peak_list_not_vector', 'r_peak_list is not a numeric vector');
    end

    if (min(r_peak_list) < 1 || max(r_peak_list) > length(ecg_vector))
        error('cvp:ecg_estimate_noise:invalid_r_peak_list', 'r_peak_list contains invalid values');
    end

    if (length(r_peak_list) < 2)
        error('cvp:ecg_estimate_noise:not_enough_r_peaks', 'Not enough R-Peaks Provided');
    end

    if ~isscalar(samplerate) || ~isnumeric(samplerate) || (fix(samplerate) ~= samplerate)
        error('cvp:ecg_estimate_noise:samplerate_not_integer', 'samplerate must be an integer');
    end


    noise_estimate = nan(1,length(ecg_vector));

    if size(ecg_vector, 1) ~= 1
        ecg_vector = reshape(ecg_vector, 1, length(ecg_vector));
    end

%% Part 1: Estimate the Noise around all R-Peaks
    segment_start_offset = floor(SECONDS_BEFORE_R * samplerate);
    segment_end_offset = floor(SECONDS_AFTER_R * samplerate);

    verified_r_peak_list = r_peak_list(r_peak_list > segment_start_offset & r_peak_list < (length(ecg_vector) - segment_end_offset));
    clear('r_peak_list');

    segment_start = reshape(verified_r_peak_list - segment_start_offset, [], 1); %Make sure it's a Mx1 Vector

    segment_end = reshape(verified_r_peak_list + segment_end_offset, [], 1); %Make sure it's a Mx1 Vector
    segment_length = segment_end(1) - segment_start(1) + 1;

    %The Matrix thrice as big as the number of rows we consider each
    %iteration (To avoid sudden changes in the calculated noise) and
    %in every iteration only the middle third is used (except in the
    %first and last one)
    window_size = floor(MATRIX_ROW_COUNT / 3);
    window_count = ceil(length(verified_r_peak_list) / window_size) - 2; %The first and last windows are special cases
    window_start = window_size + 1;
    window_end = 2*window_size;

    for j=2:(window_count+1)
        r_peak_offset = (1:MATRIX_ROW_COUNT) + (window_size * (j-2));
        if (j == (window_count+1))
            r_peak_offset(r_peak_offset > length(segment_start)) = r_peak_offset(r_peak_offset > length(segment_start)) -  window_size;
        end
%        assert(all(r_peak_offset <= length(segment_start))); %Just to avoid nasty suprises

        %Create a vector with the relevant index at the correct position
        ecg_mat_ind = segment_start(r_peak_offset) * ones(1,segment_length) + ones(1,MATRIX_ROW_COUNT)' * (0:(segment_length-1));
        %Retrieve the Values of the Matrix
        ecg_mat = ecg_vector(ecg_mat_ind);
if true
        [U,S,V] = svd(ecg_mat,'econ');

        S_noise_diag = diag(S);
        S_noise_diag(1:5) = 0;
        S_noise = diag(S_noise_diag);
else
        [U,S,V] = svd(ecg_mat);

        S_noise_diag = diag(S);
        S_noise_diag(1:5) = 0;
        S_noise = horzcat(vertcat(diag(S_noise_diag), zeros(length(U) - length(S_noise_diag), length(S_noise_diag))), zeros(length(S_noise_diag), length(V) - length(S_noise_diag)));
end
        ecg_mat_noise = U * S_noise * V';

        if (j==2)
            for ii = 1:(window_start-1)
                [noise_estimate_update, update_start, update_end] = add_to_estimate(r_peak_offset(ii), squeeze(ecg_mat_noise(ii, :)), noise_estimate, segment_start, segment_end, verified_r_peak_list, segment_start_offset);
                noise_estimate(update_start:update_end) = noise_estimate_update;
            end
        end

        for ii= window_start:window_end;
            [noise_estimate_update, update_start, update_end] = add_to_estimate(r_peak_offset(ii), squeeze(ecg_mat_noise(ii, :)), noise_estimate, segment_start, segment_end, verified_r_peak_list, segment_start_offset);
            noise_estimate(update_start:update_end) = noise_estimate_update;
        end

        if (j == (window_count+1))
            for ii= (window_end+1):MATRIX_ROW_COUNT
                [noise_estimate_update, update_start, update_end] = add_to_estimate(r_peak_offset(ii), squeeze(ecg_mat_noise(ii, :)), noise_estimate, segment_start, segment_end, verified_r_peak_list, segment_start_offset);
                noise_estimate(update_start:update_end) = noise_estimate_update;
            end
        end
    end


%% Part 2: Where there is no R-Peak around, the Signal is Noise

    noise_estimate(1:2) = nan;       %Make sure there is always some noise
    noise_estimate(end-1:end) = nan; %at the beginning and ending....
    [segments_start, segments_end, undefined_noise_segments] ...
        = enumerate_continuous_segments(isnan(noise_estimate));

    for ii = 1:length(undefined_noise_segments)
       if (ii == 1) || (undefined_noise_segments(ii))
           %No R-Peak: All Signal in the ECG Vector is noise
           ecg_noise = ecg_vector(segments_start(ii):segments_end(ii));

           noise_estimate(segments_start(ii):segments_end(ii)) = ecg_noise;
       else
           %This R-Peak Segment is adjacent to undefined segments...fade
           %over to avoid sudden changes in noise level
%           assert((ii > 1) && (undefined_noise_segments(ii-1) == true));
%           assert((ii < length(undefined_noise_segments)) && (undefined_noise_segments(ii+1) == true));

           %Fade time is half the distance to the R-Peak
           front_fade_length = floor(segment_start_offset / 2);
           tail_fade_length = floor(segment_end_offset / 2);

           front_fade_start = segments_start(ii);
           front_fade_end = segments_start(ii) + front_fade_length;

           tail_fade_end = segments_end(ii);
           tail_fade_start =  segments_end(ii) - tail_fade_length;

           assert(front_fade_end < tail_fade_start);

           front_fade_noise = noise_estimate(front_fade_start:front_fade_end);
           front_fade_ecg = ecg_vector(front_fade_start:front_fade_end);

           noise_estimate(front_fade_start:front_fade_end) = fade_signals(front_fade_ecg, front_fade_noise);

           tail_fade_noise = noise_estimate(tail_fade_start:tail_fade_end);
           tail_fade_ecg = ecg_vector(tail_fade_start:tail_fade_end);

           noise_estimate(tail_fade_start:tail_fade_end) = fade_signals(tail_fade_noise, tail_fade_ecg);
       end

    end





%% Part 3: The Result
    %Just make sure that we did not introduced any new NaNs (even though
    %the ecg_vector should not contain nans to begin with...
    noise_nans = isnan(noise_estimate);
    ecg_nans = isnan(ecg_vector);
%    assert(all(~noise_nans | ecg_nans)); % noise_nan => ecg_nan, rewrite of boolean implication

end
% Author: Jan F. Kraemer <jan.kraemer@physik.hu-berlin.de>
%

%
% Helper Functions for Part 1
%


function mixing_vec = eval_mixing_function(start_intersect_idx, stop_intersect_idx, time_eq)
    idx_offset = start_intersect_idx - 1;
    idx_length = stop_intersect_idx - start_intersect_idx + 3;
    stretch_modifier = - log(time_eq) / log(2);
    mixing_function = @(x) 0.5 - cos(((x - idx_offset) ./ idx_length) .^ stretch_modifier .* pi) / 2;
    mixing_vec = mixing_function(start_intersect_idx:stop_intersect_idx);
end


function [noise_estimate_update, update_start, update_end] = add_to_estimate(r_peak_idx, estimation_vec, noise_estimate, segment_start, segment_end, verified_r_peak_list, segment_start_offset)
    update_start = segment_start(r_peak_idx);
    update_end = segment_end(r_peak_idx);
    noise_estimate_update = noise_estimate(update_start:update_end);
    start_idx = segment_start(r_peak_idx);
    estimation_offset = (verified_r_peak_list(r_peak_idx) - segment_start_offset);

    if (r_peak_idx > 1) && (segment_end(r_peak_idx-1) > segment_start(r_peak_idx))
        start_intersect_idx = max(verified_r_peak_list(r_peak_idx-1) + 1, segment_start(r_peak_idx));
        stop_intersect_idx = min(segment_end(r_peak_idx-1), verified_r_peak_list(r_peak_idx) - 1);

%        assert(start_intersect_idx < stop_intersect_idx);

        %Calculate Point of equal uncertainty
        %(The certainty decreases linearly from the r-peak to the
        % end/start of the segment)
        % x = -(b*c - a*d) / (a-b-c+d)
        % for
        a = verified_r_peak_list(r_peak_idx-1); b= segment_end(r_peak_idx-1);
        c = verified_r_peak_list(r_peak_idx); d = segment_start(r_peak_idx);
        % and iff a + d ~= b + c
%        assert(a + d ~= b + c);
        %assert(a~=b);
        %assert(c~=d);

        time_eq = -(b*c - a*d) / (a-b-c+d);

        % The Mixing function is a cosine [0,pi] mapped to
        % [start_intersect_idx, stop_intersect_idx] adjusted to hit 0.5
        % at time_eq, see eval_mixing_function

        mixing_vec = eval_mixing_function(start_intersect_idx, stop_intersect_idx, time_eq);

        noise_estimate_update((start_intersect_idx:stop_intersect_idx)-update_start+1) = ...
            (noise_estimate(start_intersect_idx:stop_intersect_idx) .* (1-mixing_vec)) ...
            + ...
            (estimation_vec((start_intersect_idx:stop_intersect_idx) - estimation_offset + 1) .* mixing_vec);
        start_idx = stop_intersect_idx + 1;
    end
    noise_estimate_update((start_idx:segment_end(r_peak_idx)) - update_start + 1) = estimation_vec((start_idx-estimation_offset + 1):end);
%    assert(any(~isnan(noise_estimate_update)));
end

%
% Helper Functions for Part 2
%

function [mixed_signal] = fade_signals(start_signal, end_signal)
    %assert(length(start_signal) == length(end_signal))
    signal_length = length(start_signal);
    mix_length = length(start_signal) + 1;

    fade_modifier = (1:signal_length) ./ mix_length;

    mixed_signal = (start_signal .* (1 - fade_modifier)) + (end_signal .* fade_modifier);
end
