%  Example 3.4.  Generates sample functions, plots the analytic 
%  correlation function, and plots the analytic spectral density.
clear
%  Input the time sample vector.
t=0:.1:20;
%  The sample functions can be generated by passing white noise
%  through a lowpass filter.
%  The pitch stick position:
%  Input the lowpass filter parameters.
b1=sqrt(8);a1=-1;
%    Generate a discrete-time lowpass filter.
[phi1,del1]=c2d(a1,b1,0.1);
%    Generate the input to this lowpass filter.
u1=sqrt(10)*randn(size(t));
%    Perform the filtering to yield the sample function.
up=filter(del1,[1 -phi1],u1);
%  The roll stick position:
%  Input the lowpass filter parameters.
b2=1;a2=-sqrt(.25);
%    Generate a discrete-time lowpass filter.
[phi2,del2]=c2d(a2,b2,0.1);
%    Generate the input to this lowpass filter.
u2=sqrt(10)*randn(size(t));
%    Perform the filtering to yield the sample function.
ur=filter(del2,[1 -phi2],u2);
%  Plot these sample functions.
set(0,'DefaultAxesFontName','times')
set(0,'DefaultAxesFontSize',16)
set(0,'DefaultTextFontName','times')
x=[up;ur];
figure(1)
clf
plot(t,x(1,:),'r-',t,x(2,:),'r--')
xlabel('Time (sec)')
ylabel('Amplitude')
grid
%  Generate the correlation function using the analytic result.
tau=-10:0.1:10;
R11=4*exp(-abs(tau));
R22=exp(-0.5*abs(tau));
%  Plot the correlation function.
figure(2)
clf
plot(tau,R11,'r-',tau,R22,'r--')
xlabel('Time difference (sec)')
ylabel('Correlation')
grid
%  Generate the spectral density using the analytic result.
w=logspace(-2,1);
Sw11=8*ones(size(w))./(w.*w+ones(size(w)));
Sw22=ones(size(w))./(w.*w+0.25*ones(size(w)));
%  Plot the spectral density.
figure(3)
clf
loglog(w,Sw11,'r-',w,Sw22,'r--')
xlabel('Frequency (rad/sec)')
ylabel('Spectral density')
grid
