Category: Electronics

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)-------------------------------

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:

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)

[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:

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
[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: