function [x, pest] = bern( N, p, L);
%function [x, pest] = bern( N, p, L);
% puts out an NxL matrix of iid Bern(p)

if(nargin == 2), L=1;end

x = rand(N,L) < 1-p;
pest = sum(x)/N;

x1 = x(:,1); %for illustration
ind0 = find( x1 == 0);% indexes the 0's
ind1 = find( x1 == 1);% indexes the 1's
n = [1:N]';
n0 = n(ind0); x10 = x1(ind0);
n1 = n(ind1); x11 = x1(ind1);
plot( n0,x10,'g *', n1,x11,'r x');grid
title(' Bernoulli Process ');
xlabel(' Time (msec) ');
ylabel(' x ');