Chapter 6

function Factors = factorize(polyn)

format long; Factors = [];

% Use threshold of 1e-8 instead of 0 to account for

% precision effects

THRESH = 1e-8;

%

proots = roots(polyn); % get the zeroes of the polynomial

len = length(proots); % get the number of zeroes

%

while(len > 1)

if(abs(imag(proots(1))) <>

fac = [1 -real(proots(1))];

% construct the factor with proots(1) as zero

Factors = [Factors;[fac 0]];

else % if the zero has imaginary part get all zeroes whose

% imag part is -ve of imaginary part of proots(1)

negimag = imag(proots)+imag(proots(1));

% get all zeroes which have same real part as proot(1)

samereal = real(proots)-real(proots(1));

%find the complex conjugate zero

index = find(abs(negimag)

if(index) % if the complex conjugate exists

fac = [1 -2*real(proots(1)) abs(proots(1))^2];

%form 2nd order factor

Factors = [Factors;fac];

else % if the complex conjugate does not exist

fac = [1 -proots(1)];

Factors = [Factors;[fac 0]];

end

end

polyn = deconv(polyn,fac);

%deconvolve the 1st/2nd order factor from polyn

proots = roots(polyn); %determine the new zeros

len = length(polyn); %determine the number of zeros

end

% Program 6_1

% Determination of the Factored Form

% of a Rational z-Transform

%

num = input('Type in the numerator coefficients = ');

den = input('Type in the denominator coefficients = ');

K = num(1)/den(1);

Numfactors = factorize(num)

Denfactors = factorize(den)

disp('Numerator factors');disp(Numfactors);

disp('Denominator factors');disp(Denfactors);

disp('Gain constant');disp(K);

zplane(num,den)

% Program 6_2

% Determination of the Rational z-Transform

% from its Poles and Zeros

%

format long

zr = input('Type in the zeros as a row vector = ');

pr = input('Type in the poles as a row vector = ');

% Transpose zero and pole row vectors

z = zr'; p = pr';

k = input('Type in the gain constant = ');

[num, den] = zp2tf(z, p, k);

disp('Numerator polynomial coefficients'); disp(num);

disp('Denominator polynomial coefficients'); disp(den);

%Program 6_3

% Partial-Fraction Expansion of Rational z-Transform

%

num = input('Type in numerator coefficients = ');

den = input('Type in denominator coefficients = ');

[r,p,k] = residuez(num,den);

disp('Residues');disp(r')

disp('Poles');disp(p')

disp('Constants');disp(k)

% Program 6_4

% Partial-Fraction Expansion to Rational z-Transform

%

r = input('Type in the residues = ');

p = input('Type in the poles = ');

k = input('Type in the constants = ');

[num, den] = residuez(r,p,k);

disp('Numerator polynomial coefficients'); disp(num)

disp('Denominator polynomial coefficients'); disp(den)

% Program 6_5

% Power Series Expansion of a Rational z-Transform

%

% Read in the number of inverse z-transform coefficients to be computed

L = input('Type in the length of output vector = ');

% Read in the numerator and denominator coefficients of

% the z-transform

num = input('Type in the numerator coefficients = ');

den = input('Type in the denominator coefficients = ');

% Compute the desired number of inverse transform coefficients

[y,t] = impz(num,den,L);

disp('Coefficients of the power series expansion');

disp(y')

% Program 6_6

% Power Series Expansion of a Rational z-Transform

%

% Read in the number of inverse z-transform coefficients

% to be computed

N = input('Type in the length of output vector = ');

% Read in the numerator and denominator coefficients of

% the z-transform

num = input('Type in the numerator coefficients = ');

den = input('Type in the denominator coefficients = ');

% Compute the desired number of inverse transform

% coefficients

x = [1 zeros(1, N-1)];

y = filter(num, den, x);

disp('Coefficients of the power series expansion');

disp(y)

%Program_6_7

% Stability Test

%

num = input('Type in the numerator vector = ');

den = input('Type in the denominator vector = ');

N = max(length(num),length(den));

x = 1; y0 = 0; S = 0;zi = zeros(1,N-1);

for n = 1:1000

[y,zf] = filter(num,den,x,zi);

if abs(y) <>

x = 0;

S = S + abs(y);

y0 = y;zi = zf;

end

if n <>

disp('Stable Transfer Function');

else

disp('Unstable Transfer Function');

end