% simulate1fnoise(number of samples,
%                 signal amplitude ,
%                 signal frequency ,
%                 number of redundant sensors to simulate (each has independent noise),
%                 starting probability of each update)
% Outputs:
%
% pure_signal = physical signal we want to measure
% ideal_measurement = the best we could expect, after denoising (additive
%                                                       noise limited)
% additive_noise    = the additive noise
% acquisition_signal = Simulated measurement. What we would actually
% obtain.
% update_noise       = the step-like noise which we want to remove, due to
%                      the updates
% refUpdates         = time instants at which the reference is updated, and we
%                      have step-function noise

% The information we have available in the real case is refUpdates and
% acquisition_signal.


nsamp = 1000; t = 1:nsamp; nSensors = 10;
[pure_signal,ideal_measurement,additive_noise,acquisition_signal,update_noise,refUpdates] = simulate1fnoise(nsamp,4,0.01,nSensors,0.8);


figure, subplot(1,2,1), plot(t,ideal_measurement(1,:),t,pure_signal(1,:),'r'),
        title('Ideal case for measurement (no updates)'), legend({'Ideal Measurement', 'Pure Signal'})
        xlabel('time')
        subplot(1,2,2), plot(t,acquisition_signal(1,:),t,update_noise(1,:),'r'),
        title('Actual measurement'), legend({'Acquisition', 'Step noise'})
        xlabel('time')


function [pure_signal,ideal_measurement,additive_noise,acquisition_signal,update_noise,refUpdates] = simulate1fnoise(nsamp,amp,freq,nWindows,probabilityUpdate)


t = 1:nsamp;
pure_signal = amp*sin(2*pi*freq*t).*[zeros(1,floor(nsamp/4)) hamming(floor(nsamp/2))' zeros(1,nsamp-(floor(nsamp/2) + floor(nsamp/4)))] ;

additive_noise = randn(nWindows,nsamp);
ideal_measurement = pure_signal + additive_noise; % What one would ideally measure
update_noise = zeros(nWindows,nsamp);   % Noise resulting from updates (the part that we want to discover/remove from our signal)
acquisition_signal = zeros(nWindows,nsamp);    % Actual physical measurement

refUpdates = zeros(nWindows,nsamp);

for wnd = 1:nWindows
    updateError = 0;
    countNoUpdates = -1;
    for idx = 1:nsamp
        countNoUpdates = countNoUpdates+1;
        acquisition_signal(wnd,idx) = ideal_measurement(wnd,idx) + updateError;
        update_noise(wnd,idx) = updateError;
        if rand > probabilityUpdate-(countNoUpdates*0.025)
            refUpdates(wnd,idx) = 1;
            updateError = updateError+additive_noise(wnd,idx);
            countNoUpdates = 0;
        end
    end
end

end