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.

Loading…

Taking too long?

Reload document

| Open in new tab

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
[a,b] = pade(i,p,q);
imp = [1; zeros(100,1)];
h1 = filter(b,a,imp);   %Pade imp res
PadeFilter = fvtool(b,a); title('Filter Design Using Pade Method');

[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');

Matlab function for Pade Approximation (pade.m)

function [a,b] = pade(x,p,q)
%Pade Approximation
%   
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: