% This Latex file contains the MATLAB programs and code segments listed in the book 
% "Intuitive Probability and Random Processes using MATLAB" by Steven M. Kay and published
% by Springer in 2005. The MATLAB programs and code segments are included between the 
% \begin{verbatim} and\end{verbatim} commands. All other lines are Latex commands and should 
% be deleted before MATLAB is run.   The MATLAB programs and code segments can be extracted 
% from this file using any text editing program such as the MATLAB editor or Word, as examples.
% This entire file can be compiled using Latex and printed out if desired to obtain a listing 
% of the programs and code segments.  
%
% This software is provided to the user free of charge.  The author and the publisher make no 
% warranty as to its correctness or effectiveness.
%
\documentclass[11pt]{book}
\textwidth 5.5in
\textheight 8in
\begin{document}
\setcounter{chapter}{22}
\chapter{MATLAB Code}


{\bf Chapter  1}\\
\hrulefill
\begin{verbatim}
              number=0;
              for i=1:4 % set up simulation for 4 coin tosses
                if rand(1,1)<0.75 % toss coin with p=0.75
                  x(i,1)=1; % head
                else
                  x(i,1)=0; % tail
                end
              number=number+x(i,1); % count number of heads
              end
\end{verbatim}
{\bf Chapter  2}\\
\begin{verbatim}
                      for i=1:M
                        u=rand(1,1);
                        if u<=p1
                          x(i,1)=x1;
                        elseif u>p1 & u<=p1+p2
                          x(i,1)=x2;
                        elseif u>p1+p2
                          x(i,1)=x3;
                        end
                      end
\end{verbatim}
\hrulefill
\begin{verbatim}
                       for i=1:M
                         x(i,1)=randn(1,1);
                       end
\end{verbatim}
\hrulefill
\begin{verbatim}
randn('state',0)
x=randn(1000,1);
bincenters=[-3.5:0.5:3.5]';
bins=length(bincenters);
h=zeros(bins,1);
for i=1:length(x)
  for k=1:bins
    if x(i)>bincenters(k)-0.5/2 ...
     & x(i)<=bincenters(k)+0.5/2
      h(k,1)=h(k,1)+1;
    end
  end   
end
pxest=h/(1000*0.5);
xaxis=[-4:0.01:4]';
px=(1/sqrt(2*pi))*exp(-0.5*xaxis.^2);
\end{verbatim}
\hrulefill
\begin{verbatim}
       randn('state',0)
       M=100;count=0;
       x=randn(M,1);
       for i=1:M
         if x(i)>2
           count=count+1;
         end
       end
       probest=count/M
\end{verbatim}
\hrulefill
\begin{verbatim}
     randn('state',0)
     M=100;
     meanest=0;
     x=randn(M,1);
     for i=1:M
       meanest=meanest+(1/M)*x(i);
     end
     meanest
\end{verbatim}
\hrulefill
\begin{verbatim}
  randn('state',0)
  M=100;
  meanest=0;
  x=randn(M,1);
  for i=1:M
    meanest=meanest+(1/M)*x(i)^2;
  end
  meanest
\end{verbatim}
\hrulefill
\begin{verbatim}
                       A=[0.1:0.1:5]';
                       for k=1:length(A)
                         error=0;
                         for i=1:1000
                           w=randn(1,1);
                           if A(k)/2+w<=0
                             error=error+1;
                           end
                         end
                         Pe(k,1)=error/1000;
                       end
\end{verbatim}
\hrulefill
\begin{verbatim}
%  matlabexample.m 
%
%  This program computes and plots samples of a sinusoid
%  with amplitudes 1, 2, and 4.  If desired, the sinusoid can be
%  clipped to simulate the effect of a limiting device.
%  The frequency is 1 Hz and the time duration is 10 seconds. 
%  The sample interval is 0.1 seconds. The code is not efficient but
%  is meant to illustrate MATLAB statements.
%
clear all % clear all variables from workspace
delt=0.01; % set sampling time interval
F0=1; % set frequency
t=[0:delt:10]'; % compute time samples 0,0.01,0.02,...,10
A=[1 2 4]'; % set amplitudes
clip='yes'; % set option to clip
for i=1:length(A) % begin computation of sinusoid samples
   s(:,i)=A(i)*cos(2*pi*F0*t+pi/3); % note that samples for sinusoid 
                                    % are computed all at once and  
                                    % stored as columns in a matrix                                      
   if clip=='yes' % determine if clipping desired
      for k=1:length(s(:,i)) % note that number of samples given as 
                             % dimension of column using length command
         if s(k,i)>3 % check to see if sinusoid sample exceeds 3
            s(k,i)=3; % if yes, then clip
         elseif s(k,i)<-3 % check to see if sinusoid sample is less 
            s(k,i)=-3; % than -3 if yes then clip
         end
      end
   end
end
figure % open up a new figure window
plot(t,s(:,1),t,s(:,2),t,s(:,3)) % plot sinusoid samples versus time 
                                 % samples for all three sinusoids
grid % add grid to plot
xlabel('time, t') % label x-axis
ylabel('s(t)') % label y-axis
axis([0 10 -4 4]) % set up axes using axis([xmin xmax ymin ymax])
legend('A=1','A=2','A=4') % display a legend to distinguish 
                          % different sinusoids
\end{verbatim}
{\bf Chapter  3}\\
\hrulefill
\begin{verbatim}
% birthday.m
%
clear all
rand('state',0)
BD=[0:365]';
event=zeros(10000,1); % initialize to no successful events
for ntrial=1:10000
for i=1:23
   x(i,1)=ceil(365*rand(1,1)); % chooses birthdays at random  
                               % (ceil rounds up to nearest integer)
end
y=sort(x); % arranges birthdays in ascending order
z=y(2:23)-y(1:22); % compares successive birthdays to each other
w=find(z==0); % flags same birthdays
if length(w)>0
   event(ntrial)=1; % event occurs if one or more birthdays same
end
end
prob=sum(event)/10000
\end{verbatim}
{\bf Chapter  4}\\
{\bf Chapter  5}\\
\hrulefill
\begin{verbatim}
                      for i=1:M
                        u=rand(1,1);
                        if u<=0.2
                          x(i,1)=1;
                        elseif u>0.2 & u<=0.8
                          x(i,1)=2;
                        elseif u>0.8
                          x(i,1)=3;
                        end
                      end
\end{verbatim}
{\bf Chapter  6}\\
\hrulefill
\begin{verbatim}
%  PMFdata.m
%
%  This program generates the outcomes for N trials 
%  of an experiment for a discrete random variable.
%  Uses the method of Section 5.9.  
%  It is a function subprogram.
%
%  Input parameters:
%
%    N   - number of trials desired
%    xi  - values of x_i's of discrete random variable (M x 1 vector)
%    pX  - PMF of discrete random variable (M x 1 vector)
%
%  Output parameters:
%
%    x   - outcomes of N trials (N x 1 vector)
%
function x=PMFdata(N,xi,pX)
M=length(xi);M2=length(pX);
if M~=M2
   message='xi and pX must have the same dimension'
end
for k=1:M ; % see Section 5.9 and Figure 5.14 for approach used here
   if k==1
      bin(k,1)=pX(k); % set up first interval of CDF as [0,pX(1)]
  
     else
        bin(k,1)=bin(k-1,1)+pX(k); % set up succeeding intervals
                                   % of CDF
     end
  end  
  u=rand(N,1); % generate N outcomes of uniform random variable
  for i=1:N % determine which interval of CDF the outcome lies in
            % and map into value of xi
     if u(i)>0&u(i)<=bin(1)  
           x(i,1)=xi(1);
        end
        for k=2:M
           if u(i)>bin(k-1)&u(i)<=bin(k)
              x(i,1)=xi(k);
           end
        end
     end

\end{verbatim}
{\bf Chapter  7}\\
\hrulefill
\begin{verbatim}
                    
                    for m=1:M
                      u=rand(1,1);
                      if u<=1/8
                        x(m,1)=0;y(m,1)=0;
                      elseif u>1/8&u<=1/4
                        x(m,1)=0;y(m,1)=1;
                      elseif u>1/4&u<=1/2
                        x(m,1)=1;y(m,1)=0;
                      else
                        x(m,1)=1;y(m,1)=1;
                      end
                    end
\end{verbatim}
{\bf Chapter  8}\\
\begin{verbatim}
 for m=1:M
   if rand(1,1)<0.5
     y(m,1)=PMFdata(1,[1 2 3 4 5 6]',[1/6 1/6 1/6 1/6 1/6 1/6]');
   else
      y(m,1)=PMFdata(1,[2 3]',[1/2 1/2]');
   end
 end
