Matlab Code for Monson H. Hayes Statistical DSP Computer Exercise C3.4

Matlab Code for Monson H. Hayes Statistical DSP Computer Exercise C3.4 of chapter 3.

% Matlab code for the above problem
clc;
clear all; 
close all;
a(1) = 0;
a(2) = -0.81;
b(1) = 1;       % b(0) = 1 for the purpose of MATLAB indexing
x(1) = 1;       % x(0) = 1
x(2) = 0;       % x(1) = 2
x(3) = 0.81;    % x(2) = 0.81

for n= 4:24
    x(n) = a(1) * x(n-1) + a(2) * x(n-2);
end

CoefM = [1, a(1), a(2); a(1), 1, a(1); a(2), a(1), 1];  %Coefficient matrix
C0 = [1;0;0];
rx = inv(CoefM) * C0;   % rx(1), rx(2), rx(3)

for k = 4:24
    rx(k) = -a(1)*rx(k-1) - a(2)*rx(k-2);   % [rx(4), rx(24)]
end

rxHat = autocorr(x,23);      % Estimated Autocorrelation

px = fft(rx);                         % True PSD
Pxt=abs(px);                          % Magnitude of Estimated power

Pt = fft(rxHat);                      % Estimated PSD
Pxc = abs(Pt);                        % Magnitude of Estimated PSD

d = [rx(1), rx(2); rx(2), rx(1)];
s = -[rx(2); rx(3)];
g = inv(d) * s;                       % Computing a(1) and a(2) using Yule-Walker Equation
b0 = rx(1) + g(1)*rx(2) + g(2)*rx(3); % Computing b(0) using Yule-Walker Equation

syms w
pxe = b0^2/(abs(1+g(1)*exp(-1i*w) + g(2)*exp(-1i*2*w)).^2);

disp('Coefficients found Using Yule-Walker Equations:');
disp('===============================================');
disp('a(1):'),disp(g(1));
disp('a(2):'),disp(g(2));
disp('b(0):'),disp(b0);

subplot(3,2,1), stem(x);
title('x(n)')
subplot(3,2,3), stem(rx);
title('True Autocorrelation')
subplot(3,2,5),autocorr(x,'NumLags',23);
title('Estimated Autocorrelation')
subplot(3,2,2),plot(Pxc);
title('Estimated Power Spectrum  part c')
subplot(3,2,4),fplot(pxe, [-2*pi, 2*pi]);
title('Estimated Power Spectrum  part e')
subplot(3,2,6),plot(Pxt);
title('True Power Spectrum');