%% Line: Multiple Observations

% This code will reduce the noise found in observations of y by taking
% averages of k sets of data.

% I use the command 'clear' here to  make sure that there are no unwanted
% values in y, in case we have used the variable name elsewhere before
% this, so that y contains only the values we intend for it to.
clear y

% This command sets up a grid of evenly spaced values that we will use for
% x.  The first value (1) is the starting value for x and the last value
% (20) is the stopping value for x.  The middle value, separated from the
% others by colons, increments x.  Here, this command will result in the
% following values in the vector x: 1, 1.5, 2, 2.5, ..., 19.5, 20
x=1:.5:20;

% k is the number of sets of observations we are making so that we can
% average out the noise.
k=10;

% This command lets us know how many rows and columns are in the variable
% x.  Since x is a vector, m is 1.  The value that n takes depends on the
% increment we used for x.  With the default value, n=39.
[m n]=size(x);

% This set of commands generates the noisy data.  With one for loop inside
% the other, we say that the loops are "nested".  The inner loop generates
% a value of y for a given value of x.  The outer loop generates multiple
% values of y for a given value of x.
for jnd=1:k
    for ind=1:n
        y(ind,jnd)=x(ind)+normrnd(0,.2);
    end
end

% y is a k by p matrix.  Each column of the matrix represents one set of
% observations.  Each row of the matrix represents repeated values of y for
% a fixed value of x.


% This code computes the average at each x value. We use '2' as the second
% entry to the command so that the function will return a n-dimensional
% vector.  If we only typed "mean(y)", the function would return a
% k-dimensional vector, having averaged the wrong quantities.
yavg=mean(y,2);

% One observation of the original noisy data
figure(1)
plot(x,y(:,1))

% A second observation of the noisy data
figure(2)
plot(x,y(:,2))

% The averaged data
figure(3)
plot(x,yavg)


%% Sin Example: Multiple Observations 

% This code will reduce the noise found in observations of y by taking
% averages of k sets of data.  While the first example did this for a
% linear relationship between x and y, we will now consider a more
% complicated function for y.  Each part of this code does the same thing
% as what was done in the above example.  Because of this, I am not going
% to make further comments in this example.

clear y

x=1:.5:20;

k=10;

[m n]=size(x);

for jnd=1:k
    for ind=1:n
        y(ind,jnd)=.1*x(ind)+sin(x(ind))+normrnd(0,.2);
    end
end

yavg=mean(y,2);

figure(4)
plot(x,y(:,1))

figure(5)
plot(x,y(:,2))

figure(6)
plot(x,yavg)

figure(7)
plot(x,.1*x+sin(x))


%% Sin Example: Linear and Median Filters

% This example will continue to work with the definition of y that we used
% previously.  Now, however, we will try to reduce the noise by using the
% linear and median filters.

%This first section of the code generates the data, just like in the
%previous example.  
clear y
x=1:.1:20; % NOTE: This time, we incremented x by only 0.1
[m n]=size(x);
for ind=1:n
    y(ind)=.1*x(ind)+sin(x(ind))+normrnd(0,.05); %NOTE: Less noise then before
end

% This is a plot of the noisy data
figure(8)
plot(x,y)


% This portion of the code applies the linear filter to y

% This sets the first value of ylin to be equal to that of the observation
ylin(1)=y(1);

% This portion of the code takes averages for each y value that isn't
% either the first value or the last
for ind=2:(n-1)
    ylin(ind)=(ylin(ind-1)+y(ind)+y(ind+1))/3;
end

% This averages the last value of y with its preceding value.
ylin(n)=(y(n)+ylin(n-1))/2;


% Here is the data with the linear filter applied to it.
figure(9)
plot(x,ylin)

% This portion of the code applies the median filter to y.

% This sets the first value of ymed to be equal to that of the observation
ymed(1)=y(1);

% This portion of the code computes the median for each y value that isn't
% either the first value or the last
for ind=2:(n-1)
    ymed(ind)=median([ymed(ind-1) y(ind) y(ind+1)]);
end

% This computes the median of the the last value of y with its preceding value.
ymed(n)=median([ymed(n-1) y(n)]);

% Here is the data with the linear filter applied to it.
figure(10)
plot(x,ymed)







