# Introduction to the z transform by ghkgkyyt

VIEWS: 5 PAGES: 20

• pg 1
```									                             Introduction to the Z-transform
University of Penn
November 9, 2007

Abstract
The following is an article from the University of Penn’s website, that introduces the Z-transform to
students. The original document can be viewed at http://www.ling.upenn.edu/courses/ling525/
z.html. The following, is an attempt to convert it to L TEX, interpreting the ASCII-art formulae into
A
The content is the same – with some changes in notation – namely the usage of j rather than i
as the complex “imaginary” constant (since i is used for current in Electrical Engineering), and the
usage of the asterisk (*) rather than the hash (#) to indicate convolution – consistent with Wikipedia.
This document was typeset by Stuart Longland, and is provided as-is. The author of this document
is unknown.
The L TEX source used to generate this document is available at http://dev.gentoo.org/~redhatter/
A
misc/z-transform.7z1 .

1
7-zip archive: see http://www.7-zip.org

1
1     Introducing the Z-transform
After reading this section, you may want to look at Chapter 1, “Signal Processing Basics,” in the User’s
Guide for the Matlab Signal Processing Toolbox. There are a few things that we have not covered, and
will not cover – for instance, the state space and the partial fraction expansion ways of representing
linear systems. Some other topics, such as lattice form representations, will be covered later. However,
you should now be able to follow much of the exposition, and beneﬁt from seeing the range of available
signal processing functions.
Note that the User’s Guide is a rather large (about 6 MB) .pdf ﬁle, so you will probably prefer to read
it via a fast network connection. Although the whole Guide is 720 pages long, the cited chapter is only
about 45 pages, so it is plausible to print it out.

1.1    Background
In this segment, we will be dealing with the properties of sequences made up of integer powers of some
complex number:

x[n] = z n ∀n ∈ (−∞..∞) ; z = some complex number
You should start with a clear graphical intuition about what such sequences are like. If the number z
happens to be one or zero, we will get a sequence of constant values. If z is a positive real number, we
will get a sampled exponential ramp, that is either rising or falling depending on whether z is less than 1
or greater than 1:

n = -50:50;   z = .97;   x = z.^n;           z = 1.03;   x = z.^n;
plot(n,x,’go’); title(’z = .97’);            plot(n,x,’go’); title(’z = 1.03’);

2
If the number z is a negative real number, we will get a sequence that alternates between positive
and negative values. Depending on whether z is less than, equal to, or greater than -1, these values will
increase exponentially, remain constant, or decrease exponentially in magnitude with increasing n.

z = -1.03; x = z.^n;                         z = -1; x = z.^n;
plot(n,x,’go’); title(’z = -1.03’);          plot(n,x,’go’); title(’z = -1’);

z = -.97; x = z.^n;
plot(n,x,’go’); title(’z = -.97’);

3
If z is a complex number that happens to lie on the unit circle (in the complex plane, with the x-axis
representing the real part of z and the y-axis representing the imaginary part), then z n will be a sinusoid
sampled at intervals of z radians.
This can easily be shown by considering that in this case z = ejω for some real number ω, which is
equivalent to cos(ω) + j sin(ω), so that z n will be ejnω , or cos(nω) + j sin(nω).

z = exp(i*pi/10); x = z.^n;                   z = exp(i*pi/20); x = z.^n;
plot3(n, real(x), imag(x), ’o’);              plot3(n, real(x), imag(x), ’o’);
hold on;                                      hold on;
plot3(n, real(x), imag(x), ’-’);              plot3(n, real(x), imag(x), ’-’);
title(’z = exp(i*pi/10)’);                    title(’z = exp(i*pi/20)’);
hold off; view(-15,15);                       hold off; view(-30,30);

If z is a complex number (with a non-zero imaginary part) whose magnitude is slightly less than
or greater than 1, then z n will be a sinusoidal spiral like those we have just seen, whose magnitude is
exponentially decreasing or increasing depending on whether z has a magnitude less than or greater than
1.

z = .97*exp(i*pi/10); x = z.^n;               z = 1.03*exp(i*pi/10); x = z.^n;
plot3(n, real(x), imag(x), ’o’);              plot3(n, real(x), imag(x), ’o’);
hold on;                                      hold on;
plot3(n, real(x), imag(x), ’-’);              plot3(n, real(x), imag(x), ’-’);
title(’z = .97*exp(i*pi/10)’);                title(’z = 1.03*exp(i*pi/10)’);
hold off; view(-15,15);                       hold off; view(-15,15);

4
Now suppose a linear time-invariant (LTI) system with impulse response h(n) has as its input one of
these exponential sequences

x[n] = z n
for some (arbitrarily chosen) complex number z.
The output will be the convolution sum h ∗ x, or

y[n] =           h[k]x[n − k]
k
=          h[k]z n−k
k
= z    n
h[k]z −k
k

Thus the input was z n , and the output is z n multiplied by a constant that depends on the value of z
and on the impulse response h. If we write that constant as H(z) = k h[k]z −k then we can rewrite the
system equation as (switching back to MATLAB)

y[n] = H(z) ∗ z n
We can see another way to say this if we express the LTI system with impulse response h in the more
general form of the “convolution matrix” Mh – that is, a matrix with a set of shifted time-reversed copies
of the impulse response h in its rows, as discussed in an earlier lecture. Now since

y[n] = Mh z n
= H(z)z n

we can see that the complex exponential z n is an eigenvector of Mh , with H(z) as the eigenvalue, since
Mh times z n equals the constant H(z) times z n .
Because of the superposition property of linear systems, this “eigenrelationship” makes it convenient
to express a signal as a linear combination of complex exponentials. If

n
x[n] =               ak zk
k

then the system output for an LTI system with impulse response h will be

n
y[n] =         ak H(zk )zk
k

That is, the output is also a linear combination of the same set of complex exponential sequences, with
each coeﬃcient being the product of the input coeﬃcient ak and the system’s eigenvalue H(zk ) for the
n
eigenfunction zk .
The expression for the system eigenvalues in terms of z

H(z) =                h[n]z −n
n

is known as the “Z-transform” of h (for n from −∞ to ∞). In nicer notation:
∞
H(z) =                  h[n]z −n
n=−∞

This equation is closely related to that for the DFT. Recall that the DFT is
2π
X(k) =          x[n]e−j N kn
n

5
If we rewrite the exponential on the right hand side slightly as
2π    −n
ejk N
each value of k can be seen as just picking a diﬀerent complex number
jk2π
z=e           N

to serve as the basis for a complex exponential series.
However, the correct analogy is with the DTFT, not the DFT.
Remember that the DFT relates x[n], a periodic function of a discrete variable x in the time domain,
to X[k], a periodic function of a discrete variable k in the frequency domain.
The “discrete time fourier transform” (DTFT) relates x[n], a nonperiodic function of a discrete variable
x in the time domain, to X(ω), a periodic function of a continuous variable ω in the frequency domain.
As a result, the DTFT shares with the Z-transform the fact that the transform equation deﬁnes a
function of a continuous (complex) variable, and that the input side of the equation sums over all n rather
than a ﬁnite set:

∞
DTFT: X(ω) =                      x[n]e−jωn
n=−∞
∞
=               x[n](ejω )−n
n=−∞

Thus the DTFT is exactly the Z-transform for z = ejω .
Since ejx = cos(x)+j sin(x), restricting z to imaginary powers of e is the same as requiring that |z| = 1,
i.e. that (in the complex plane) z must fall on a circle with radius one.
We can look at this another way. Let’s express the complex number z in polar form as rejω . Then the
Z-transform of x[n]

X(z) =             x[n]z −n
n

becomes

X(rejω ) =            x[n](rejω )−n
n
or

X(rejω ) =           x[n]r−n e−jωn
n

The right-hand side of this equation is just the DTFT of the sequence x[n] multiplied by a real
exponential r−n . Thus

X(rejω ) = X(z) = DT F T (x[n]r−n )
where r is the magnitude of z. Where r = 1, X(z) = DT F T (x).

6
1.2    History of the Z-transform
The z transform is unusual, in being named after a letter of the alphabet rather than a famous math-
ematician. The Fourier transform is named after Baron Jean Baptiste Joseph Fourier (1768-1830); the
Walsh-Hadamard transform is named after J.L. Walsh (?) and Jacques Salomon Hadamard (1865-1963);
we haven’t discussed the Laplace and Hilbert transforms yet, but we will (at least brieﬂy), and they are
named after Pierre-Simon de Laplace (1749-1827) and David Hilbert (1862-1943) respectively.
Laplace transforms have long been used in solving (continuous-time) linear constant-coeﬃcient diﬀer-
ential equations.
According to p 420 of Contemporary Linear Systems (Strum and Kirk 1994),
A method for solving linear, constant-coeﬃcient diﬀerence equations by Laplace transforms was
introduced to graduate engineering students by Gardner and Barnes in the early 1940s. They
applied their procedure, which was based on jump functions, to ladder networks, transmission
lines, and applications involving Bessel functions. This approach is quite complicated and in a
separate attempt to simplify matters, a transform of a sampled signal or sequence was deﬁned
in 1947 by W. Hurewicz as

∞
Z [f (kT )] =         f (kT )z −k
k=0

which was later denoted in 1952 as a “Z-transform” by a sampled-data control group at
Columbia University led by professor John R. Raggazini and including L.A. Zadeh, E.I. jury,
R.E. Kalman, J.E. Bertram, B. Friedland, and G.F. Franklin.
The Hurewicz equation is not expressed in the same way as the Z-transform we have introduced – it
is one-sided, and it is expressed as a function of the sampled data sequence f rather than the complex
number z – but the relationship is clear, and the applications were similar from the beginning. So perhaps
the Z-transform should really be called the “HurewicZ-transform” – but it is too late to change.
In any case, it is presumably not an accident that the Z-transform was invented at about the same
time as digital computers.

1.3    Z-transform of some Simple Signals
Consider x[n] = an u[n]. Its Z-transform is
∞                      ∞
X(z) =            an u[n]z −n =         (az −1 )n
n=−∞                   n=0

This will converge (the sum will be ﬁnite) just in case |z| > |a|.
Since for any geometric series
∞
c
crn =
n=0        1−r
the Z-transform sum above will be
1         z
X(z) =             =     ROC:|z| > |a|
1 − az −1   z−a
This equation for the Z-transform of x[n] = an u[n]
z
X(z) =
z−a
is a “rational function”, that is, a ratio of polynomials. We can characterize it by its zeros (the roots
of the numerator) and its poles (the roots of the denominator).
In this case there is one zero (z = 0) and one pole (z = a).
We also need to know the “region of convergence” (ROC) for the Z-transform (here |z| > |a|). For the
next section of exposition, we will neglect the ROC. This is not a good idea in general.

7
1.4     Properties of the Z-transform
It follows directly from the deﬁnition
X(z) =        x[n]z −n
n

that the Z-transform is LINEAR: that is, if

x[n] = cx1[n]

then

X(z) =            cx1[n]z −n
n

= c         x1[n]z −n
n
= cX1(z)

and if

x[n] = x1[n] + x2[n]
then

X(z) =             (x1[n] + x2[n])z −n
n
=         x1[n]z −n +         x2[n]z −n
n                   n
= X1(z) + X2(z)

Thus if
x[n] = an u[n] + bn u[n]
then

z         z
X(z) =           +
(z − a) (z − b)
z(z − b) + z(z − a)
=
(z − a)(z − b)
z(2z − (b + a))
=
(z − a)(z − b)

In addition to linearity, Z-transforms have a number of other properties that make them a useful tool
in analyzing LTI systems:

1.4.1     The convolution property of Z-transforms
As with the special case of the Fourier transform, convolution in the time domain is multiplication in the
z domain:
If

x[n] = x1[n] ∗ x2[n]
then

X(z) = X1(z)X2(z)

8
1.4.2       The shift property of Z-transforms
If

x[n] = x1[n − k]
then

X(z) = z −k X1(z)

1.4.3       The time reversal property of Z-transforms
If

x[n] = x1[−n]
then

X(z) = X1(1/z)

1.5     Z-transform of arbitrary Linear Constant-Coeﬃcient Diﬀerence Equa-
tions
The convolution property of the Z-transform means that for an arbitrary LTI system with impulse response
h[n], the system equation

y[n] = h[n] ∗ x[n]
implies that

Y (z) = H(z)X(z)
and therefore

H(z) = Y (z)/X(z)
We can also express the same LTI system as a linear constant-coeﬃcient diﬀerence equation, if the
system is causal:

a0 y[n] + a1 y[n − 1] + · · · + aM y[n − M ] = b0 x[n] + b1 x[n − 1] + · · · + bN x[n − N ]
If this concept is not entirely clear to you, you may want to review the lecture notes on Digital Filters
as Linear Constant-Coeﬃcient Diﬀerence Equations.
We can take the Z-transform of both sides of the above equation, representing the Z-transform by the
notation Z {· · ·}.

Z {a0 y[n] + a1 y[n − 1] + · · · + aM y[n − M ]} = Z {b0 x[n] + b1 x[n − 1] + · · · + bN x[n − N ]}

Since the Z-transform is linear, we can take it term-wise:

Z {a0 y[n]} + Z {a1 y[n − 1]} + · · · + Z {aM y[n − M ]} = Z {b0 x[n]} + Z {b1 x[n − 1]} + · · · + Z {bN x[n − N ]}

Likewise we can pull the b and a constants outside each transformed term:

a0 Z {y[n]} + a1 Z {y[n − 1]} + · · · + aM Z {y[n − M ]} = b0 Z {x[n]} + b1 Z {x[n − 1]} + · · · + bN Z {x[n − N ]}

9
Now the shift property of the Z-transform lets us replace (for all k) Z {x[n − k]} with z −k Z {x[n]} and
Z {y[n − k]} with z −k Z {y[n]}:

a0 Z {y[n]} + a1 z −1 Z {y[n]} + · · · + aM z −M Z {y[n]} = b0 Z {x[n]} + b1 z −1 Z {x[n]} + · · · + bN z −N Z {x[n]}

Now all the Z {· · ·} expressions are either Z {x[n]} Z {[y]} and so we can replace them all with X(z)
or Y (z):

a0 Y (z) + a1 z −1 Y (z) + · · · + aM z −M Y (z) = b0 X(z) + b1 z −1 X(z) + · · · + bN z −N X(z)
Factoring out Y (z) and X(z):

Y (z) a0 + a1 z −1 + · · · + aM z −M = X(z) b0 + b1 z −1 + · · · + bN z −N
The convolution property of the Z-transform told us that H(z), the Z-transform of the system’s impulse
response, is equal to Y (z)/X(z), so let’s solve for Y (z)/X(z) in our equation:

b0 + b1 z −1 + · · · + bN z −N
Y (z)/X(z) =
a0 + a1 z −1 + · · · + aM z −M
By convention, a0 is 1. We can multiply through by 1/b0 , replacing the b coeﬃcients with c coeﬃcients
such that cn = bn /b0 :

1 1 + c1 z −1 + · · · + cN z −N
Y (z)/X(z) =
b0 1 + a1 z −1 + · · · + aM z −M
If a0 is already 1, then the c coeﬃcients are just the same as the b coeﬃcients. Thus the Z-transform of
the impulse response of such a system — ANY system described by a linear constant-coeﬃcient diﬀerence
equation — is a ratio of polynomials in z −1 , where the coeﬃcients in the numerator come from the x
(input) coeﬃcients in the diﬀerence equation, and the coeﬃcients in the denominator come from the y
(output) coeﬃcients in the diﬀerence equation.

1.6    Z-transform poles and zeros as a Characterization of a System
Let’s go back to our expression for H(z):

1 1 + c1 z −1 + · · · + cN z −N
Y (z)/X(z) =
b0 1 + a1 z −1 + · · · + aM z −M

The fundamental theorem of algebra tells us that we can express the numerator and denominator
polynomials uniquely as a product of ﬁrst order terms in z −1 = 1/z, with (complex) coeﬃcients q1 , ..., qN
and p1 , ..., pM :

1 (1 − q1 /z)(1 − q2 /z) · · · (1 − qN /z)
H(z) =
b0 (1 − p1 /z)(1 − p2 /z) · · · (1 − pM /z)
In this equation, each of the qn and pn is a complex number. If z equals one of the qn , the numerator
becomes zero, and so does the whole function H(z). These values are called “zeros” of the system. As z
approaches one of the pn , the denominator approaches zero, and the value of the function H(z) approaches
inﬁnity. These values are called “poles” of the system; graphically, you can imagine putting the complex
plane on the ground, so that a pole sticks up into the sky at the point where z = pn .
Starting with our original linear constant-coeﬃcient diﬀerence equation

a0 y[n] + a1 y[n − 1] + · · · + aM y[n − M ] = b0 x[n] + b1 x[n − 1] + · · · + bN x[n − N ]

we have given a mechanical procedure for deriving an expression for the (factored form of the) Z-
transform of the system impulse resonse h[n]:

10
1 (1 − q1 /z)(1 − q2 /z) · · · (1 − qN /z)
H(z) =
b0 (1 − p1 /z)(1 − p2 /z) · · · (1 − pM /z)
If all of the a coeﬃcients are zero, then the output is just a moving weighted average of the input, and
the denominator of our expression for H(z) will be 1. This is a ﬁlter with only zeros and no poles, an
“all-zero” ﬁlter.
Conversely, if all of the b coeﬃcients are 0 (except of course for b0 which must be 1, or the system
would ignore the input!), then the current output sample is basically predicted as a weighted combination
of the earlier output samples, and the numerator of our expression for H(z) will be 1. This is a ﬁlter with
only poles and no zeros, an “all-pole” ﬁlter.
If the diﬀerence equation has both a and b coeﬃcients, then the ﬁlter has both poles and zeros.
Some important properties of poles and zeros follow from this algebraic background. For example,
if we are dealing with real-valued signals – as we normally are – than the coeﬃcients of the LCCDE
must obviously also be real. It follows that the roots of the numerator and denominator polynomials in
the Z-transform will either be real, or will come in complex-conjugate pairs. Thus the zeros and poles
will likewise either be real, or will come in complex-conjugate pairs, since they are just the roots of the
numerator and denominator polynomials.

1.7    From LCCDE coeﬃcients to poles and zeros and back
This part is easy – given a computer program to do the calculations – because it is just a matter of factoring
or expanding polynomials. In Matlab, the function roots() will take us from polynomial coeﬃcients to
the roots of the polynomial, while poly() will take us from the roots back to the polynomial coeﬃcients.
The a or b coeﬃcients (such as are used in filter()) are polynomial coeﬃcients.
Thus if p is a vector containing the z-values of the poles of a system (the roots of the denominator
polynomial in its Z-transform), then the a coeﬃcients will be a = poly(p), while we can go the other way
with p = roots(a). Likewise for the zeros (the roots of the numerator polynomial).
See the next section for a concrete example.

1.8    From the poles and zeros of the Z-transform to the spectrum
Since we know that the Z-transform reduces to the DTFT for z = ejω , and we know how to calculate the
Z-transform of any causal LTI (i.e. the Z-transform of its impulse response) from the coeﬃcients of the
diﬀerence equation, we can write down an expression for its spectrum (i.e. the ratio of the spectrum of
the output to the spectrum of the input) by simply setting z = ejω in the Z-transform.
This expression will be a periodic, continuous function of the variable ω, which is frequency in radians.
The period is obviously 2π. We can get N evenly-spaced samples of one period of the spectrum (which is
n                                    n
also what the DFT gives us) by setting (in Matlab-ish) ω to N 2π for n = 0 to N − 1, i.e. z = ej N 2π for
n = 0 to N − 1.

11
This is quite easy to calculate. It is also especially easy to think about graphically, if we represent the
system in pole-zero form:

>> circle = exp(i*(0:63)*2*pi/64); # 64 points around the unit circle
>> plot(real(circle),imag(circle),’o’);
>> r = .95*exp(i*.5*pi); # complex root at frequency .5, amplitude .95
>> axis([-2 2 -2 2]); axis(’equal’); hold on
>> plot(real(r),imag(r),’bx’, real(r’),imag(r’),’bx’); # plot complex conjugate
poles
>> plot(0,0,’bo’); # plot zero

12
Now we just evaluate the factored form of the Z-transform of a system with these poles (and the single,
necessary, degenerate zero). This is the product of the diﬀerences between the z-values (here the points
on the unit circle) and the roots r and r . The sampled spectrum is just the inverse of this product of
distances:

>>   distances = (abs(circle-r) .* abs(circle-r’));
>>   plot(distances);
>>   plot(1:64,log(1./distances),’bx:’);
>>   axis([0 32 -1 4]);

Of course, an alternative method of calculation – generally simpler in practice – is simply to take the
DFT of the impulse response of the ﬁlter.

>> A = poly([r r’])
A =
1.0000   -0.0000    0.9025
>> impulse = zeros(256,1); impulse(15)=1;
>> impresp = filter([1 0 0],A,impulse);
>> plot(impresp);
>> lspecplot(fft(impresp),1);

13
1.9    Relationship between pole amplitude and damping
1.10     Stability of recursive ﬁlters
[TK]

1.11     From a recorded signal to estimated AR ﬁlter coeﬃcients: LPC analysis
We’ve mentioned that the recursive part of the LCCDE ﬁlter – the output side, which contributes the a
coeﬃcients – is sometimes called an autoregressive ﬁlter. It expresses each sample as a linear combination
of Na previous samples, and we can think of this as being like predicting each sample in the sequence by
linear regression on the previous Na samples.
Suppose we are given a signal that was actually generated by such an autoregressive ﬁlter – one in
which all the b coeﬃcients except for b0 are 0. Will each sample actually be exactly the inner product of
the ﬁlter’s a coeﬃcients with the Na previous samples? No, not exactly, because the input to the system
will still have an eﬀect, being fed in through the b0 coeﬃcient. However, for many combinations of inputs
and system values, we will still be able to estimate the a coeﬃcients in a useful way. One ubiquitous
method of speech analysis – linear predictive coding , or LPC – is based on this assumption.
Assume that within some (fairly short) region of speech, M samples long, we are predicting each sample
as a linear combination of N previous samples, i.e.

p(i) =            a(j) ∗ s(i − j), 1 ≤ j ≤ N, N ≤ i < M
j

In matrix notation,
p = Sa
where p is an M × 1 matrix (column vector) of predicted values, S is an M × N matrix whose ith
row contains speech samples s(i) to s(i + N − 1), and a is an N × 1 matrix (column vector) of prediction
weights.
Thus for a third-order predictor, with input samples s1 · · · sM and predictor coeﬃcients [a(1)a(2)a(3)],
p = Sa means
                                                     
p(4)                  s(1)       s(2)       s(3)

   p(5)    


           s(2)       s(3)       s(4)    
         

   p(6)



           s(3)       s(4)       s(5)

     a(1)
=                                          ×  a(2) 
                                                              

   p(7)                s(4)       s(5)       s(6)    
     .
.
              .
.          .
.          .
.
     a(3)
.                     .          .          .
                                                     
                                                     
p(M )               s(M − 3) s(M − 2) s(M − 1)
Notice that this is closely related to the previously-discussed concept of a recursive ﬁlter, and exempliﬁes
why such ﬁlters are commonly called “autoregressive”.
In this case, the a(i) are just the coeﬃcients of the left-hand (output) side of a linear constant-coeﬃcient
diﬀerence equation — ignoring the a(0) coeﬃcient in of an LDCE, which is 1 by convention.
We know S in p = Sa, but not a or p. What should we do? Suppose that we want to choose the
weights a so as to minimize the prediction error

norm(Sa − s)
This is what is called a “least squares” problem: ﬁnd a vector x providing the “best” solution to an
overdetermined system of equations Ax = b, where “best” means minimizing norm(Ax − b).
There are a variety of methods of solving such problems, depending on the properties of A and b.
MATLAB oﬀers us two solutions, x = A \ b and x = pinv(A)*b. We’ve seen the “backslash” A \ b
solution before; we’ll discuss the pseudoinverse solution pinv(A)*b at greater length later.
For now, we’ll just observe that if A has more rows than columns and is not of full rank, then “choose
x to minimize norm(Ax − b)” does not have a unique solution. The solution x = pinv(A)*b give us the

14
smallest x (x minimizing norm(x)) while the solution x = A \ b gives the x with the fewest possible
nonzero components.
We’ll examine the this question as we come to it in practice, but if there is a diﬀerence, the pseudo-
inverse solution is probably the one that we want.
Let’s get a chunk of speech to work with:

>> S1 = S(2246:2347);
>> plot(S1);
>> length(S1)
ans =
102

This is just a convenient piece of the middle of a vowel from the audio ﬁle we have used before. It
represents as close to 3 pitch periods as we can get in this sampled signal:

Now let’s get things set up as speciﬁed.
If we had a vector of input samples [1 2 3 4 5 6 7 8 9 10], and we were going to construct a third-order
predictor, and we wanted to avoid making hypotheses about samples outside the range of input samples
we would need Ax = b to come out as
                                   
1   2   3                    4
2   3   4                    5
                                   
                                   
                                 

   3   4   5   
     x1     
   6    

4   5   6    ×  x2  =    7
                                 
                                     
                                   

   5   6   7   
     x3     
   8    


   6   7   8   


   9    

7   8   9                    10
In this situation, A is pretty close to what is called a “Hankel” matrix (see help hankel in Matlab),
and so we can get A set up by

function A = makeA(S,N)
%
% Set up matrix for linear prediction calculation
% S is a vector of length M whose ith sample will be predicted
% as a linear combination of the N previous samples.
% Thus we want to set up a matrix A,
% with M-N rows and N columns,
% whose jth row consists of the samples from S(j) to S(j+N-1).
% Obviously M must be greater than N.
AA = hankel(S,S(1:N));
A = AA(1:(length(S)-N),:);

15
Thus we get

>> makeA(1:10,3)
ans =
1     2       3
2     3       4
3     4       5
4     5       6
5     6       7
6     7       8
7     8       9

The role of b in our equation Ax = b will be played by the input samples from S(N + 1) to S(M ).

>> A = makeA(S1, 14);
>> b = S1(15:length(S1));

A is full rank

>> rank(A)
ans =
14

and so it doesn’t matter much what method we use to get x:

>> norm(A\b - pinv(A)*b)
ans =
5.0799e-15

>> x = pinv(A)*b;
>> x’
ans =
Columns 1 through 7
-0.1660    0.3545   -0.4917        0.5188     -0.7772       1.0481   -1.3567
Columns 8 through 14
1.1888   -1.5503    1.4317       -1.0071      0.8159    -1.2254      1.4033

How much of the variance of the input have we accounted for?

>> sum((A*x-mean(A*x)).^2) / sum((b-mean(b)).^2)
ans =
0.9294

16
Let’s look at the (windowed version of the) original signal, and its DFT spectrum:

>>   S2 = hamming(length(S1)).*S1;
>>   plot(S2);
>>   q = zeros(512,1);
>>   q(1:length(S2)) = S2;
>>   lspecplot(fft(q), 8000);

17
Now let’s look at the impulse response of the recursive ﬁlter whose coeﬃcients are x — well, OK, whose
coeﬃcients are 1 followed by -1 times x backwards... we have to obey the conventions derived from the
LCDE formulation:

>> impulse = zeros(256,1); impulse(15)=1;
>> fcoef = [1 -x(14:-1:1)’]
fcoef =
Columns 1 through 7
1.0000   -1.4033    1.2254   -0.8159    1.0071   -1.4317    1.5503
Columns 8 through 14
-1.1888    1.3567   -1.0481    0.7772   -0.5188    0.4917   -0.3545
Column 15
0.1660
>> impresp = filter([1 0 0 0 0 0 0 0 0 0 0 0 0 0 0], fcoef, impulse);
>> plot(impresp);

and the (log) amplitude spectrum of this impulse response:

>> lspecplot(fft(impresp),8000);

18
Now let’s ﬁgure out what the pole frequencies and amplitudes are:

>> angle(roots(fcoef))*4000/pi
ans =
1.0e+03 *
3.2923
-3.2923
2.9825
-2.9825
2.3412
-2.3412
1.4361
-1.4361
1.4294
-1.4294
0.4731
-0.4731
0.6518
-0.6518
>> abs(roots(fcoef))
ans =
0.9653
0.9653
0.8436
0.8436
0.8628
0.8628
0.9703
0.9703
0.8048
0.8048
0.9803
0.9803
0.7576
0.7576

In tabular form, the seven complex poles, corresponding to the seven complex-conjugate roots of the
predictor polynomial, have frequencies and amplitudes of
Frequency (Hz) Amplitude (0-1)
473.1        .9803
651.8        .7576
1,429.4        .8048
1,436.1        .9703
2,341.2        .8628
2,982.5        .8436
3,292.3        .9653
Note that by performing a 14-th order analysis, we’ve guaranteed that we will get seven complex poles
(well, we might have substituted a couple of real poles for one of them). If we want to use this method to
ﬁnd the vowel formants, we have to decide which poles correspond to formats and which do not. Here, it
is clear (because of the amplitude) that F 1 is 473 Hz., and F 2 is 1,436 Hz. (though of course the apparent
estimation accuracy is spurious). It is less clear what F 3 should be.

19
1.12     Inverse ﬁltering
Suppose we have a causal LTI system S1, with impulse response h, whose z transform, H(z), is a ratio of
polynomials Bpoly /Apoly .
Now suppose we have another causal LTI system S2, with impulse response g, whose Z-transform,
G(z), happens to be exactly the inverse, namely Apoly /Bpoly .
If we convolve some input x with h, and then convolve the result with g – x ∗ h ∗ g – this is equivalent
to convolving x with the convolution of h with g – x ∗ (h ∗ g).
By the convolution rule for the Z-transform, the Z-transform of (h ∗ g) will be the product H(z)G(z).
Given the way we constructed the two systems, this product will be 1 for all values. Therefore the combined
system (h ∗ g) will do nothing to its input.
More interestingly, if we already have the output S1(x) = x ∗ h, we can apply the second system to
the result to get the original input back: S2(S1(x)) = x. This process is called inverse ﬁltering .
Let’s try this with a simple one-pole recursive ﬁlter with a center frequency half the sampling rate,
excited by an impulse:

r = .95*exp(i*.5*pi);
A= poly([r r’]);
impulse=zeros(256,1); impulse(15)=1;
out1 = filter([1 0 0],A,impulse);
out1(1:30)’
ans =
Columns 1 through 7
0         0         0         0         0                       0           0
Columns 8 through 14
0         0         0         0         0                       0           0
Columns 15 through 21
1.0000    0.0000   -0.9025   -0.0000    0.8145                   0.0000   -0.7351
Columns 22 through 28
-0.0000    0.6634    0.0000   -0.5987   -0.0000                 0.5404      0.0000
Columns 29 through 30
-0.4877   -0.0000
out2 = filter(A,[1 0 0],out1);
out2(1:30)’
ans =
Columns 1 through 12
0     0     0     0     0     0     0     0                 0       0      0        0
Columns 13 through 24
0     0     1     0     0     0     0     0                 0       0      0        0
Columns 25 through 30
0     0     0     0     0     0

Success! we got the impulse back
Now let’s try the recusive ﬁlter estimated by the LPC modeling in the previous section. Here we have
seven poles – because we used a 14th-order model. If we apply this ﬁlter to a known input – say an impulse
– and then inverse ﬁlter, we’ll get our known input back:
[TK]
But it’s more interesting to inverse ﬁlter the original speech. here we don’t actually know what the
excitation was, since we are starting from a real-world signal. In fact, the input wasn’t really generated by
this type of ﬁlter at all, but rather by a physical process that is something like – can be usefully modeled
as – such a ﬁlter. So the result of inverse ﬁltering in this case will be a hypothetical signal – what the
input would have been if the output really were created by our modeled ﬁlter.
[TK]

20

```
To top