\end{verbatim}
\hrulefill
\begin{verbatim}
                    
                    for m=1:M
                      ux=rand(1,1);
                      uy=rand(1,1);
                      if ux<=1/4; % Refer to px[i]
                        x(m,1)=0;
                        if uy<=1/2 % Refer to py|x[j|0]
                          y(m,1)=0;
                        else
                          y(m,1)=1;
                        end
                      else
                        x(m,1)=1;  % Refer to px[i]
                        if uy<=1/3 % Refer to py|x[j|1]
                          y(m,1)=0;
                        else
                          y(m,1)=1;
                        end
                      end  
                    end
\end{verbatim}
{\bf Chapter  9}\\
\hrulefill
\begin{verbatim}
%  covexample.m
clear all % clears out all previous variables from workspace
rand('state',0); % sets random number generator to initial value
M=1000;
for m=1:M % generate realizations of X (see Section 7.11)
   u=rand(1,1);
   if u<=0.25
      x(1,m)=-8;x(2,m)=0;
   elseif u>0.25&u<=0.5
      x(1,m)=0;x(2,m)=-8;
   elseif u>0.5&u<=0.75
      x(1,m)=2;x(2,m)=6;
   else
      x(1,m)=6;x(2,m)=2;
   end
end
meanx=[0 0]'; % estimate mean vector of X
for m=1:M
   meanx=meanx+x(:,m)/M;
end
meanx
CX=zeros(2,2);
for m=1:M % estimate covariance matrix of X
   xbar(:,m)=x(:,m)-meanx;
   CX=CX+xbar(:,m)*xbar(:,m)'/M;
end
CX
A=[1/sqrt(2) -1/sqrt(2);1/sqrt(2) 1/sqrt(2)];
for m=1:M % transform random vector X
   y(:,m)=A*x(:,m);
end
meany=[0 0]'; %estimate mean vector or Y
for m=1:M
   meany=meany+y(:,m)/M;
end
meany
CY=zeros(2,2);
for m=1:M % estimate covariance matrix of Y
   ybar(:,m)=y(:,m)-meany;
   CY=CY+ybar(:,m)*ybar(:,m)'/M;
end
CY
\end{verbatim}
{\bf Chapter  10}\\
\hrulefill
\begin{verbatim}
                        for i=1:M
                        u=rand(1,1);
                          if u<=0.75
                            x(i,1)=1; % head mapped into 1
                          else
                            x(i,1)=0; % tail mapped into 0
                          end
                        end
\end{verbatim}  
\hrulefill
\begin{verbatim}
         %  Assume outcomes are in x, which is M x 1 vector
         M=1000;
         bincenters=[-4:0.5:4]'; % set binwidth = 0.5
         bins=length(bincenters);
         h=zeros(bins,1);
         for i=1:length(x)  % count outcomes in each bin
           for k=1:bins
             if x(i)>bincenters(k)-0.5/2. ...
               & x(i)<bincenters(k)+0.5/2
               h(k,1)=h(k,1)+1;
             end
           end   
         end
         pxest=h/(M*0.5); % see (10.38)
\end{verbatim}
\hrulefill
\begin{verbatim}
%  Q.m
%
%  This program computes the right-tail probability 
%  (complementary cumulative distribution function) for 
%  a N(0,1) random variable.
%
%  Input Parameters:
%
%    x - Real column vector of x values
%
%  Output Parameters:
%
%    y - Real column vector of right-tail probabilities
%
%  Verification Test Case:
%
%  The input x=[0 1 2]'; should produce y=[0.5 0.1587 0.0228]'.
%
   function y=Q(x)
   y=0.5*erfc(x/sqrt(2));  % complementary error function
\end{verbatim}
\hrulefill\\
\begin{verbatim}
%  Qinv.m
%  
%  This program computes the inverse Q function or the value 
%  which is exceeded by a N(0,1) random variable with a 
%  probability of x.
%
%  Input Parameters:
%
%    x - Real column vector of right-tail probabilities
%        (in interval [0,1])
%
%  Output Parameters:
%
%    y - Real column vector of values of random variable
%
%  Verification Test Case:
%  
%  The input x=[0.5 0.1587 0.0228]'; should produce 
%  y=[0 0.9998 1.9991]'.
%
   function y=Qinv(x)
   y=sqrt(2)*erfinv(1-2*x);  % inverse error function
