4
Frequency Analysis
Frequency analysis of any given signal involves the transformation of a time-domain
signal into its frequency components. The need for describing a signal in the frequency
domain exists because signal processing is generally accomplished using systems that are
described in terms of frequency response. Converting the time-domain signals and
systems into the frequency domain is extremely helpful in understanding the character-
istics of both signals and systems.
In Section 4.1,the Fourier series and Fourier transform will be introduced. The
Fourier series is an effective technique for handling periodic functions. It provides a
method for expressing a periodic function as the linear combination of sinusoidal
functions. The Fourier transform is needed to develop the concept of frequency-domain
signal processing. Section 4.2 introduces the z-transform,its important properties,and
its inverse transform. Section 4.3 shows the analysis and implementation of digital
systems using the z-transform. Basic concepts of discrete Fourier transforms will be
introduced in Section 4.4,but detailed treatments will be presented in Chapter 7. The
application of frequency analysis techniques using MATLAB to design notch filters and
analyze room acoustics will be presented in Section 4.5. Finally,real-time experiments
using the TMS320C55x will be presented in Section 4.6.
4.1 Fourier Series and Transform
In this section,we will introduce the representation of analog periodic signals using
Fourier series. We will then expand the analysis to the Fourier transform representation
of broad classes of finite energy signals.
4.1.1 Fourier Series
Any periodic signal, x(t),can be represented as the sum of an infinite number of
harmonically related sinusoids and complex exponentials. The basic mathematical
representation of periodic signal x(t) with period T
0
(in seconds) is the Fourier series
defined as
Real-Time Digital Signal Processing. Sen M Kuo,Bob H Lee

1
T
0

T
0
xte
ÀjkV
0
t
dt: 4:1:2
This integral can be evaluated over any interval of length T
0
. For an odd function,it is
easier to integrate from 0 to T
0
. For an even function,integration from ÀT
0
=2toT
0
=2
is commonly used. The term with k  0 is referred to as the DC component because
c
0

1
T
0

T

T
0
2
À
T
0
2
Ae
ÀjkV
0
t
dt 
A
T
0
e
ÀjkV
0
t
ÀjkV
0




t
2
À
t
2

0
,the effect of decreasing t is to
spread the signal power over the frequency range. On the other hand,when t is fixed but
the period T
0
increases,the spacing between adjacent spectral lines decreases.
t
x(t)
A
−
t
2
t
2
0
2
2
−T
0
−T
0
T
0
T
0
Figure 4.1 Rectangular pulse train
128
FREQUENCY ANALYSIS
A periodic signal has infinite energy and finite power,which is defined by Parseval's
theorem as

2
represents the power of the kth harmonic component of the signal,the total
power of the periodic signal is simply the sum of the powers of all harmonics.
The complex-valued Fourier coefficients, c
k
,can be expressed as
c
k
 c
k
jje
jf
k
: 4:1:6
A plot of jc
k
j versus the frequency index k is called the amplitude (magnitude) spectrum,
and a plot of f
k
versus k is called the phase spectrum. If the periodic signal x(t) is real
valued,it is easy to show that c
0
is real valued and that c
k
and c
Àk
are complex
conjugates. That is,
c
k

0
.
Example 4.2: Consider the output of an ideal oscillator as the perfect sinewave
expressed as
xtsin 2pf
0
t, f
0

V
0
2p
:
We can then calculate the Fourier series coefficients using Euler's formula
(Appendix A.3) as
sin2pf
0
t
1
2j
e
j2pf
0
t
À e
Àj2pf
0
t
ÀÁ


 2p=T
