Categories

## Matlab Code for Computer Exercise C8.6 Monson H. Hayes Statistical DSP

In this post, we will discuss Matlab Code for Computer Exercise C8.6 Monson H. Hayes Statistical Digital Signal Processing book. You can also refer the question below. The output of matlab code simulation is available at the end of this page.

Matlab Code for computer exercise C8.6 is shown below.

``````%% Part (a)
clc;
clear;
close all;
w1=pi/2;
w2=1.1*pi/2;
phi1=0:0.1:(2*pi)+0.1;
phi2=0:0.1:(2*pi)+0.1;
n=0:1:length(phi1)-1;
b=[1 2 0 -2 -1];
w=filter(b,1,randn(64,1)); %data generation with gaussian noise - 100 samples
[Pxx,F] = pwelch(w);

plot(F,10*log10(Pxx));hold on
title('power spectrum of w(n)');
xlabel('Frequency');
ylabel('psd in dB');
grid on;
x=2.*cos(n.*w1+phi1)+2.*cos(n.*w2+phi2)+w';
[Pxx_x,F_x] = pwelch(x);
plot(F_x,10*log10(Pxx_x));
title('power spectrum of x(n)');
xlabel('Frequency');
ylabel('psd in dB');
legend('psd of w(n)','psd of x(n)');
grid on;
% end part (a)-------------------------------------------

%% Part (b)
clc;
clear;
close all;
w1=pi/2;
w2=1.1*pi/2;
phi1=0:0.1:(2*pi)+0.1;
phi2=0:0.1:(2*pi)+0.1;
n=0:1:length(phi1)-1;
b=[1 2 0 -2 -1];
w=filter(b,1,randn(64,1));
va = var(w);
disp('The variance is:'),disp(va);
[Pxx,F] = pwelch(w);
mag=abs(Pxx);
x1=2.*cos(n.*w1+phi1);
x2=2.*cos(n.*w2+phi2);
[Pxx_x1,F_x1] = pwelch(x1);
mag_x1=abs(Pxx_x1);
[Pxx_x2,F_x2] = pwelch(x2);
mag_x2=abs(Pxx_x2);
plot(F,10*log10(Pxx));hold on
plot(F_x1,10*log10(Pxx_x1));hold on
plot(F_x2,10*log10(Pxx_x2));hold on
title('power spectrum');
xlabel('Frequency');
ylabel('psd in dB');
legend('psd of w(n)','psd of x1(n)','psd of x2(n)');
grid on;
% end part (b) ------------------------------------------

%% Part (c)
clc;
clear;
close all;
w1=pi/2;
w2=1.1*pi/2;
phi1=0:0.1:(2*pi)+0.1;
phi2=0:0.1:(2*pi)+0.1;
n=0:1:length(phi1)-1;
b=[1 2 0 -2 -1];
w=filter(b,1,randn(64,1));

%Resolution
delw = 0.89*(2*pi/length(n));
disp('Resolution is:'),disp(delw);
disp('The values of x(n) to resolve frequencies in to two sinusoids:');
disp(length(n));

x1=2.*cos(n.*w1+phi1);
x2=2.*cos(n.*w2+phi2);
x=x1+x2+w';

pxx= periodogram(x);

pxx_x1 = periodogram(x1);
pxx_x2= periodogram(x2);

plot(10*log10(pxx));hold on
plot(10*log10(pxx_x1));hold on
plot(10*log10(pxx_x2));hold on

title('power spectrum');
xlabel('Frequency');
ylabel('psd in dB');
legend('psd of x(n)','psd of x1(n)','psd of x2(n)');
grid on;
%%% The number of values of x(n) to solve the frequencies of the two
%%% sinusoids are equal to length of phi.
%end part (c)--------------------------------------------------------

%% Part (d)
clc;
clear;
close all;
w1=pi/2;
w2=1.1*pi/2;
phi1=0:0.1:(2*pi)+0.1;
phi2=0:0.1:(2*pi)+0.1;
n=0:1:length(phi1)-1;
b=[1 2 0 -2 -1];
w=filter(b,1,randn(100,1));

x1=2.*cos(n.*w1+phi1);
x2=2.*cos(n.*w2+phi2);

pxx= periodogram(w);

pxx_x1 = periodogram(x1);
pxx_x2= periodogram(x2);

plot(10*log10(pxx));hold on
plot(10*log10(pxx_x1));hold on
plot(10*log10(pxx_x2));hold on

title('power spectrum');
xlabel('Frequency');
ylabel('psd in dB');
legend('psd of w(n)','psd of x1(n)','psd of x2(n)');
grid on;
%%% The number of values of x(n) to solve the frequencies of the two
%%% sinusoids are equal to length of phi.
%end part (d)---------------------------------------------------------

%% Part (e)
clc;
clear;
close all;
w1=pi/2;
w2=1.1*pi/2;
phi1=0:0.1:(2*pi)+0.1;
phi2=0:0.1:(2*pi)+0.1;
n=0:1:length(phi1)-1;
b=[1 2 0 -2 -1];

