function [signal_power, noise_power, snr, noise_color, color_corrected_snr] = calculate_signal_properties (signal_vector, noise_vector, samplerate, stride_length, window_radius)
%CALCULATE_SIGNAL_PROPERTIES Calculates Signal and Noise power,
%Signal-to-Noise Ratio and the Color of the noise at each step in a given
%radius

%% Part 1: Check Inputs
    if ~isvector(signal_vector) || ~isnumeric(signal_vector)
        error('cvp:calculate_signal_properties:signal_not_vector', 'signal_vector is not a numeric vector');
    end

    if ~isvector(noise_vector) || ~isnumeric(noise_vector)
        error('cvp:calculate_signal_properties:noise_not_vector', 'noise_vector is not a numeric vector');
    end

    if size(signal_vector) ~=  size(noise_vector)
        error('cvp:calculate_signal_properties:signal_noise_not_equal_size', 'signal_vector and noise_vector must be of equal size');
    end

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

    if ~isscalar(stride_length) || ~isnumeric(stride_length) || (fix(stride_length) ~= stride_length)
        error('cvp:calculate_signal_properties:stride_length_not_integer', 'stride_length must be an integer');
    end

    if ~isscalar(window_radius) || ~isnumeric(window_radius) || (fix(window_radius) ~= window_radius)
        error('cvp:calculate_signal_properties:window_radius_not_integer', 'window_radius must be an integer');
    end

%% Part 2: Initializations

    signal_length = length(signal_vector);

    step_count = floor(signal_length / stride_length);
    step_start_delta = floor( (stride_length-1) / 2);
    step_end_delta = (stride_length - 1) - step_start_delta;

    signal_power = nan(step_count, 1);
    noise_power = nan(step_count, 1);
    snr = nan(step_count, 1);
    noise_color = nan(step_count, 1);
    color_corrected_snr = nan(step_count, 1);

    edge_length = window_radius + 1;

    first_step = 1 + edge_length;
    last_step = step_count - edge_length;

    window_size = (1 + 2 * window_radius) * stride_length;
    fft_calc_size = 2 ^ nextpow2(window_size);

    %Helper Vectors
    fft_window = hann(window_size,  'periodic')';
    freqency_helper = log(((1:(fft_calc_size / 2)) ./ (fft_calc_size)) * samplerate);
    fft_helper =  1:(fft_calc_size / 2);





%% Part 3:    Calculate the properties


    for ii=first_step:last_step
        step_idx = stride_length * (ii-1) - step_start_delta;
        segment_start_idx = stride_length * (ii-1) - step_start_delta - (window_radius * stride_length);
        segment_end_idx   = stride_length * (ii-1) + step_end_delta   + (window_radius * stride_length);

if false
        signal_segment = detrend(signal_vector(segment_start_idx:segment_end_idx));
        noise_segment = detrend(noise_vector(segment_start_idx:segment_end_idx));
else
        signal_segment = signal_vector(segment_start_idx:segment_end_idx) - mean(signal_vector(segment_start_idx:segment_end_idx));
        noise_segment = noise_vector(segment_start_idx:segment_end_idx) - mean(noise_vector(segment_start_idx:segment_end_idx));
end

        noise_fft = fft(noise_segment .* fft_window, fft_calc_size);
        noise_frequency_power = noise_fft(fft_helper) .* conj(noise_fft(fft_helper)) ./ (fft_calc_size / 2);


        noise_coeff = polyfit(freqency_helper(2:end), log(noise_frequency_power(2:end)), 1);

        noise_beta = -1 * noise_coeff(1);


        noise_power_val = sum(abs(noise_segment) .^ 2) / (window_size / samplerate);
        signal_power_val = sum(abs(signal_segment) .^ 2) / (window_size / samplerate);

        noise_color(ii) = noise_beta;

        %TODO: Vectorize ... should be straight forward.
        signal_power(ii) = signal_power_val;
        noise_power(ii) = noise_power_val;

        if (signal_power_val < noise_power_val) || isnan(signal_power_val) || (signal_power_val < eps)
            snr_val = 1;
        elseif (signal_power_val > 1000 * noise_power_val) || isnan(noise_power_val)
            snr_val = 1000;
        else
            snr_val = max(1, min(signal_power_val / noise_power_val, 1000));
        end
        snr_val = max(0, min(10*log10(snr_val), 30));
        snr(ii) = snr_val;


        %Noise is not equally "disturbing" for any noise color. If the
        %noise is dominant in the low frequencies (beta = 2) we do not care
        %much, since the noise disturbes mostly the large features.
        %If it is very strong in the higher frequencies (beta = 0) all
        %interesting features will be hard to discern.
        %Thus we want to correct for it.
        %This is done rather empirical in the moment.
        %
        %The snr is taken as "True Value(TM)" for beta = 1. Any higher and we
        %want to reduce the impact of noise, any lower and we want to increase it.
        %We define upper and lower multiplier plateu values
        multiplier_plateu_high = 5;
        multiplier_plateu_low  = 0.2;
        %That are valid from the following noise colors
        color_cutoff_high = 1.75;
        color_cutoff_low = 0.25;

        %In between the values lie on a sigmoid like curve (a [0, pi] cos
        %in this case), with is scaled to hit both plateus and (1, 1)

        bounded_beta = max(0, min(noise_beta, 2));
        %correction_factor  = nan;
        if (bounded_beta > color_cutoff_high)
            correction_factor = multiplier_plateu_high;
        elseif (bounded_beta < color_cutoff_low)
            correction_factor = multiplier_plateu_low;
        else
            correction_factor = ...
                (multiplier_plateu_high - multiplier_plateu_low) * ...
                (1/2 - 1/2 * cos( ...
                    pi * (...
                            (bounded_beta - color_cutoff_low) / ...
                            (color_cutoff_high - color_cutoff_low)...
                         ) ^ ( ...
                                 (-1 * log( (1 - color_cutoff_low) / (color_cutoff_high - color_cutoff_low))) / ...
                                 log(2) ...
                             ) ...
                    )) + ...
                multiplier_plateu_low;
        end

        color_corrected_snr(ii) = snr_val * correction_factor;
    end
%% Part 4: The result
% Well...nothing else to do...

end
% Author: Jan F. Kraemer <jan.kraemer@rwth-aachen.de>
%

