Discrete-time stochastic processes, Approximate a signal by Pade\'s method and P

Question

Discrete-time stochastic processes, Approximate a signal by Pade's method and Prony's method ( in MATLAB).

This is a programming work. You need write a m-file which includes:

1. Create a zero-mean and unit variance white noise w[n], n=1:100.

2. Filter w[n] to obtained a ARMA(3,3) signal x[n]. Design the filter with 3 zeros and 3 poles:
Use the first 3 digits of your UIN as 3 zeros.
Use the inverse of the last 3 digits that is bigger than 1 as the poles.
[example, if UIN is 381936315, then your system has zeros 3, 8, and 1, has poles 1/5, 1/3, and 1/6]

3. Apply Pade's method to approximate x[n] using ARMA(10,10) model.
[i.e, find H(z)=B(z)/A(z) so that h[n] approximate x[n],
where both B(z) and A(z) have order 10]

4. Apply Prony's method to approximate x[n] using ARMA(10,10) model

5. Compare the results in step (3) and step (4)

Explanation / Answer

Given a transfer function G, the script computes an approximate reduced order model of the system using the technique of moment matching.
If
R=Pade_Approximation(G,r)

If the power series expansion of G is
G=c0+c1*s+c2*s^2+...+c_(2r)s^(2r)+c_(2r+1)s^(2r+1)+...
and then that of R would be
R=c0+c1*s+c2*s^2+...+c_(2r)s^(2r)+.....
with the first 2r coefficients of both series being matched.

matlab program
function PadeApprox=Pade_Approximation(G,r)
%Function PadeApprox=Pade_Approximation(G,r)
%
% Computes the r-th order Pade Approximation of a given n-th order
% transfer function G, with 1<=r<=n.
% Example
% G=tf([1 2],[1 3 4 5])
% r=2;
% R=Pade_Approximation(G,r)
%
%gives the output
%Transfer function:
% 0.04762 s + 0.8571
%----------------------
% s^2 + 0.7619 s + 2.143
%
% S. Janardhanan
% janas@ee.iitd.ac.in
%
% Code last updated on 18-Sep-2008

error(nargchk(2,2,nargin));
error(nargoutchk(1,1,nargout));
if ~isa(G,'tf') || ~isscalar(G)
error('Input needs to be a SISO transfer function');
end

if ~isreal(r) || (fix(r)~=r) || (r<1) || (r>n)
error('Invalid value of reduced model order')
end
[num,den]=tfdata(G,'v');
D_fact=num(1)/den(1);
num=num-D_fact*den;
num1=num(end:-1:1)/den(1);
den1=den(end:-1:1)/den(1);
n=length(den1)-1;
%For Pade Approximation
%Finding c_i's
c(1)=num1(1)/den1(1);
for i=2:min(n,2*r)
c(i)=(num1(i)-sum(den1(2:(i)).*c(end:-1:1)))/den1(1);
end
if (2*r)>n
for i=n+1:(2*r)
c(i)=-sum(den1(2:n+1).*c(end:-1:(end-n+1)))/den1(1);
end
end

%Finding Coefficients
C1=c(repmat((r+1:-1:2),r,1)+repmat((0:(r-1))',1,r));
b=-inv(C1)*c(1:r)';
a(r)=0;
for i=1:r
a(i)=c(i:-1:1)*b(1:i);
end
b=[1 b(end:-1:1)'];
a=a(end:-1:1);

[a,b]=tfdata(tf(a,b)+D_fact,'v');
b=b.*(abs(b)>1e-6);
a=a.*(abs(a)>1e-6);
%Pade Approximant

PadeApprox=tf(a,b);

matlab code for Prony's method.

function [iapp,ai,a_list,tau_list,omega_list,SUB_IND,energy,p]=applyprony(test_time,test_data,N,SUB_N,criteria_val)
% Implements MATLAB Signal Processing Tooltbox prony function
% performprony GUI

test_t_increment=test_time(2)-test_time(1);
fs=1/test_t_increment; % Sampling frequency=5 e8
[inumz,idenz]=prony(test_data,N,N);
% Compute and plot the impulse response of the Prony approximating system
[iapp,tapp]=impz(inumz,idenz,length(test_time),fs);
[r,p,k]=residuez(inumz,idenz);
a_list=r(:);
spoles=log(p(:))*fs;
tau_list=1./real(spoles);
omega_list=imag(spoles);
n_size=size(test_time,1);
n=0:1:n_size;
for ct=1:N
energy(ct)=abs(a_list(ct)^2)*(sum(abs(p(ct).^n).^2));
end

if (criteria_val==1)
[Mag, ISort]=sort(abs(a_list));
else
[En,ISort]=sort(energy);
end
ISort=ISort(end:-1:1);
FULL_IND=[ISort(:)]' ;
ai=zeros(size(test_time));
SUB_IND=FULL_IND(1:SUB_N);
for i=SUB_IND
ai=ai+a_list(i)*exp(spoles(i)*tapp);
end
ai=real(ai);
% Remove bug, PA should be same if modes are equal to model order
if(SUB_N==N)
ai(1)=iapp(1);
end
% Update the SUB_N for prony fit when no. of modes is not specified
if(SUB_N==0)
SUB_IND=FULL_IND;
end

Related Questions

Navigate

• Browse (All)
• Browse D
• Subjects

Integrity-first tutoring: explanations and feedback only — we do not complete graded work. Learn more.