Attachment 'Sampling.m'
Download 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
Attached Files
To refer to attachments on a page, use attachment:filename, as shown below in the list of files. Do NOT use the URL of the [get] link, since this is subject to change and can break easily.You are not allowed to attach a file to this page.