1 %% Create example data
   2 x = 0:0.01:8*pi; % abscissa values
   3 y = x'/5 + sin(10*x)';
   4 plot(x,y);
   5 whos
   6 
   7 % Downsample by a factor 10
   8 x2 = x(1:10:end);
   9 y2 = x2'/5 + sin(10*x2)';
  10 figure;
  11 plot(x2,y2);
  12 whos
  13 
  14 % Downsample by a factor 100
  15 x3 = x(1:100:end);
  16 y3 = x3'/5 + sin(10*x3)';
  17 figure;
  18 plot(x3,y3);
  19 whos
  20 
  21 % Plot for PPT
  22 % figure;
  23 % plot(x,y, 'r');
  24 % hold on;
  25 % ph = plot(x3,y3,'xb');
  26 % set(ph, 'MarkerSize', 10);
  27 % hold on;
  28 % ph = plot(x3,-1,'xk');
  29 % set(ph, 'MarkerSize', 10);
  30 % hold on;
  31 % plot(x3,y3, 'b');
  32 
  33 %% Missing data points, interpolation
  34 % create faulty data
  35 y4 = y;
  36 y4(1000:1020) = 0;
  37 figure;
  38 plot(x,y4);
  39 
  40 % quick fix: replace by average of edges
  41 y4(1000:1020) = mean( y4([999, 1021]) );
  42 figure;
  43 plot(x,y4);
  44 
  45 % linear interpolation
  46 y4(1000:1020) = y4(999) + [1:21]*(y4(1021)-y4(999))/20;
  47 figure;
  48 plot(x,y4);
  49 
  50 % create "more faulty" data
  51 y4(1000:1100) = 0;
  52 figure;
  53 plot(x,y4);
  54 
  55 % linear interpolation
  56 y4(1000:1100) = y4(999) + [1:101]*(y4(1101)-y4(999))/100;
  57 figure;
  58 plot(x,y4);
  59 
  60 %% Signal-to-noise (SNR)
  61 clear all;  % let's start afresh
  62 close all;
  63 
  64 % Make a signal
  65 x = 0:0.01:4*pi; % abscissa values
  66 y1 = sin( x );
  67 y1(1:end/2) = 0; % include "baseline"
  68 plot(x,y1);
  69 
  70 % Make a slightly different signal
  71 y2 = 1.1 * y1;
  72 hold on;
  73 plot(x,y2, 'r')
  74 title('y1 and y2');
  75 
  76 % Add some noise
  77 lx = length(x);
  78 y1 = y1 + 0.05*randn(1,lx);
  79 y2 = y2 + 0.05*randn(1,lx);
  80 figure;
  81 plot(x,y1);
  82 hold on;
  83 plot(x,y2, 'r');
  84 title('y1 and y2 with noise');
  85 
  86 % compute SNR (power)
  87 vary1 = var(y1(1:end/2));
  88 vary2 = var(y2(1:end/2));
  89 figure;
  90 plot(x,y1.^2/vary1);
  91 hold on;
  92 plot(x,y2.^2/vary2, 'r'); % "zoom" into baseline to check it fluctuates around 1
  93 title('SNR, power');
  94 % Note: the power "rectifies" the signal, i.e. it's all positive
  95 
  96 % compute SNR (standard deviation)
  97 stdy1 = std(y1(1:end/2));
  98 stdy2 = std(y2(1:end/2));
  99 figure;
 100 plot(x,y1/stdy1);
 101 hold on;
 102 plot(x,y2/stdy2, 'r'); % "zoom" into baseline to check it fluctuates around 1
 103 title('SNR, std');
 104 
 105 % Root-mean-square (RMS)
 106 rmsy1 = sqrt( mean ( y1(1:end/2).^2 ) );
 107 rmsy2 = sqrt( mean ( y2(1:end/2).^2 ) );
 108 figure;
 109 plot(x,y1/rmsy1);
 110 hold on;
 111 plot(x,y2/rmsy2, 'r'); % "zoom" into baseline to check it fluctuates around 1
 112 title('SNR, RMS');
 113 
 114 
 115 %% Error propagation
 116 % subtraction
 117 diffy = y1-y2;
 118 figure;
 119 plot(x, diffy);
 120 title('Diff');
 121 
 122 var( y1(1:600) )
 123 var( y2(1:600) )
 124 var( diffy(1:600) )
 125 
 126 addy = y1+y2;
 127 var( addy(1:600) )
 128 
 129 % Error propagation can be disastrous
 130 figure;
 131 y1ori = sin( x );
 132 plot(x, y1ori); % plot sine curve
 133 
 134 hold on;
 135 ydiff = diff(y1ori);
 136 plot(x(2:end), 100*ydiff, 'r'); % plot "clean" derivate of sine curve
 137 title('Diff2')
 138 
 139 plot(x(2:end), 100*diff(y1), 'g'); % plot noisy derivative of sine curve
 140 
 141 %% Smoothing (get previous figure 'Diff' back)
 142 % command "smooth" is easy but not necessarily the best way for data filtering
 143 figure;
 144 sy = smooth(diffy, 10);
 145 plot(x, sy, 'r');
 146 title('smooth')
 147 
 148 sy = smooth(diffy, 100);    % smooth more
 149 hold on;
 150 plot(x, sy, 'k');
 151 
 152 % What could possibly go wrong?
 153 
 154 % Edge effects of filtering:
 155 figure;
 156 plot(x(2:end), diff(sy));
 157 
 158 % Oversmoothing:
 159 sy = smooth(diffy, 1000);   % smooth too much
 160 plot(x, sy, 'g');
 161 
 162 
 163 
 164 
 165 
 166 
 167 %% The END