Chapter 4 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
1
=============================================================== % % 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];
2
[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 =========================================================
3
Chapter 5 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
4