x1=2.*cos(n.*w1+phi1);
x2=2.*cos(n.*w2+phi2);
x_n1 =filter(x1,1,randn(256,1));
x_n2 =filter(x2,1,randn(256,1));

pxx_x1 = periodogram(x_n1);
pxx_x2= periodogram(x_n2);
figure(1)
subplot(2,1,1),plot(10*log10(pxx_x1));hold on
subplot(2,1,2),plot(10*log10(pxx_x2));hold on

title('power spectrum');
xlabel('Frequency');
ylabel('psd in dB');
%legend('psd of x1(n)','psd of x2(n)');
grid on;
%%% The two spectrums are visible
%% For 20 different realizations
xn1 =filter(x1,1,randn(256,20));
xn2 =filter(x2,1,randn(256,20));
figure(2)

pxx_xn1 = periodogram(xn1);
pxx_xn2= periodogram(xn2);

subplot(2,1,1);plot(10*log10(pxx_xn1));
title('power spectrum of x1(n) for 20 realizations');
xlabel('Frequency');
ylabel('psd in dB');
grid on;
subplot(2,1,2);plot(10*log10(pxx_xn2));

title('power spectrum of x2(n) for 20 realizations');
xlabel('Frequency');
ylabel('psd in dB');
grid on;
%end part (e)----------------------------------------

%% Part (f)
clc;
clear;
close all;
w1=pi/2;
w2=1.1*pi/2;
phi1=0:0.1:(2*pi)+0.1;
phi2=0:0.1:(2*pi)+0.1;
n=0:1:length(phi1)-1;
b=[1 2 0 -2 -1];
K = 8;
N = 256;
L = N/K;

x1=2.*cos(n.*w1+phi1);
x2=2.*cos(n.*w2+phi2);
x_n1 =filter(x1,1,randn(L,1));
x_n2 =filter(x2,1,randn(L,1));

pxx_x1 = bartlett(L);
pxx_x2= bartlett(L);
figure(1)
plot(10*log10(pxx_x1));hold on
plot(10*log10(pxx_x2));hold on

title('power spectrum');
xlabel('Frequency');
ylabel('psd in dB');
legend('psd of x1(n)','psd of x2(n)');
grid on;
%%% The two spectrums are visible
%% For 20 different realizations
xn1 =filter(x1,1,randn(L,20));
xn2 =filter(x2,1,randn(L,20));
%end part(f)-------------------------------
``````
Categories

## 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));
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');
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');
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:

Categories

## Matlab Code for Monson H. Hayes Statistical DSP Computer Exercise C4.2

Matlab Code for Monson H. Hayes Statistical DSP Computer Exercise C4.2 (Chapter 4 computer exercise) and its corresponding output images are shown below. The question (C4.2) is shown below.

Matlab Code:

``````clc;
clear;
close all;
%% Initialization for the problem

b = [1, -0.9, 0.81];
a = [1, -1.978, 2.853, -1.877, 0.9036];
H = tf(b,a,0.1,'Variable','z^-1');

%% Part (a)

imp = [1; zeros(100,1)];
h = filter(b,a,imp);
stem(0:100,h),
title('The first 100 samples of the unit sample response h(n)');

%% Part (b)

p = input('Please enter the number of poles, p: ');
q = input('Please enter the number of zeros, q: ');
n = input('Please enter the number of iteration, n: ');
[a,b,err] = ipf(h,p,q,n,1);
disp('--------------------------------------------')
disp('a coefficients using iterative prefiltering'),disp(a);
disp('b coefficients using iterative prefiltering'),disp(b);
disp('Error using iterative prefiltering'),disp(err);
disp('---------------------------------------------------');

%% Part (c)

variance = 0.0001;
v = sqrt(variance)*randn(size(h));  % white gaussian noise
SNR = snr(h,v);
y = awgn(h,SNR);
[ac,bc,errc] = ipf(y,p,q,n,1);
disp('a coefficients using iterative prefiltering for noisy measurements:');
disp(ac);
disp('b coefficients using iterative prefiltering for noisy measurements:')
disp(bc);
disp('Error using iterative prefiltering for noisy measurements:');
disp(errc);
disp('-----------------------------------------------------------');

%% Part (d)

disp('a coefficients using Prony method for exact measurements:'), disp(ad);
disp('b coefficients using Prony method for exact measurements:'), disp(bd);
disp('Error using Prony method for exact measurements::'), disp(min(errd));
disp('--------------------------------------------------------');
disp('a coefficients using Prony method for noisy measurements:'), disp(ad2);
disp('b coefficients using Prony method for noisy measurements:'), disp(bd2);
disp('Error using Prony method for Noisy measurements::'), disp(min(errd2));
``````

Matlab function for Prony Approximation (prony.m)

``````function [a,b,err] = prony(x,p,q)
%Prony
%
x = x(:);
N = length(x);
if p+q>=length(x), error('Model order too large'), end
X = convm(x,p+1);
Xq = X(q+1:N+p-1,1:p);
a = [1;-Xq\X(q+2:N+p,1)];
b = X(1:q+1,1:p+1)*a;
err = x(q+2:N)' *X(q+2:N,1:p+1)*a;
end
``````

