MATLAB and C Programs - Digital Signal Processing: A Practitioner's Approach - page 140

omega=2*pi*f;
b=omega*omega;
a=[1 2*delta*omega omega*omega];
Ts=1/(10*f);
sysc=tf(b,a);
sysd=c2d(sysc,Ts,'tustin');
[b1 a1]=tfdata(sysd,'v');
y1=filter(b1,a1,u);
%pause;
[x,y]=pdf_y(y1,50);
[sy1,f1]=spk(y1,1);
gain=mean(sy1(1:150));
N=length(f1); m=abs(freqz(b1,a1,N));
sigma=1.5;
%sigma=12;
xx=(x.*x)/(2*sigma);
yhat=exp(-xx)/(sqrt(2*pi*sigma));
scale=sum(yhat); yhat=yhat/scale;
subplot(222);plot(x,y,'o',x,y,x,yhat,'.-');grid;
ylabel(' PDF of yk ');title(' B ');
%title('Fig 2.12')
subplot(221);plot(f1,sy1,'.',f1,m*gain,'.-');grid;
ylabel(' Coloured Noise yk ');title(' A ' );
%title('Fig 2.13')
%xlabel(' Coloured Gaussian Noise - Ouput')
%print -depsc f21_2;pause;
%print -depsc f21_3
A.2.7 Program f 21_67
%Digital Signal Processing:A Practitioner's Approach
%Dr.Kaluri Venkata Ranga Rao, [e-mail]
%Generates figures 2.17,2.18
clear;close
sigma_x=0.2;
sigma_y=0.2;
x_val=-5;y_val=7;
row_max=8;
A=fix((rand([row_max 2])-0.5)*20);
x1=x_val-randn( [1 row_max])*sigma_x;
x2=y_val-randn([1 row_max])*sigma_y;
c1=A(:,1); p1=c1.*x1';
c2=A(:,2); p2=c2.*x2';