\end{verbatim}
{\bf Chapter  11}\\
\begin{verbatim}
rand('state',0) % sets random number generator to 
                % initial value
M=10000;gamma=5;% change M for different estimates
u=rand(M,1);    % generates M U(0,1) realizations
x=-log(1-u);    % generates M exp(1) realizations
k=0;
for i=1:M       % computes estimate of P[X>gamma]
   if x(i)>gamma
      k=k+1;
      y(k,1)=(1/sqrt(2*pi))*exp(-0.5*x(i)^2+x(i)); % computes weights 
                                                   % for estimate
    end
end
Qest=sum(y)/M  % final estimate of P[X>gamma]
\end{verbatim}
{\bf Chapter  12}\\
\begin{verbatim}
randn('state',0) % set random number generator to initial value
G=[1 0;0.9 sqrt(1-0.9^2)]; % define G matrix
M=2000; % set number of realizations
for m=1:M
  x=randn(1,1);y=randn(1,1); % generate realizations of two independent 
                             % N(0,1) random variables
  wz=G*[x y]'+[1 1]'; % transform to desired mean and covariance
  WZ(:,m)=wz; % save realizations in 2 x M array
end
Wmeanest=mean(WZ(1,:)); % estimate mean of W
Zmeanest=mean(WZ(2,:)); % estimate mean of Z
WZbar(1,:)=WZ(1,:)-Wmeanest; % subtract out mean of W
WZbar(2,:)=WZ(2,:)-Zmeanest; % subtract out mean of Z
Cest=[0 0;0 0];
for m=1:M
  Cest=Cest+(WZbar(:,m)*WZbar(:,m)')/M; % compute estimate of 
                                        % covariance matrix
end
Wmeanest % write out estimate of mean of W
Zmeanest % write out estimate of mean of Z
Cest % write out estimate of covariance matrix
\end{verbatim}
{\bf Chapter  13}\\
\begin{verbatim}
randn('state',0) % set random number generator to initial value
rho=0.9;
M=500; % set number of realizations to generate
for m=1:M
  x(m,1)=randn(1,1); % generate realization of N(0,1) random 
                     % variable (Step1)
  ygx(m,1)=rho*x(m)+sqrt(1-rho^2)*randn(1,1); % generate 
                                              % Y|(X=x) (Step 2)
end
\end{verbatim}
{\bf Chapter  14}\\
\begin{verbatim}
        C=[1 2/3 1/3;2/3 1 2/3;1/3 2/3 1];
        G=chol(C)'; % perform Cholesky decomposition
                    % MATLAB produces C=A'*A so G=A'
        M=200;
        for m=1:M % generate realizations of x
          u=[randn(1,1) randn(1,1) randn(1,1)]';
          x(:,m)=G*u; % realizations stored as columns of 3 x 200 matrix
        end
\end{verbatim}
{\bf Chapter  15}\\
\begin{verbatim}
%  This program demonstrates the central limit theorem.  It determines the
%  PDF for the sum S_N of N IID random variables. Each marginal PDF is assumed
%  to be nonzero over the interval (0,1). The repeated convolution integral is 
%  implemented using a discrete convolution.  The plots of the PDF of S_N as
%  N increases are shown successively (press carriage return for next plot).
%
% clt_demo.m
clear all
delu=0.005;
u=[0:delu:1-delu]'; % p_X defined on interval [0,1)
p_X=ones(length(u),1); % try p_X=abs(2-4*u) for really strange PDF
x=[u;u+1]; % increase abcissa values since repeated 
           % convolution increases nonzero width of output
p_S=zeros(length(x),1);
N=12; % number of random variables summed
for j=1:length(x) % start discrete convolution approximation
                  % to continuous convolution
  for i=1:length(u)
    if j-i>0&j-i<=length(p_X)
      p_S(j)=p_S(j)+p_X(i)*p_X(j-i)*delu;
    end  
  end
end
plot(x,p_S) % plot results for N=2 
grid
axis([0 N 0 1]) % set axes lengths for plotting
xlabel('x')
ylabel('p_S')
title('PDF for S_N')
text(0.75*N,0.85,'N = 2') % label plot with the 
                          % number of convolutions
for n=3:N  
  pause
  x=[x;u+n-1]; % increase abcissa values since 
               % repeated convolution increases
               % nonzero width of output
  p_S=[p_S;zeros(length(u),1)];
  g=zeros(length(p_S),1);
  for j=1:length(x) % start discrete convolution
  for i=1:length(u)
    if j-i>0
      g(j,1)=g(j,1)+p_X(i)*p_S(j-i)*delu;
    end  
  end
end
p_S=g; % plot results for N=3,4,...,12
plot(x,p_S) 
grid
axis([0 N 0 1])
xlabel('x')
ylabel('p_S')
title('PDF for S_N')
text(0.75*N,0.85,['N = ' num2str(n)])
end
\end{verbatim}
{\bf Chapter  16}\\
\begin{verbatim}
randn('state',0)
N=51;
x=randn(N,1)+0.1*[0:N-1]'; % for Figure 16.7a
y=sqrt(0.95.^[0:50]').*randn(N,1); % for Figure 16.7b
\end{verbatim}
\hrulefill\\
\begin{verbatim}
randn('state',0)
u=randn(21,1);
for i=1:21
  if i==1
    x(i,1)=0.5*(u(1)+randn(1,1)); % needed to initialize sequence
  else
    x(i,1)=0.5*(u(i)+u(i-1));
  end
end
\end{verbatim}
\hrulefill\\
\begin{verbatim}
clear all
rand('state',0)
n=[0:31]';
nreal=50;
for i=1:nreal
   x(:,i)=cos(2*pi*0.1*n+2*pi*rand(1,1));
end
plot(n,x(:,1),'.')
grid
hold on
for i=2:nreal
   plot(n,x(:,i),'.')
end
axis([0 31 -1.5 1.5])
\end{verbatim} 
\hrulefill\\
\begin{verbatim}
clear all
randn('state',0)
years=[1895:2002]';
N=length(years);
n=[0:N-1]';
A=[N sum(n);sum(n) sum(n.^2)]; % precompute matrix (see (16.9))
B=inv(A); % invert matrix
for i=1:20
  xn=9.76+sqrt(10.05)*randn(N,1); % generate realizations
  baest=B*[sum(xn);sum(n.*xn)]; % estimate a and b using (16.9)
  aest=baest(2);best=baest(1);
  meanest(:,i)=aest*n+best; % determine mean sequence estimate
end
figure % plot mean sequence estimates and overlay
plot(n,meanest(:,1))
grid
xlabel('n')
ylabel('Estimated mean')
axis([0 107 5 15])
hold on
for i=2:20
  plot(n,meanest(:,i))
end
\end{verbatim}
{\bf Chapter  17}\\
\begin{verbatim}
clear all
randn('state',0)
a1=0.25;a2=0.98;
varu1=1-a1^2;varu2=1-a2^2;
varx1=varu1/(1-a1^2);varx2=varu2/(1-a2^2); % this is r_X[0]
x1(1,1)=sqrt(varx1)*randn(1,1); % set initial condition X[-1] 
                                % see Problems 17.17, 17.18
x2(1,1)=sqrt(varx2)*randn(1,1);
for n=2:31
   x1(n,1)=a1*x1(n-1)+sqrt(varu1)*randn(1,1);
   x2(n,1)=a2*x2(n-1)+sqrt(varu2)*randn(1,1);
end
\end{verbatim}
\hrulefill\\
\begin{verbatim}
n=[0:30]';N=length(n);
a1=0.25;a2=0.98;
varu1=1-a1^2;varu2=1-a2^2;
r1true=(varu1/(1-a1^2))*a1.^n; % see (17.21)
r2true=(varu2/(1-a2^2))*a2.^n;
for k=0:N-1
   r1est(k+1,1)=(1/(N-k))*sum(x1(1:N-k).*x1(1+k:N));
   r2est(k+1,1)=(1/(N-k))*sum(x2(1:N-k).*x2(1+k:N));
end
\end{verbatim}
\hrulefill\\
\begin{verbatim}
Nfft=1024; % set FFT size
Pav1=zeros(Nfft,1);Pav2=Pav1; % set up arrays with desired dimension
f=[0:Nfft-1]'/Nfft-0.5; % set frequencies for later plotting 
                        % of PSD estimate
for i=0:I-1
   nstart=1+i*L;nend=L+i*L; % set up beginning and end points 
                            % of ith block of data
   y1=x1(nstart:nend);
   y2=x2(nstart:nend);
% take FFT of block, since FFT outputs samples of  Fourier 
% transform over frequency range [0,1), must shift FFT outputs 
% for [1/2,1) to [-1/2, 0), then take complex magnitude-squared, 
% normalize by L and average
Pav1=Pav1+(1/(I*L))*abs(fftshift(fft(y1,Nfft))).^2; 
Pav2=Pav2+(1/(I*L))*abs(fftshift(fft(y2,Nfft))).^2;
end
\end{verbatim}
{\bf Chapter  18}\\
\begin{verbatim}
clear all
randn('state',0)
a=0.9;varu=0.5;vars=varu/(1-a^2);varw=1;N=50; % set up parameters
for n=0:N-1 % generate signal realization
   nn=n+1;
   if n==0 % use Gaussian random processes
      s(nn,1)=sqrt(vars)*randn(1,1); % initialize first sample
                                     % to avoid transient
   else
      s(nn,1)=a*s(nn-1)+sqrt(varu)*randn(1,1);
   end
end
x=s+sqrt(varw)*randn(N,1); % add white Gaussian noise
Nfft=1024; % set up FFT length
Ps=varu./(abs(1-a*exp(-j*2*pi*[0:Nfft-1]'/Nfft)).^2); % compute PSD of
                                                      % signal, frequency 
                                                      % interval is [0,1]
Hf=Ps./(Ps+varw); % form Wiener smoother
sestf=Hf.*fft(x,Nfft);  % signal estimate in frequency domain, 
                        % frequency interval is [0,1]
sest=real(ifft(sestf,Nfft)); % inverse Fourier transform
\end{verbatim}
\hrulefill\\
\begin{verbatim}
N=length(xseg); % xseg is the data shown in Figure 18.11
Nfft=1024; % set up FFT length for Fourier transforms
freq=[0:Nfft-1]'/Nfft-0.5; % PSD frequency points to be plotted
P_per=(1/N)*abs(fftshift(fft(xseg,Nfft))).^2; % compute periodogram
p=12; % dimension of autocorrelation matrix
for k=1:p+1 % estimate ACS for k=0,1,...,p (MATLAB indexes must start at 1)
   rX(k,1)=(1/N)*sum(xseg(1:N-k+1).*xseg(k:N));
end
r=rX(2:p+1); % fill in right-hand-side vector
for i=1:p % fill in autocorrelation matrix
   for j=1:p
      R(i,j)=rX(abs(i-j)+1);
   end
end
a=inv(R)*r; % solve linear equations to find AR filter parameters
varu=rX(1)-a'*r; % find excitation noise variance
den=abs(fftshift(fft([1;-a],Nfft))).^2; % compute denominator of AR PSD
P_AR=varu./den; % compute AR PSD
\end{verbatim}
{\bf Chapter  19}\\
\begin{verbatim}
% assume realizations are x[n], y[n] for n=1,2,...,N
for k=0:M % compute zero and positive lags, see (19.46)
  rxypos(k+1,1)=(1/(N-k))*sum(x(1:N-k).*y(1+k:N)); % compute values 
                                                   % for k=0,1,...,M
end
for k=1:M % compute negative lags, see (19.46)
  rxyneg(k+1,1)=(1/(N-k))*sum(x(k+1:N).*y(1:N-k)); % compute values 
                                                   % for k=-M,-(M-1),...,-1
end
rxy=[flipud(rxyneg(2:M+1,1));rxypos]; % arrange values from k=-M to k=M
\end{verbatim}
{\bf Chapter  20}\\
\begin{verbatim}
Fs=100;  % set sampling rate for later plotting
L=50;I=20; % L = length of block, I = number of blocks 
n=[0:I*L-1]'; % set up time indices
Nfft=1024; % set FFT length for Fourier transform
Pav=zeros(Nfft,1);
f=[0:Nfft-1]'/Nfft-0.5; % set discrete-time frequencies
for i=0:I-1
   nstart=1+i*L;nend=L+i*L; % set start and end time indices
                            % of block
   y=x(nstart:nend); % extract block of data
Pav=Pav+(1/(I*L))*abs(fftshift(fft(y,Nfft))).^2; 
                                        % compute periodogram 
                                        % and add to average 
                                        % of periodograms
end
F=f*Fs; % convert to continuous-time (analog) frequency in Hz
Pest=Pav/Fs; % convert discrete-time PSD to continuous-time PSD
plot(F,Pest)
\end{verbatim}
\hrulefill\\
\begin{verbatim}
clear all
rand('state',0)
t=[0:0.01:0.99]'; %  set up transmit pulse time interval
F0=10;
s=cos(2*pi*F0*t); %  transmit pulse
ss=[s;zeros(1000-length(s),1)];  % put transmit pulse in receive window
tt=[0:0.01:9.99]';  % set up receive window time interval
x=zeros(1000,1);  
for i=1:100  % add up all echos, one for each 0.1 sec interval
  tau=round(10*i+10*(rand(1,1)-0.5));  % time delay for each 0.1 sec interval
                                       % is uniformly distributed - round
                                       % time delay to integer 
  x=x+rand(1,1)*shift(ss,tau);
end
\end{verbatim}
\hrulefill\\
\begin{verbatim}
% shift.m
%
function y=shift(x,Ns)
%
%  This function subprogram shifts the given sequence by Ns points.
%  Zeros are shifted in either from the left or right.
%
%  Input parameters:
%    x  - array of dimension Lx1
%    Ns - integer number of shifts where Ns>0 means a shift to the
%         right and Ns<0 means a shift to the left and if Ns=0, then
%         the sequence is not shifted
%
%   Output parameters:
%    y  - array of dimension Lx1 containing the 
%         shifted sequence
L=length(x);
if abs(Ns)>L
   y=zeros(L,1);
   else  
if Ns>0
y(1:Ns,1)=0;
y(Ns+1:L,1)=x(1:L-Ns);
elseif Ns<0
   y(L-abs(Ns)+1:L,1)=0;
   y(1:L-abs(Ns),1)=x(abs(Ns)+1:L);
else
   y=x;
end
end
\end{verbatim}
{\bf Chapter  21}\\
\begin{verbatim}
clear all
rand('state',0)
lambda=2; % set arrival rate
T=5; % set time interval in seconds
for i=1:1000
   z(i,1)=(1/lambda)*log(1/(1-rand(1,1))); % generate interarrival times
   if i==1 % generate arrival time
      t(i,1)=z(i);
      else
         t(i,1)=t(i-1)+z(i,1);
      end
      if t(i)>T % test to see if desired time interval has elapsed
         break
      end
   end
   M=length(t)-1; % number of arrivals in interval [0,T]
   arrivals=t(1:M); % arrival times in interval [0,T]
\end{verbatim}
{\bf Chapter  22}\\
\begin{verbatim}
clear all
rand('state',0)
N=1000;  % set number of samples desired
p0=[1/3 1/3 1/3]'; % set initial probability vector
P=[4/8 3/8 1/8;3/8 2/8 3/8;1/8 3/8 4/8]; % set transition prob. matrix
xi=[0 1 2]'; % set values of PMF
X0=PMFdata(1,xi,p0); % generate X[0] (see Appendix 6B for PMFdata.m 
                     % function subprogram)
i=X0+1; % choose appropriate row for PMF
X(1,1)=PMFdata(1,xi,P(i,:)); % generate X[1]
i=X(1,1)+1; % choose appropriate row for PMF
for n=2:N % generate X[n]
   i=X(n-1,1)+1; % choose appropriate row for PMF
   X(n,1)=PMFdata(1,xi,P(i,:)); 
end
\end{verbatim}
\hrulefill\\
\begin{verbatim}
% sierpinski.m
%
clear all
rand('state',0)
r(:,1)=[0 0]'; % set up reference points
r(:,2)=[100 0]';
r(:,3)=[50 100]';
x0=[10 80]'; % set initial state
plot(x0(1),x0(2),'.') % plot state outcome as point
axis([0 100 0 100])
hold on
xn_1=x0;
for n=1:10000 % generate states
   j=floor(3*rand(1,1)+1); % choose at random one of three 
                           % reference points
   xn=round(0.5*(r(:,j)+xn_1)); % generate new state
   plot(xn(1),xn(2),'.') % plot state outcome as point
   xn_1=xn; % make current state the previous one for 
            % next transition
end
grid
hold off
\end{verbatim}
\end{document}