%% Create example data
x = 0:0.01:8*pi; % abscissa values
y = x'/5 + sin(10*x)';
plot(x,y);
whos

% Downsample by a factor 10
x2 = x(1:10:end);
y2 = x2'/5 + sin(10*x2)';
figure;
plot(x2,y2);
whos

% Downsample by a factor 100
x3 = x(1:100:end);
y3 = x3'/5 + sin(10*x3)';
figure;
plot(x3,y3);
whos

% Plot for PPT
% figure;
% plot(x,y, 'r');
% hold on;
% ph = plot(x3,y3,'xb');
% set(ph, 'MarkerSize', 10);
% hold on;
% ph = plot(x3,-1,'xk');
% set(ph, 'MarkerSize', 10);
% hold on;
% plot(x3,y3, 'b');

%% Missing data points, interpolation
% create faulty data
y4 = y;
y4(1000:1020) = 0;
figure;
plot(x,y4);

% quick fix: replace by average of edges
y4(1000:1020) = mean( y4([999, 1021]) );
figure;
plot(x,y4);

% linear interpolation
y4(1000:1020) = y4(999) + [1:21]*(y4(1021)-y4(999))/20;
figure;
plot(x,y4);

% create "more faulty" data
y4(1000:1100) = 0;
figure;
plot(x,y4);

% linear interpolation
y4(1000:1100) = y4(999) + [1:101]*(y4(1101)-y4(999))/100;
figure;
plot(x,y4);

%% Signal-to-noise (SNR)
clear all;  % let's start afresh
close all;

% Make a signal
x = 0:0.01:4*pi; % abscissa values
y1 = sin( x );
y1(1:end/2) = 0; % include "baseline"
plot(x,y1);

% Make a slightly different signal
y2 = 1.1 * y1;
hold on;
plot(x,y2, 'r')
title('y1 and y2');

% Add some noise
lx = length(x);
y1 = y1 + 0.05*randn(1,lx);
y2 = y2 + 0.05*randn(1,lx);
figure;
plot(x,y1);
hold on;
plot(x,y2, 'r');
title('y1 and y2 with noise');

% compute SNR (power)
vary1 = var(y1(1:end/2));
vary2 = var(y2(1:end/2));
figure;
plot(x,y1.^2/vary1);
hold on;
plot(x,y2.^2/vary2, 'r'); % "zoom" into baseline to check it fluctuates around 1
title('SNR, power');
% Note: the power "rectifies" the signal, i.e. it's all positive

% compute SNR (standard deviation)
stdy1 = std(y1(1:end/2));
stdy2 = std(y2(1:end/2));
figure;
plot(x,y1/stdy1);
hold on;
plot(x,y2/stdy2, 'r'); % "zoom" into baseline to check it fluctuates around 1
title('SNR, std');

% Root-mean-square (RMS)
rmsy1 = sqrt( mean ( y1(1:end/2).^2 ) );
rmsy2 = sqrt( mean ( y2(1:end/2).^2 ) );
figure;
plot(x,y1/rmsy1);
hold on;
plot(x,y2/rmsy2, 'r'); % "zoom" into baseline to check it fluctuates around 1
title('SNR, RMS');


%% Error propagation
% subtraction
diffy = y1-y2;
figure;
plot(x, diffy);
title('Diff');

var( y1(1:600) )
var( y2(1:600) )
var( diffy(1:600) )

addy = y1+y2;
var( addy(1:600) )

% Error propagation can be disastrous
figure;
y1ori = sin( x );
plot(x, y1ori); % plot sine curve

hold on;
ydiff = diff(y1ori);
plot(x(2:end), 100*ydiff, 'r'); % plot "clean" derivate of sine curve
title('Diff2')

plot(x(2:end), 100*diff(y1), 'g'); % plot noisy derivative of sine curve

%% Smoothing (get previous figure 'Diff' back)
% command "smooth" is easy but not necessarily the best way for data filtering
figure;
sy = smooth(diffy, 10);
plot(x, sy, 'r');
title('smooth')

sy = smooth(diffy, 100);    % smooth more
hold on;
plot(x, sy, 'k');

% What could possibly go wrong?

% Edge effects of filtering:
figure;
plot(x(2:end), diff(sy));

% Oversmoothing:
sy = smooth(diffy, 1000);   % smooth too much
plot(x, sy, 'g');






%% The END