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.

Loading…

Taking too long?

Reload document

| Open in new tab

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)

[ad,bd,errd] = prony(h,p,q);
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));
[ad2,bd2,errd2] = prony(y,p,q);
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: