Z-transform and its applications in signal processing

% m-file to compute the first 5 values of the
% inverse z-transform using the power series method
% (Program 4D.1, p235; program name: prog4d1.m)
%
n = 5; % number of power series points
N1 = [1 -1.122346 1]; D1 = [1 -1.433509 0.85811];
N2 = [1 1.474597 1]; D2 = [1 -1.293601 0.556929];
N3 = [1 1 0]; D3 = [1 -0.612159 0];
B = [N1; N2; N3]; A = [D1; D2; D3];
[b,a] = sos2tf([B A]);
b = [b zeros(1,n-1)];
[h,r] = deconv(b,a); %perform long division
disp(h);
============================================================
%
% m-file for finding the partial fraction expansion of a z-transform
% (Program 4D.2, p237; program name: prog4d2.m)
%
N1=[1 -1.122346 1];
N2=[1 -0.437833 1];
N3=[1 1 0];
D1=[1 -1.433509 0.85811];
D2=[1 -1.293601 0.556929];
D3=[1 -0.612159 0];
sos=[N1 D1; N2 D2; N3 D3];
[b, a] = sos2tf(sos);
[r, p, k]=residue(b, a)

================================================
%
% A script illustrating cascade to parallel realization
% (Program 4B.3, p237; program name: prog4d3.m)
%
nstage=2;
N1 = [1 0.481199 1];
N2 = [1 1.474597 1];
D1 = [1 0.052921 0.83173];
D2 = [1 -0.304609 0.238865];
sos = [N1 D1; N2 N2];
[b, a] = sos2tf(sos);
[c, p, k] = residue(b, a);
m = length(b);
b0 = b(m)/a(m);
j=1;
for i=1:nstage
bk(j)=c(j)+c(j+1);
bk(j+1)=-(c(j)*p(j+1)+c(j+1)*p(j));
ak(j)=-(p(j)+p(j+1));
ak(j+1)=p(j)*p(j+1);
j=j+2;
end
b0
ak
bk
c
p
k

===============================================================
%
% A simple script illustrating basics of cascade to parallel
% conversion for a 4th order filter (cprealization.m)
%
b=[1, -2, 1, 0, 0];
a=[1, -0.41421, 0.08579, 0.292895, 0.5];
[c, p, k] = residuez(b, a);
c
p
k
b0 = b(3)/a(5)
b01=c(1)+c(2)
b11=-(c(1)*p(2)+c(2)*p(1))
a11=-(p(1)+p(2))
a12=p(1)*p(2)
b02=c(3)+c(4)
b12=-(c(3)*p(4)+c(4)*p(3))
a12=-(p(3)+p(4))
a22=p(3)*p(4)
=================================================================
%
% A simple script illustrating inverse z transform by power series codes
% File-name: pseries.m
b1 = [1 0.481199 1];
b2 = [1 1.474597 1];
a1 = [1 0.052921 0.83173];
a2 = [1 -0.304609 0.238865];
sos = [b1 a1; b2 a2];
[b, a] = sos2tf(sos);
m = length(b);
r = b;
for n=1:m
[h(n), r] = deconv(r, a);
h
r
r = [r(2:m) 0];
end

=================================================
function IZT
%program IZT (izt.m) is for:
%(1) computing the inverse z-transform via the power series or partial
% fraction expansion method
%(2) converting a transfer function, H(z), in cascade form to an equivalent
% transfer function in parallel, via partial fraction expansion. The
% basic building block is the second order biquad

IZT_Mode = 1; %1 - use power series to perform the IZT
%2 - use partial fraction/residue to perform the IZT
%3 - cascade to parallel conversion
n = 5; % number of power series points
b1 = [1 0.481199 1]; a1 = [1 0.052921 0.83173];
b2 = [1 1.474597 1]; a2 = [1 -0.304609 0.238865];
B = [b1; b2]; A = [a1; a2];
[b,a] = sos2tf([B A]);
if IZT_Mode==1 %compute the IZT by power series
b = [b zeros(1,n-1)];
[h,r] = deconv(b,a); %perform long division
disp('The result of the inverse Z-transform by power series is:'); disp(h);
else
[res,poles,rem] = residuez(b,a);
disp('The poles of the transform function are:'); disp(poles');
disp('The partial fraction coefficients are:'); disp(res');
if IZT_Mode==3
i = 1 ;
for j=1:size(B,1)
[b,a] = residuez(res(i:i+1),poles(i:i+1), 0);
fprintf('Numerator and Denominator for stage%d:\n',j);disp(b);disp(a);
i = i + 2;
end
end
end
=========================================================
%
% m-file to illustrate the computation of frequency response
% of an IIR filter using FFT (freqrespex1.m).
%
Fs=500; % sampling frequency
b1=[1 -1.6180 1]; b2=[]; % numerator/denominator filter coefficients
a1=[1 -1.5161 0.878]; a2=[];
B=[b1; b2]; A=[a1; a2];
[b,a]=sos2tf([B A]);
n=256; % number of frequency points
dx=0;
if dx
freqz(b,a,n,Fs); % Frequency response evaluation by FFT
else
f=(0:n-1)*(Fs/2)/n;
freqz(b,a,f,Fs); % frequency response by direct evaluation.
end
=========================================================
Correlation and Convolution

function Correlation
%Program to compute normalised/unnormalised auto- or crosscorrelation
crosscorrelation = 1 ;
normalized = 0 ;
if crosscorrelation == 0
x = [1 3 2 -1 -2]; %for autocorrelation
if normalized==0
R11 = xcorr(x,'biased') % unnormalized autocorrelation
else
R11 = xcorr(x,'coeff') % normalized autocorrelation
end
else
x1 = [4 2 -1 3 -2 -6 -5 4 5]; %for crosscorrelation
x2 = [-4 1 3 7 4 -2 -8 -2 1]; %for crosscorrelation
if normalized==0
R12 = xcorr(x1,x2,'biased') % unnormalized crosscorrelation
else
R12 = xcorr(x1,x2,'coeff') % normalized crosscorrelation
end
end

if crosscorrelation
subplot(3,1,1), plot(1:length(x1),x1), ylabel('x1(n)');
subplot(3,1,2), plot(1:length(x2),x2), ylabel('x2(n)');
subplot(3,1,3), plot(1:length(R12),R12), ylabel('R12');
else
subplot(2,1,1), plot(1:length(x),x), ylabel('x(n)');
subplot(2,1,2), plot(1:length(R11),R11), ylabel('R11');
end