Matlab function for iterative prefiltering (ipf.m)

``````function [a,b,err] = ipf(x,p,q,n,a)
%Iterative Prefiltering
%
x = x(:);
N = length(x);
if p+q>=length(x), error('Model order too large'), end
if nargin < 5
a = prony(x,p,q); end
delta = [1; zeros(N-1,1)];
for i=1:n
f = filter(1,a,x);
g = filter(1,a,delta);
u = convm(f,p+1);
v = convm(g,q+1);
ab = -[u(1:N,2:p+1) -v(1:N,:) ]\u(1:N,1);
a = [1; ab(1:p)];
b = ab(p+1 :p+q+1);
err = norm( u(1:N,1) + [u(1:N,2:p+1) -v(1:N, :)]*ab);
end``````

Output of Simulation for p = 4, q = 5, n = 2:

Categories

## Matlab Code for Monson H. Hayes Statistical DSP Computer Exercise C4.1

Matlab Code for Monson H. Hayes Statistical DSP Computer Exercise C4.1 (Chapter 4 problems) is shown below. The output of the matlab code is also available at the end of this page. Computer exercise C4.1 taken from the problem of chapter 4 of Monson H. Hayes Statistical Digital Signal Processing is shown below.

Matlab code for the above problem C4.1 is shown below:

``````% Matlab code for problem C4.1
clc;
%clear;
close all;
wc = input('Please enter the cutoff frequency, wc: ');
p = input('Please enter the number of poles, p: ');
q = input('Please enter the number of zeros, q: ');
nd = input('Please enter the delay of the filter, nd: ');
N = p + q;
for n = 0:99
if n == nd
i(n+1) = wc/pi;
else
i(n+1) = (sin((n-nd)*wc))/((n-nd)*pi);
end
end
imp = [1; zeros(100,1)];
h1 = filter(b,a,imp);   %Pade imp res

[a2,b2,err] = prony(i,p,q);
h2 = filter(b2,a2,imp);   %Prony imp res
for n=1:100
e1(n) = i(n)-h1(n); %error for PAde
e2(n) = i(n)-h2(n); %error for Prony
end
disp('---------------------------------------------');
disp('Minimum error for Prony:'), disp(err);

subplot(2,1,1), stem(e1), title('The Error for Pade Approximation');
subplot(2,1,2), stem(e2), title('The Error for Prony Approximation');
PronyFilter = fvtool(b2,a2); title('Filter Design Using Prony Method');

wn = pi/2;
[z,p,k] = ellip(10,0.5,20,wn/pi);
[sos] = zp2sos(z,p,k);      % Convert to SOS form
EllipticFilter = fvtool(sos); title('10th Order Eliptic Low Pass Filter');
``````

``````function [a,b] = pade(x,p,q)
%
x = x(:);
if p+q>=length(x), error('Model order too large'), end
X = convm(x,p+1);
Xq = X(q+2:q+p+1,2:p+1);
a = [1;-Xq\X(q+2:q+p+1,1)];
b = X(1:q+1,1:p+1)*a;
end
``````

Matlab function for Prony Approximation (prony.m)

``````function [a,b,err] = prony(x,p,q)
%Prony
%
x = x(:);
N = length(x);
if p+q>=length(x), error('Model order too large'), end
X = convm(x,p+1);
Xq = X(q+1:N+p-1,1:p);
a = [1;-Xq\X(q+2:N+p,1)];
b = X(1:q+1,1:p+1)*a;
err = x(q+2:N)' *X(q+2:N,1:p+1)*a;
end
``````

Output of simulation:

Categories

## 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. The problem, matlab code and output of simulation is shown below.

``````% 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');
``````
Categories

## Monson H. Hayes Statistical Digital Signal Processing Problem C3.3

Monson H. Hayes statistical Digital Signal Processing Problem C3.3 (Chapter 3 Computer Exercise): In this exercise we will look at how accurately the sample autocorrelation is in estimating the autocorrelation of white noise. The full problem is given below:

``````%MATLAB Code for the above problem
clc;
clear;
close all;
m = 0;
v = 1;
x = m + sqrt(v) .* randn(1000,1); % Generates 1000 samples of zero mean, m=0 and unit variance, v=1
rxb = autocorr(x,100);
[acf,lags] = autocorr(x); %returns the lag numbers that MATLAB uses to compute acf

for k = 1:100
for m = 1:9
for n = 1:99
rx(k) = (1/1000) * sum(x(n+100*m) * x(n-k+100*m));    %segment white noise sequence into 10 d/t sequences
end
end
end

xd = m + sqrt(v) .* randn(10000,1);
rxd = autocorr(xd,100);

subplot(4,1,1),stem(x);
title('1000 Samples of zero mean unit variance white Gaussian noise');
ylabel('x(n)');
xlabel('n');
subplot(4,1,2),autocorr(x,'NumLags',100);
subplot(4,1,3),autocorr(rx,'NumLags',99);
subplot(4,1,4),autocorr(xd,'NumLags',100);

``````