0
. The
number of frequency components increases as T
0
is increased,whereas the envelope of
the magnitude of the spectral components remains the same. If we increase the period
without limit (i.e., T
0
3I),the line spacing tends toward 0. The discrete frequency
components converge into a continuum of frequency components whose magnitudes
have the same shape as the envelope of the discrete spectra. In other words,when the
period T
0
approaches infinity,the pulse train shown in Figure 4.1 reduces to a single
pulse,which is no longer periodic. Thus the signal becomes non-periodic and its
spectrum becomes continuous.
In real applications,most signals such as speech signals are not periodic. Consider the
signal that is not periodic (V
0
3 0orT
0
3I),the number of exponential components
in (4.1.1) tends toward infinity and the summation becomes integration over the entire
continuous range (ÀI,I. Thus (4.1.1) can be rewritten as
xt
1
2p

I

I
ÀI
e
Àat
ute
ÀjVt
dt


I
0
e
ÀajVt
dt

1
a  jV
:
The Fourier transform XV is also called the spectrum of the analog signal x(t). The
spectrum XV is a complex-valued function of frequency V,and can be expressed as
XV


XV


e
jfV
, 4:1:12
where jXVj is the magnitude spectrum of x(t),and fV is the phase spectrum of x(t).

k
can be expressed in terms of XV using (4.1.2) and (4.1.10) as
c
k

1
T
0
XkV
0
: 4:1:14
For a given finite interval function,its Fourier transform at a set of equally spaced
points on the V-axis is specified exactly by the Fourier series coefficients. The distance
between adjacent points on the V-axis is 2p=T
0
radians.
If x(t) is a real-valued signal,we can show from (4.1.9) and (4.1.10) that
FT xÀt  X
Ã
V and XÀVX
Ã
V:4:1:15
It follows that
jXÀVj  jXVj and fÀVÀfV: 4:1:16
Therefore the amplitude spectrum jXVj is an even function of V,and the phase
spectrum is an odd function.
If the time signal x(t) is a delta function dt,its Fourier transform can be calculated as
XV

I

sinV
0
t jpdV  V
0
ÀdV À V
0

cosV
0
t pdV  V
0
dV À V
0

sgnt
1, t ! 0
À1, t < 0
&
2
jV
Table 4.2 Useful properties of the Fourier transform
Time function xt Property Fourier transform XV
a
1
x
1
ta
2
x
2

xt sinV
0
t Modulation
1
2j
XV À V
0
ÀXV  V
0

xt cosV
0
t Modulation
1
2
XV  V
0
XV À V
0

e
Àat
xt Frequency shifting XV  a
132
FREQUENCY ANALYSIS
Example 4.4: Find the Fourier transform of the time function
yte
Àajtj
, a > 0:
This equation can be written as

series
Xz

I
nÀI
xnz
Àn
, 4:2:1
where Xz represents the z-transform of xn. The variable z is a complex variable,and
can be expressed in polar form as
z  re
jy
, 4:2:2
where r is the magnitude (radius) of z,and y is the angle of z. When r  1, jzj1is
called the unit circle on the z-plane. Since the z-transform involves an infinite power
series,it exists only for those values of z where the power series defined in (4.2.1)
THE Z-TRANSFORM
133
converges. The region on the complex z-plane in which the power series converges is
called the region of convergence (ROC).
As discussed in Section 3.1,the signal xn encountered in most practical applications
is causal. For this type of signal,the two-sided z-transform defined in (4.2.1) becomes a
one-sided z-transform expressed as
Xz

I
n0
xnz
Àn
: 4:2:3

j < 1:
The equivalent condition for convergence (or ROC) is
jzj > jaj:
Thus we obtain Xz as
Xz
z
z À a
, jzj > jaj:
There is a zero at the origin z  0 and a pole at z  a. The ROC and the pole±zero
plot are illustrated in Figure 4.2 for 0 < a < 1,where `Â' marks the position of the
pole and `o' denotes the position of the zero. The ROC is the region outside
the circle with radius a. Therefore the ROC is always bounded by a circle since the
convergence condition is on the magnitude jzj. A causal signal is characterized by
an ROC that is outside the maximum pole circle and does not contain any pole.
The properties of the z-transform are extremely useful for the analysis of discrete-time
LTI systems. These properties are summarized as follows:
1. Linearity (superposition). The z-transform is a linear transformation. Therefore the
z-transform of the sum of two sequences is the sum of the z-transforms of the
individual sequences. That is,
134
FREQUENCY ANALYSIS
|z| = a
|z| = 1
Re[z]
Im[z]
Figure 4.2 Pole,zero,and ROC (shaded area) on the z-plane
ZTa
1
x
1

2
z are the z-transforms of the
signals x
1
n and x
2
n,respectively. This linearity property can be generalized for
an arbitrary number of signals.
2. Time shifting. The z-transform of the shifted (delayed) signal ynxn À k is
YzZTxn À k  z
Àk
Xz, 4:2:5
2. where the minus sign corresponds to a delay of k samples. This delay property states
that the effect of delaying a signal by k samples is equivalent to multiplying its
z-transform by a factor of z
Àk
. For example,ZTxn À 1  z
À1
Xz. Thus the unit
delay z
À1
in the z-domain corresponds to a time shift of one sampling period in the
time domain.
3. Convolution. Consider the signal
xnx
1
nÃx
2
n, 4:2:6
2. where à denotes the linear convolution introduced in Chapter 3,we have

Àan
cz
z À e
Àa
sin!
0
n
z sin!
0

z
2
À 2z cos!
0
1
cos!
0
n
zz À cos!
0

z
2
À 2z cos!
0
1
4.2.2 Inverse z-transform
The inverse z-transform can be expressed as
xnZT
À1

m0
a
m
z
Àm
, 4:2:9
where Az and Bz are expressed in either descending powers of z,or ascending powers
of z
À1
. Dividing Bz by Az obtains a series of negative powers of z if a positive-time
sequence is indicated by the ROC. If a negative-time function is indicated,we express
Xz as a series of positive powers of z. The method will not work for a sequence defined
136
FREQUENCY ANALYSIS
in both positive and negative time. In addition,it is difficult to obtain a closed-form
solution of the time-domain signal xn via the long-division method.
The long-division method can be performed recursively. That is,
xn b
n
À

n
m1
xn À ma
m
45D
a
0
, n  1,2, FFF 4:2:10
where

2
À x1a
1
À x0a
2
=a
0
 3:6439,
FFF
This yields the time domain signal xnf1,3,3:6439, FFFg obtained from long
division.
The partial-fraction-expansion method factors the denominator of Xz if it is
not already in factored form,then expands Xz into a sum of simple partial fractions.
The inverse z-transform of each partial fraction is obtained from the z-transform
tables such as Table 4.3,and then added to give the overall inverse z-transform. In
many practical cases,the z-transform is given as a ratio of polynomials in z or z
À1
as
shown in (4.2.9). If the poles of Xz are of first order and M  L À 1,then Xz can be
expanded as
Xzc
0


LÀ1
l1
c
l
1 À p
l

l





zp
l
: 4:2:13
THE Z-TRANSFORM
137
If the order of the numerator B(z) is less than that of the denominator A(z) in (4.2.9),
that is L À 1 < M,then c
0
 0. If L À 1 > M,then X(z) must be reduced first in order to
make L À 1 M by long division with the numerator and denominator polynomials
written in descending power of z
À1
.
Example 4.7: For the z-transform
Xz
z
À1
1 À 0:25z
À1
À 0:375z
À2
we can first express X(z) in positive powers of z,expressed as
Xz
z

1
z  0:5




z0:75
 0:8
and
c
2

1
z À 0:75




zÀ0:5
À0:8:
Thus we have
Xz
0:8z
z À 0:75
À
0:8z
z  0:5
:
The overall inverse z-transform x(n) is the sum of the two inverse z-transforms.
From entry 3 of Table 4.3,we obtain


j
for an mth order pole at z  p
l
. The
coefficients g
j
may be obtained with
138
FREQUENCY ANALYSIS
g
j

1
m À j3
d
mÀj
dz
mÀj
z À p
l

m
Xz
z
!
zp
l
: 4:2:14
Example 4.8: Consider the function


d
dz
z  1




z1
 1,
g
2

z À 1
2
Xz
z





z1
 z  1j
z1
 2:
Thus
Xz
z
z À 1

Thus the inversion integral in (4.2.8) can be easily evaluated using Cauchy's residue
theorem expressed as
xn
1
2pj

c
Xzz
nÀ1
dz


residues of Xzz
nÀ1
at poles of Xzz
nÀ1
within C:
4:2:16
THE Z-TRANSFORM
139
The residue of Xzz
nÀ1
at a given pole at z  p
l
can be calculated using the formula
R
zp
l

d

Example 4.9: Given the following z-transform function:
Xz
1
z À 1z À 0:5
,
we have
Xzz
nÀ1

z
nÀ1
z À 1z À 0:5
:
This function has a simple pole at z  0 when n  0,and no pole at z  0 for
n ! 1. For the case n  0,
Xzz
nÀ1

1
zz À 1z À 0:5
:
The residue theorem gives
xnR
z0
 R
z1
 R
z0:5
 zXzz
nÀ1

z0:5
 2 À 20:5
nÀ1
 21À0:5
nÀ1
hi
, n ! 1:
We have discussed three methods for obtaining the inverse z-transform. A limitation
of the long-division method is that it does not lead to a closed-form solution. However,
it is simple and lends itself to software implementation. Because of its recursive nature,
care should be taken to minimize possible accumulation of numerical errors when the
number of data points in the inverse z-transform is large. Both the partial-fraction-
expansion and the residue methods lead to closed-form solutions. The main disadvan-
tage with both methods is the need to factor the denominator polynomial,which is done
by finding the poles of X(z). If the order of X(z) is high,finding the poles of X(z) may be
140
FREQUENCY ANALYSIS
a difficult task. Both methods may also involve high-order differentiation if X(z)
contains multiple-order poles. The partial-fraction-expansion method is useful in gen-
erating the coefficients of parallel structures for digital filters. Another application of z-
transforms and inverse z-transforms is to solve linear difference equations with constant
coefficients.
4.3 Systems Concepts
As mentioned earlier,the z-transform is a powerful tool in analyzing digital systems. In
this section,we introduce several techniques for describing and characterizing digital
systems.
4.3.1 Transfer Functions
Consider the discrete-time LTI system illustrated in Figure 3.8. The system output is
computed by the convolution sum defined as ynxnÃhn. Using the convolution
property and letting ZTxn  Xz and ZT yn  Yz,we have

systems,as illustrated in Figure 4.4. In the cascade (series) interconnection,the output
of the first system, y
1
n,is the input of the second system,and the output of the second
system, y(n),is the overall system output. From Figure 4.4(a),we have
Y
1
zXzH
1
z and YzY
1
zH
2
z:
Thus
YzXzH
1
zH
2
z:
Therefore the overall transfer function of the cascade of the two systems is
Hz
Yz
Xz
 H
1
zH
2
z: 4:3:3
Since multiplication is commutative, H

2
z: 4:3:6
x(n)
x(n)
H
1
(z)
H
1
(z)
H(z)
H(z)
H
2
(z)
H
2
(z)
y
1
(n)
X(z)
y(n)
y(n)
Y(z)
Y
1
(z)
y
1

ÀÁ
 H
1
zH
2
z:
Thus the overall system H(z) can be realized as the cascade of the first-order
system H
1
z1 À z
À1
and the second-order system H
2
z1 À z
À1
À z
À2
.
4.3.2 Digital Filters
The general I/O difference equation of an FIR filter is given in (3.1.16). Taking the
z-transform of both sides,we have
Yzb
0
Xzb
1
z
À1
XzÁÁÁb
LÀ1
z

LÀ1
l0
b
l
z
À1
: 4:3:8
The signal-flow diagram of the FIR filter is shown in Figure 3.6. FIR filters can be
implemented using the I/O difference equation given in (3.1.16),the transfer function
defined in (4.3.8),or the signal-flow diagram illustrated in Figure 3.6.
Similarly,taking the z-transform of both sides of the IIR filter defined in (3.2.18)
yields
Yzb
0
Xzb
1
z
À1
XzÁÁÁ b
LÀ1
z
ÀL1
XzÀa
1
z
À1
YzÀÁÁÁ Àa
M
z
ÀM

b
l
z
Àl
1 

M
m1
a
m
z
Àm

Bz
1  Az
, 4:3:10
where Bz

LÀ1
l0
b
l
z
Àl
and Az

M
m1
a
m

− a
1
− a
2
− a
M
y(n−M)
y(n−2)
y(n−1)
z
−1
z
−1
b
1
b
2
b
L−1
x(n−L+1)
z
−1
z
−1
Figure 4.6 Detailed signal-flow diagram of IIR filter
144
FREQUENCY ANALYSIS
Hz
b
0

m1
z À p
m


b
0
z À z
1
z À z
2
ÁÁÁz À z
M

z À p
1
z À p
2
ÁÁÁz À p
M

: 4:3:12
The roots of the numerator polynomial are called the zeros of the transfer function H(z).
In other words,the zeros of H(z) are the values of z for which Hz0,i.e.,Bz0.
Thus H(z) given in (4.3.12) has M zeros at z  z
1
, z
2
, FFF, z
M

Hz
Yz
Xz

1
L

LÀ1
l0
z
Àl

1
L
1 À z
ÀL
1 À z
À1
!
: 4:3:13
This equation can be rearranged as
Yzz
À1
Yz
1
L
XzÀz
ÀL
Xz
ÂÃ

zplane(b, a),which displays the pole±zero diagram of H(z).
Example 4.12: Consider the IIR filter with the transfer function
Hz
1
1 À z
À1
 0:9z
À2
:
We can plot the pole±zero diagram using the following MATLAB script:
b  [1];a [1, À1, 0.9];
zplane(b, a);
Similarly,we can plot Figure 4.7 using the following MATLAB script:
b  [1, 0, 0, 0, 0, 0, 0, 0, 1];a [1, À1];
zplane(b, a);
As shown in Figure 4.7,the system has a single pole at z  1,which is at the same
location as one of the eight zeros. This pole is canceled by the zero at z  1. In this case,
the pole±zero cancellation occurs in the system transfer function itself. Since the system
Re[z]
Im[z]
|z|=1
Zero
Pole
Figure 4.7 Pole±zero diagram of the moving-averaging filter, L  8
146
FREQUENCY ANALYSIS
output YzXzHz,the pole±zero cancelation may occur in the product of system
transfer function H(z) with the z-transform of the input signal Xz. By proper selection
of the zeros of the system transfer function,it is possible to suppress one or more poles of
the input signal from the output of the system,or vice versa. When the zero is located

results in 1 À a L À 1=L,which is slightly less than 1. When L is large,i.e.,a longer
window,the pole is closer to the unit circle.
x(n)
(
•
)
2
w(n) = x
2
(n)
H(z)
y(n) = P
x
(n)
ˆ
Figure 4.8 Block diagram of recursive power estimator
Re[z]
Im[z]
|z| = 1
Zero
Pole
Figure 4.9 Pole±zero diagram of the recursive power estimator
SYSTEMS CONCEPTS
147
An LTI system Hz is stable if and only if all the poles are inside the unit circle. That is,
jp
m
j < 1 for all m: 4:3:16
In this case,lim
n3I

, n ! 0.
4.3.4 Frequency Responses
The frequency response of a digital system can be readily obtained from its transfer
function. If we set z  e
j!
in H(z),we have
Hz


z  e
j!


I
nÀI
hnz
Àn


z  e
j!


I
nÀI
hne
Àj!n
 H!: 4:3:17
Thus the frequency response of the system is obtained by evaluating the transfer
function on the unit circle jzjje

is a first-order FIR filter. Taking the z-transform of both sides and re-arranging
the terms,we obtain
Hz
1
2
1  z
À1
ÀÁ
:
From (4.3.17),we have
H!
1
2
1  e
Àj!
ÀÁ

1
2
1  cos ! À j sin !,
jH!j
2
fReH!g
2
fImH!g
2

1
2
1  cos !,

!
2
hi
À
!
2
:
As discussed earlier,MATLAB is an excellent tool for analyzing signals in the
frequency domain. For a given transfer function, H(z),expressed in a general form in
(4.3.10),the frequency response can be analyzed with the MATLAB function
[H, w ] freqz(b, a, N);
which returns the N-point frequency vector w and the N-point complex frequency
response vector H,given its numerator and denominator coefficients in vectors b and
a,respectively.
SYSTEMS CONCEPTS
149
Example 4.15: Consider the difference equation of IIR filter defined as
ynxnyn À 1À0:9yn À 2: 4:3:19a
This is equivalent to the IIR filter with the transfer function
Hz
1
1 À z
À1
 0:9z
À2
: 4:3:19b
The MATLAB script to analyze the magnitude and phase responses of this IIR
filter is listed (exam 4_15.m in the software package) as follows:
b  [1];a [1, À1, 0.9];
[H, w ] freqz(b, a, 128);

2
 b
1
z  b
2
z
2
 a
1
z  a
2
: 4:3:20
The roots of the characteristic equation
z
2
 a
1
z  a
2
 0 4:3:21
are the poles of the filter,which may be either real or complex. For complex poles,
p
1
 re
jy
and p
2
 re
Àjy
, 4:3:22

q
q
Im[z]
Re[z]
|z| = 1
Figure 4.10 A second-order IIR filter with complex-conjugated poles
Re[z]
Im[z]
z = e
jw
z = −1
p
1
p
2
V
2
V
1
U
1
U
2
z
1
z
2
Figure 4.11 Geometric evaluation of the magnitude response from the pole±zero diagram
Similarly,we can obtain two zeros,z
1

e
j!
À z
1
e
j!
À z
2

e
j!
À p
1
e
j!
À p
2

: 4:3:26
Assuming that b
0
 1,the magnitude response of the system can be shown as
jH!j 
U
1
U
2
V
1
V

151