% ECE 248, Kalman and Adaptive Filtering, Fall 2010, UCSB
% written by: Katie Byl, Sep 26, 2010
% edit list: 
%    Oct. 10, 2010 - minor changes, to post on the web
%
% Description:
%
% Toward justifying Gaussian noise is a "good" description in
% general, we use a toy example...  A ball is dropped and goes
% through a series of binary processes (random walk).

DOPLOT = true;   % Set to "false" to avoid animation sequence...
%DOPLOT = false;
Nlevel = 9;  % number of binary decisions (+1 vs -1)
Nbin = Nlevel+1;  % number of bins at bottom

hF = figure(10); % handle to figure
set(hF,'Position',[10 20 450 600])  % make window "taller"
clf; hold on
set(hF,'NumberTitle','off')
set(hF,'Name','A Random Walk vs Gaussian Noise')
set(gcf,'Color',[1 1 1]) % set from gray to white, as background
axis off % get rid of tick marks and axes, etc, to look "cleaner"

baseval = -15;  % delta y value for plotting
R = .1; % radius of the ball
triR = .3;
dy = 1; % delta y for one level...
dx = 1;
xfinal = dx*(-(Nbin-1)/2:(Nbin-1)/2);
avals= 2*pi*[0:1/30:1];
xb = R*cos(avals); % a ball; x values
yb = R*sin(avals); % ball y values (to plot)
if DOPLOT
for n=1:length(xfinal)
    xoff = xfinal(n);
    plot(xoff+[-1 -1 1 1]*R,[.6 baseval baseval .6],'k-','LineWidth',2,...
        'Color',[1 1 1]*.6)
end

for nlev = 1:Nlevel
    yoff = (Nlevel+1-nlev)*dy;
    for nbin = 1:nlev
        xoff = -(nlev-1)/2 + (nbin-1)*dx;
        patch(xoff+[0 -1 1 0]*triR,yoff+[0 -1 -1 0]*.25,'b-','LineWidth',2,...
            'FaceColor',[0 .7 0])
        plot(xoff-R+[-.05 0],yoff+[.7 .3],'k-','LineWidth',2,...
            'Color',[1 1 1]*.6)
        plot(xoff+R-[-.05 0],yoff+[.7 .3],'k-','LineWidth',2,...
            'Color',[1 1 1]*.6)
    end
end
end

bin_ct = 0*[1:Nbin];
for Ndo = 1:200  % Drop 200 balls
    whatway = (-1).^[round(rand(1,Nlevel))]; %  "left" or "right" decisions (50/50)
    dxgo = .5*dx*cumsum(whatway); % x pts on ball path
    dygo = -dy*[1:length(whatway)]; % y pts on ball path
    dygo = dygo-dygo(end)+R+.05;
    xtot = zeros(1,2*length(dxgo));
    ytot = xtot;
    xtot(1:2:end) = dxgo;  % Actual pts to draw a "pretty" path
    xtot(2:2:end) = dxgo;
    ytot(1:2:end) = dygo+.65; % pts in y
    ytot(2:2:end) = dygo;
    
    xtot(end+1) = xtot(end);
    ytot(end+1) = baseval+bin_ct(1+dxgo(end)-xfinal(1))*2*R + R;
    xtot = [0 0 xtot];
    ytot = [ytot(1)+[1.5 .4], ytot];
    
    if DOPLOT  % Animate each ball as it falls and collects in a bin at bottom
        %plot(xtot,ytot,'r.-','LineWidth',2)
        p1 = patch(xtot(1)+xb,ytot(1)+yb,'r-');
        axis image
        doway = .25;  % The smaller this is, the slower the animation.
        if Ndo>40  % After a while "fast-forward" the frames
            doway = 3;
        elseif Ndo>=6 % (intermediate frame speed value)
            doway = 1;
        end
        for nway = [1:doway:length(xtot)-eps, length(xtot)]
            set(p1,'XData',interp1(1:length(xtot),xtot,nway)+xb,...
                'YData',interp1(1:length(ytot),ytot,nway)+yb);
            drawnow;
            if Ndo<=3
                pause(.02);
            end
        end
    end
    
    % Now, add a new BALL to this BIN:
    bin_ct(1+dxgo(end)-xfinal(1)) = bin_ct(1+dxgo(end)-xfinal(1)) + 1;
    
end

DoTot = 1e3;
DoMore = DoTot - sum(bin_ct);
QuickDo = (-1).^round(rand(DoMore,Nlevel));
QuickEnd = .5*sum(QuickDo,2);

% Now, show the results, compared with expected values:
hF2 = figure(2); clf
clf; hold on
set(hF2,'NumberTitle','off')
set(hF2,'Name','After Fast-forwarding on the random walk...')
pbase = (1/2)^Nlevel;
for n=1:length(bin_ct)
    bin_ct(n) = bin_ct(n) + sum(QuickEnd==xfinal(n));
    if n>1
        pbase(n) = pbase(n-1)*(Nlevel+2-n)/(n-1);
    end
    patch(xfinal(n)+[-1 -1 1 1 -1]*.4,[0 1 1 0 0]*pbase(n),'r')
    patch(xfinal(n)+[-1 -1 1 1 -1]*.2,[0 1 1 0 0]*bin_ct(n)*(1/DoTot),'b')
end
title('Expected Probability MASS Fn (PMF) (red); Experimental trials (blue)')
A = axis;
tt = text(.9*A(1)+.1*A(2), .1*A(3)+.9*A(4),'Hit enter to continue...');

input('Hit enter to see corresponding PDF (DENSITY function). ')
set(tt,'String',' ')
xc=xfinal;
xcp = min(xfinal)*1.5:.01:max(xfinal)*1.5;
ycp = normpdf(xcp,0,sqrt(Nlevel/8)*sqrt(2));
pp=plot(xcp,ycp,'k--','LineWidth',4,'Color',[0 .7 0]);
%set(pp,'YData',ycp)
figure(hF2) % make sure figure "jumps" into view again, as active...
