%
% This function illustrates the Levy's construction of the standard
% Brownian motion with the "midpoint displacement technique". One iteration
% of the function, doubles the set of points at which the BM is evaluated.
%
% Example of use:
%
% x = randn(1);
% figure(1);
% clf;
% gif_file = 'BM_simulation.gif';
% for i =1:15,
%  t = [0:length(x)]/length(x);
%  if i>1,
%    %delete(h);
%    set(h,'Color',[0.8,.8,0.8]);
%    set(h2,'Color',[0.8,.8,0.8]);
%  else
%    subplot(2,1,1);hold off; subplot(2,1,2);hold off;
%  end;
%  subplot(2,1,1);
%  %plot(t,[0;x],'.b'); hold on;
%  h = plot(t,[0;x],'r'); hold on;
%  set(gca,'fontsize',14);
%  title('Standard Brownian motion: Levy midpoint displacement');
%  xlabel('time');
%  subplot(2,1,2);
%  h2 = plot(t,[0,cumsum(randn(1,length(x)))]/sqrt(length(x)),'r');
%   set(gca,'fontsize',14);
%  title('Standard Brownian motion: Cumulative sums of white noise');
%  xlabel('time');
%  hold on;
%  pause;
%  drawnow
%  frame = getframe(1);
%  im = frame2im(frame);
%  [imind,cm] = rgb2ind(im,256);
%  if i ==1,
%    imwrite(imind,cm,gif_file,'gif','Loopcount',inf);
%  else
%   imwrite(imind,cm,gif_file,'gif','WriteMode','Append');
%  end;
%  x = Levy_BM(x);
% end;
% 
function y = Levy_BM(x);

n = length(x);
y = zeros(2*n,1);
y(2:2:(2*n)) = x;
y(1:2:(2*n)) = ([0;x(1:(n-1))]+x)/2+ randn(n,1)/sqrt(2*n*2);
