Matlab Code for Computer Exercise C7.3 Monson H. Hayes Statistical DSP

In this post, Matlab Code for Computer Exercise C7.3 Monson H. Hayes Statistical Digital Signal Processing book and its simulation output is discussed. Please refer the question (C7.3) below.

Matlab code for computer exercise 7.3 (a)

clc;
clear;
close all;
%Define the length of the simulation.
N = input('Please enter the number of iteration for discrete Kalman Filter, N = ');  

%Transition matrix 
A = 0.8;  % State Transition Matrix 
C = 1.0;  % Obserbation Matrix

%Define the noise covariances.
Q = 0.36; % Variance of process noise
R = 1; % Variance of measurement noise
%% Generate AR(1) process
e = randn(100,1); 
b = [1 0]; 
X = filter(b,0.8,e);

%Calculate the process and measurement noise.
w1 = randn(1,N); %This line and the next can be commented out after running
v1 = randn(1,N); %the script once to generate multiple runs with identical
%noise (for better comparison).
w = w1*sqrt(Q);
v = v1*sqrt(R);

%Initial condition on the state, x.
x_0 = X(1);

%Initial guesses for state and a posteriori covariance.
p_0 = 1;

%Calculate the first estimates for all values based upon the initial guess
%of the state and the a posteriori covariance. The rest of the steps will
%be calculated in a loop.
 x(1)= A*x_0 + w(1); 
 y(1)= C*x(1)+ v(1);
 
%Predictor equations
x_pred(1) = A*x_0;
p_pred(1) = A*A*p_0 + Q;
residual(1) = y(1)- C*x_pred(1);

%Corrector equations
K(1)= C*p_pred(1)/(C*C*p_pred(1)+R);
p_est(1)= p_pred(1)*(1-C*K(1));
x_est(1)= x_pred(1)+K(1)*residual(1);

%Calculate the rest of the values.
for j=2:N
    %Calculate the state and the output

    x(j)= A*x(j-1)+ w(j);
    y(j) = C*x(j)+ v(j);
    
    %Predictor equations
    x_pred(j)= A*x_est(j-1);
    residual(j)= y(j)- C*x_pred(j);
    p_pred(j) = A*A*p_est(j-1)+ Q;
    %Corrector equations
    K(j) = C*p_pred(j)/(C*C*p_pred(j)+ R);
    p_est(j) = p_pred(j)*(1-C*K(j));
    x_est(j) = x_pred(j)+ K(j)*residual(j);
    x_pred(j) = x_est(j);
    
end
j=1:N;

disp('The Kalman gain for the given number of iterations is:');
disp('=====================================================');
disp(K(j));
plot(K);

Matlab code for computer exercise 7.3 (b)

clc;
clear;
close all;
%Define the length of the simulation.
nlen = 11;

%Define the system.  
A = 1;  % a=1 for a constant,
C = 1;

%Define the noise covariances.
Q = 1;
R = 0.64;
%% Generate AR(3) process
e2 = randn(100,1); 
b1 = [1 0]; 
y = filter(b1,[-0.1 -0.09 0.648],e2);
% initial estimates (guesses) for x
x_0 = y(1);
%Initial guesses for error covariance.
p_0 = 0.6;

%Calculate the process and measurement noise.
w1 = randn(1,nlen); %This line and the next can be commented out after running
v1 = randn(1,nlen); %the script once to generate multiple runs with identical
%noise (for better comparison).
w = w1*sqrt(Q);
v = v1*sqrt(R);

%Calculate the first estimates for all values based upon the initial guess
%of the state and the a posteriori covariance. The rest of the steps will
%be calculated in a loop.

%Calculate the state and the output
x(1) = A*x_0 + w(1);
z(1) = C*x(1) + v(1);

%Predictor equations
x_pred(1) = A*x_0;
p_pred(1) = A*A*p_0 + Q;
residual(1) = z(1)-C*x_pred(1);

%Corrector equations
k(1) = C*p_pred(1)/(C*C*p_pred(1) + R);
p_est(1) = p_pred(1)*(1-C*k(1));
x_est(1) = x_pred(1)+ k(1)*residual(1);

%Calculate the rest of the values.
for j = 2:nlen
    %Calculate the state and the output
    x(j)= A*x(j-1)+ w(j);
    z(j)= C*x(j)+ v(j);
    %Predictor equations
    x_pred(j)= A*x_est(j-1);
    residual(j) = z(j)-C*x_pred(j);
    p_pred(j) = A*A*p_est(j-1)+ Q;
    %Corrector equations
    k(j) = C*p_pred(j)/(C*C*p_pred(j)+ R);
    p_est(j) = p_pred(j)*(1-C*k(j));
    x_est(j) = x_pred(j)+k(j)*residual(j);
    x_pred(j) = x_est(j);
    x_pred(j) = x_est(j);
    %disp(k(j));
    end

%j=1:nlen;
disp('The Kalman gain for n = 0 to n = 10 is:');
disp('=======================================');
disp(k);
%% plot Kalman gain for 10 iterations
plot(k);
title('Kalman Gain for 11 iterations');
xlabel('Number of iterations');
ylabel('Gain');

Matlab code for computer exercise 7.3 (c)

clc;
clear all;
%Define the length of the simulation.
nlen = 10;

%Define the system.  
A = 0.8;  
H = 1;

%Define the noise covariances.
Q = 1;
R = 0.64;
%% Generate AR(3) process
phi = 0;
e = randn(100,1); 
b = [1 phi]; 
y = filter(b,[1.2728 -0.09 0.648],e);
% initial estimates (guesses) for
% 1) the state.
x_0 = y(1);

% 2) the  error covariance.
p_0 = 0;

%Calculate the process and measurement noise.
w1 = randn(1,nlen); %This line and the next can be commented out after running
v1 = randn(1,nlen); %the script once to generate multiple runs with identical
%noise (for better comparison).
w = w1*sqrt(Q);
v = v1*sqrt(R);

%Calculate the first estimates for all values based upon the initial guess
%of the state and the a posteriori covariance. The rest of the steps will
%be calculated in a loop.
%
%Calculate the state and the output
x(1) = A*x_0 + w(1);
z(1) = H*x(1) + v(1);
%Predictor equations
x_pred(1) = A*x_0;
p_pred(1) = A*A*p_0 + Q;
residual(1) = z(1)-H*x_pred(1);

%Corrector equations
k(1) = H*p_pred(1)/(H*H*p_pred(1)+ R);
p_est(1) = p_pred(1)*(1-H*k(1));
x_est(1) = x_pred(1)+k(1)*residual(1);

%Calculate the rest of the values.
for j=2:nlen
    %Calculate the state and the output
    x(j) = A*x(j-1) + w(j);
    z(j) = H*x(j)+ v(j);
    %Predictor equations
    x_pred(j) = A*x_est(j-1);
    residual(j) = z(j)- H*x_pred(j);
    p_pred(j) = A*A*p_est(j-1)+ Q;
    %Corrector equations
    k(j) = H*p_pred(j)/(H*H*p_pred(j)+ R);
    p_est(j) = p_pred(j)*(1-H*k(j));
    x_est(j) = x_pred(j) + k(j)*residual(j);
    x_pred(j) = x_est(j);
    
    end
disp(k);
j=1:nlen;

%% plot Kalman gain for 10 iterations
plot(k(j));
title('Steady state value of gain')
xlabel('Number of iterations');
ylabel('Gain');

Matlab code for computer exercise 7.3 (d)

clc;
clear;
close all;
%Define the length of the simulation.
nlen = 10;

%Define the system.  
a = 0.8;  % a=1 for a constant,
h = 2;

%Define the noise covariances.
Q = 1;
R = 0.64;
%% Generate AR(3) process
phi =0;
e = randn(100,1); 
b = [1 phi]; 
y = filter(b,[-0.1 -0.09 0.648],e);
% initial estimates (guesses) for
% 1) the state.
x_0 = y(1);

% 2) the  error covariance.
p_0 = 0.8;

%Calculate the process and measurement noise.
w1 = randn(1,nlen); %This line and the next can be commented out after running
v1 = randn(1,nlen); %the script once to generate multiple runs with identical
%noise (for better comparison).
w = w1*sqrt(Q);
v = v1*sqrt(R);

%Calculate the first estimates for all values based upon the initial guess
%of the state and the a posteriori covariance. The rest of the steps will
%be calculated in a loop.
%
%Calculate the state and the output
x(1) = a*x_0 + w(1);
z(1) = h*x(1) + v(1);
%Predictor equations
x_pred(1) = a*x_0;
residual(1) = z(1)- h*x_pred(1);
p_pred(1) = a*a*p_0 + Q;
%Corrector equations
k(1) = h*p_pred(1)/(h*h*p_pred(1)+ R);
p_est(1) = p_pred(1)*(1-h*k(1));
x_est(1) = x_pred(1)+ k(1)*residual(1);

%Calculate the rest of the values.
for j=2:nlen
    %Calculate the state and the output
    x(j) = a*x(j-1) + w(j);
    z(j) = h*x(j) + v(j);
    %Predictor equations
    x_pred(j) = a*x_est(j-1);
    p_pred(j) = a*a*p_est(j-1) + Q;
    residual(j) = z(j)-h*x_pred(j);
    
    %Corrector equations
    k(j) = h*p_pred(j)/(h*h*p_pred(j)+ R);
    p_est(j) = p_pred(j)*(1-h*k(j));
    x_est(j) = x_pred(j)+k(j)*residual(j);
    x_pred(j) = x_est(j);
    
    end
disp(k);
 j=1:nlen;

%% plot Kalman gain for 10 iterations
figure(1);
plot(k(j));
legend('kalman gain');
title('Steady state value of gain')
xlabel('Number of iterations');
ylabel('Gain');
figure(2);
plot(x(j)); hold on;
plot(x_est(j)); hold on;
%plot(z(j)); hold on;
title('True Measured and estimated value')
xlabel('Number of iterations');
ylabel('values');
legend('True','Estimated');

Matlab code for computer exercise 7.3 (e)

clc;
clear;
close all;
%Define the length of the simulation.
nlen = 11;

%Define the system.  Change these to change the system.
A = -0.95;  
H = 1;

%Define the noise covariances.
Q = 1; % Variance of process noise
R = 0.8; % Variance of measurement noise

%% Generate AR(3) process
phi = 0;
e1 = randn(100,1); 
b = [1 phi]; 
y = filter(b,[-0.95 -0.9025],e1);

% initial estimates (guesses) for
% 1) the state.
x_0 = y(1);

% 2) the  error covariance.
p_0 = 0.6;

%Calculate the process and measurement noise.
w1 = randn(1,nlen); %This line and the next can be commented out after running
v1 = randn(1,nlen); %the script once to generate multiple runs with identical
%noise (for better comparison).
w = w1*sqrt(Q);
v = v1*sqrt(R);

%Calculate the first estimates for all values based upon the initial guess
%of the state and the a posteriori covariance. The rest of the steps will
%be calculated in a loop.
%
%Calculate the state and the output
x(1) = A*x_0 + w(1);
z(1) = H*x(1) + v(1);
%Predictor equations
x_pred(1) = A*x_0;
p_pred(1) = A*A*p_0 + Q;
residual(1) = z(1)- H*x_pred(1);

%Corrector equations
k(1) = H*p_pred(1)/(H*H*p_pred(1)+ R);
p_est(1) = p_pred(1)*(1-H*k(1));
x_est(1) = x_pred(1)+k(1)*residual(1);

%Calculate the rest of the values.
for j= 2:nlen
    %Calculate the state and the output
    x(j) = A*x(j-1)+ w(j) - w(j-1);
    z(j) = H*x(j) + v(j);
    %Predictor equations
    x_pred(j) = A*x_est(j-1);
    residual(j) = z(j)- H*x_pred(j);
    p_pred(j) = A*A*p_est(j-1)+ Q;
    %Corrector equations
    k(j) = H*p_pred(j)/(H*H*p_pred(j)+ R);
    p_est(j) = p_pred(j)*(1 - H*k(j));
    x_est(j) = x_pred(j)+ k(j)*residual(j);
    
    
end
disp(k);

%% plot Kalman gain for 10 iterations
j= 1:nlen;
figure(1);
plot(k(j));
legend('kalman gain');
title('Steady state value of gain')
xlabel('Number of iterations');
ylabel('Gain');
figure(2);
plot(x(j)); hold on;
plot(x_est(j)); hold on;
%plot(z(j)); hold on;
title('True Measured and estimated value')
xlabel('Number of iterations');
ylabel('values');
legend('True','Estimated');

Matlab Code Simulation Output: