DOS Matlab
DOS Matlab
SIGNAL
PROCESSING
with selected topics
ADAPTIVE SYSTEMS
TIME-FREQUENCY ANALYSIS
SPARSE SIGNAL PROCESSING
Ljubiša Stanković
2015
2
ISBN-13: 978-1514179987
ISBN-10: 1514179989
All right reserved. Printed and bounded in the United States of America.
To
my parents
Božo and Cana,
my wife Snežana,
and our
Irena, Isidora, and Nikola.
4
Contents
I Review 19
Chapter 1 Continuous-Time Signals and Systems 21
1.1 Continuous-Time Signals 21
1.2 Periodic Signals and Fourier Series 24
1.2.1 Fourier Series of Real-Valued Signals 28
1.2.2 Linear Systems 33
1.3 Fourier Transform 35
1.3.1 Fourier Transform and Linear Time-Invariant
Systems 37
1.3.2 Properties of the Fourier Transform 37
1.3.3 Relationship Between the Fourier Series and
the Fourier Transform 40
1.4 Fourier Transform and Stationary Phase Method 42
1.5 Laplace Transform 48
1.5.1 Linear Systems Described by Differential Equa-
tions 51
1.5.2 Table of the Laplace Transform 52
1.6 Butterworth Filter 53
5
6 Contents
Index 811
T
HIS
teaching and research in signal processing. It is written for students
and engineers as a first book in digital signal processing, assuming
that a reader is familiar with the basic mathematics, including integrals, dif-
ferential calculus, and linear algebra. Although a review of continuous-time
analysis is presented in the first chapter, a prerequisite for the presented
content is a basic knowledge about continuous-time signal processing.
The book consists of three parts. After an introductory review part, the
basic principles of digital signal processing are presented within Part two of
the book. This part starts with Chapter two which deals with basic defini-
tions, transforms, and properties of discrete-time signals. The sampling the-
orem, providing essential relation between continuous-time and discrete-
time signals, is presented in this chapter as well. Discrete Fourier transform
and its applications to signal processing are the topic of the third chapter.
Other common discrete transforms, like Cosine, Sine, Walsh-Hadamard,
and Haar are also presented in this chapter. The z-transform, as a power-
ful tool for analysis of discrete-time systems, is the topic of Chapter four.
Various methods for transforming a continuous-time system into a corre-
sponding discrete-time system are derived and illustrated in Chapter five.
Chapter six is dedicated to the forms of discrete-time system realizations.
Basic definitions and properties of random discrete-time signals are given
in Chapter six. Systems to process random discrete-time signals are consid-
ered in this chapter as well. Chapter six concludes with a short study of
quantization effects.
The presentation is supported by numerous illustrations and exam-
ples. Chapters within Part two are followed by a number of solved and
unsolved problems for practice. Theory is explained in a simple way with
a necessary mathematical rigor. The book provides simple examples and
13
14 Preface
London,
July 2013 - July 2015.
Author
Introduction
S
IGNAL
physical or symbolic representation of an information. Signal theory
and processing are the areas dealing with the efficient generation,
description, transformation, transmission, reception, and interpretation of
information. In the beginning, the most common physical processes used
for these purposes were the electric signals, for example, varying current
or electromagnetic waves. Signal theory is most commonly studied within
electrical engineering. Signal theory theory are strongly related to the ap-
plied mathematics and information theory. Examples of signals include
speech, music, image, video, medical, biological, geophysical, sonar, radar,
biomedical, car engine, financial, and molecular data. In terms of signal gen-
eration, the main topics are in sensing, acquisition, synthesis, and reproduc-
tion of information. Various mathematical transforms, representations, and
algorithms are used for describing signals. Signal transformations are a set
of methods for decomposition, filtering, estimation and detection. Modu-
lation, demodulation, detection, coding, and compression are the most im-
portant aspects of the signal transmission. In the process of interpretation,
various approaches may be used, including adaptive and learning-based
tools and analysis.
Mathematically, signals are presented by functions of one or more vari-
ables. Examples of one-dimensional signals are speech and music signals.
A typical example of a two-dimensional signal is an image while video se-
quence is a sample of a three-dimensional signal. Some signals, for example,
geophysical, medical, biological, radar, or sonar, may be represented and
interpreted as one-dimensional, two-dimensional, or multidimensional.
Signals may be continuous functions of independent variables, for ex-
ample, functions of time and/or space. Independent variables may also
be discrete, with the signal values being defined only over an ordered set
15
16 Introduction
xd(n)
x(n)
x(t)
8 1000
7 0111
0.4 0.4 6 0110
5 0101
4 0100
0.2 0.2 3 0011
2 0010
1 0001
0 0 0000
0 5 10 15 0 5 10 15 0 5 10 15
t n n
Figure 1 Illustration of a continuous signal and its discrete-time and digital version.
ANALOG SYSTEM
x(t) y(t)
ha(t)
DIGITAL SYSTEM
Figure 2 Illustration of an analog and a digital system used to process an analog signal.
Review
19
Chapter 1
Continuous-Time Signals and Systems
M
OST
time signals. In many applications, the result of signal process-
ing is presented and interpreted in the continuous-time domain.
Throughout the course of digital signal processing, the results will be dis-
cussed and related to the continuous-time forms of signals and their param-
eters. This is the reason why the first chapter is dedicated to a review of
signals and transforms in the continuous-time domain. This review will be
of help in establishing proper correspondence and notation for the presen-
tation that follows in the next chapters.
In the Heaviside function definition, the value of u(0) = 1/2 is also used.
Note that the independent variable t is continuous, while the signal itself is
not a continuous function. It has a discontinuity at t = 0.
The boxcar signal (rectangular window) is formed as b(t) = u(t +
1/2) − u(t − 1/2), that is, b(t) = 1 for −1/2 ≤ t < 1/2 and b(t) = 0 else-
where. A signal obtained by multiplying the unit-step signal by t is called
the ramp signal, with notation R(t) = tu(t).
21
22 Continuous-Time Signals and Systems
"∞ "t
u(t) = δ(τ )u(t − τ )dτ = δ(τ )dτ
−∞ −∞
or
du(t)
= δ ( t ). (1.4)
dt
A sinusoidal signal, with amplitude A, frequency Ω0 , and initial phase
ϕ, is a signal of the form
x ( t + T ) = x ( t ). (1.6)
is also periodic with period T = 2π/Ω0 . Fig. 1.1 depicts basic continuous-
time signals.
Ljubiša Stanković Digital Signal Processing 23
u(t) 1 1
δ(t)
0 0
-1 (a) -1 (b)
-4 -2 0 2 4 -4 -2 0 2 4
1 1
sin(πt)
b(t)
0 0
-1 (c) -1 (d)
-4 -2 0 2 4 -4 -2 0 2 4
t t
Figure 1.1 Continuous-time signals: (a) unit-step signal, (b) impulse signal, (c) boxcar signal,
and (d) sinusoidal signal.
n =0
Example 1.2. Find the periods of signals: x1 (t) = sin(2πt/36), x2 (t) = cos(4πt/15 +
2), x3 (t) = exp( j0.1t), x4 (t) = x1 (t) + x2 (t), and x5 (t) = x1 (t) + x3 (t).
⋆Periods are calculated according to (1.6). For x1 (t) the period follows
from 2πT1 /36 = 2π, as T1 = 36. Similarly, T2 = 15/2 and T3 = 20π. The
period of x4 (t) is the smallest interval containing T1 and T2 . It is T4 = 180 (5
periods of x1 (t) and 24 periods of x2 (t)). For signal x5 (t), when the periods of
components are T1 = 36 and T3 = 20π, there is no common interval T5 such
that the periods T1 and T3 are contained an integer number of times. Thus,
the signal x5 (t) is not periodic.
• Signal energy
"∞
Ex = | x (t)|2 dt, (1.9)
−∞
"T
1
PAV = lim | x (t)|2 dt.
T →∞ 2T
−T
if the Dirichlet conditions are met: (1) the signal x (t) has a finite number of
discontinuities within the period T; (2) it has a finite average value in the
period T; and (3) the signal has a finite number of maxima and minima.
Since the signal analysis deals with real-world physical signals, rather than
with mathematical generalizations, these conditions are almost always met.
Ljubiša Stanković Digital Signal Processing 25
# $ 1 T/2"
e j2πmt/T , e j2πnt/T = e j2πmt/T e− j2πnt/T dt
T
− T/2
%
1 for m = n
= sin(π (m−n)) .
π (m−n)
= 0 for m ̸= n
It means that the inner product of any two different basis functions is zero
(orthogonal set), while the inner product of a function with itself is 1 (normal
set). In the case of orthonormal set of basis functions, it is easy to show that
the weighting coefficients Xn can be calculated as projections of x (t) onto
the basis functions e j2πnt/T ,
# $ 1 T/2
"
Xn = x (t), e j2πnt/T = x (t)e− j2πnt/T dt. (1.12)
T
− T/2
This relation follows after a simple multiplication of the right and left sides
& T/2
of (1.11) by e− j2πmt/T and an integration within the period T1 −T/2 (·) dt.
Example 1.3. Show that the Fourier series coefficients Xn of a periodic signal x (t)
can be obtained by minimizing the mean square error between the signal and
∑nN=− N Xn e j2πnt/T within the period T.
⋆The mean square value of error
N
e(t) = x (t) − ∑ Xn e j2πnt/T ,
n=− N
From ∗
∂I/∂Xm = 0 follows
T/2
"
( )
N
1
e− j2πmt/T x (t) − ∑ Xn e j2πnt/T dt = 0
T n=− N
− T/2
T/2
"
1
Xm = x (t)e− j2πmt/T dt. (1.13)
T
− T/2
26 Continuous-Time Signals and Systems
∂F (z) ∂ | f (z)|2
∗ = = 2e− jα f (z).
∂z ∂z∗
In our case
' '
' '
| f (z)|2 = 'a + jb + e jα ( x + jy)'
= ( a + x cos α − y sin α)2 + (b + x sin α + y cos α)2
For the minimization of a function of two variables x and y we need partial
derivatives
∂ | f (z)|2
= 2 cos α( a + x cos α − y sin α)+ (1.14)
∂x
2 sin α(b + x sin α + y cos α)
= 2 Re{e− jα f (z)}
and
∂ | f (z)|2
= 2 Im{e− jα f (z)}. (1.15)
∂y
Therefore, all calculations with two real-valued equations (1.14) and (1.15)
are the same as using one complex valued relation
* +
∂ | f (z)|2 ∂ | f (z)|2 ∂ ∂ ∂F (z)
+j = +j F ( x, y) = .
∂x ∂y ∂x ∂y ∂z∗
Since the signal and the basis functions are periodic with period T, in
all previous integrals, we can use
T/2
" " +Λ
T/2
1 − j2πnt/T 1
x (t)e dt = x (t)e− j2πnt/T dt (1.16)
T T
− T/2 − T/2+Λ
Ljubiša Stanković Digital Signal Processing 27
Example 1.5. Calculate the Fourier series coefficients of a periodic signal x (t)
defined as
∞
x (t) = ∑ x0 (t + 2n)
n=−∞
with
x0 (t) = u(t + 1/4) − u(t − 1/4). (1.18)
x(t) Xn
-1 -1/4 1/4 1
t 0 n
Figure 1.2 Periodic signal (left) and its Fourier series coefficients (right).
1.5 1.5
1 1
x (t)
x (t)
0.5 0.5
1
0 0
(a) (b)
-0.5 -0.5
-2 -1 0 1 2 -2 -1 0 1 2
t t
1.5 1.5
1 1
x (t)
x (t)
0.5 0.5
30
6
0 0
(c) (d)
-0.5 -0.5
-2 -1 0 1 2 -2 -1 0 1 2
t t
Figure 1.3 Illustration of signal reconstruction by using a finite Fourier series with: (a)
coefficients Xn within −1 ≤ n ≤ 1, (b) coefficients Xn within −2 ≤ n ≤ 2, (c) coefficients Xn
within −6 ≤ n ≤ 6, and (d) coefficients Xn within −30 ≤ n ≤ 30.
For a real-valued signal x (t) the Fourier series coefficients can be written as
T/2
" T/2
"
1 2πnt 1 2πnt
Xn = x (t) cos( )dt − j x (t) sin( )dt
T T T T
− T/2 − T/2
An − jBn
= . (1.20)
2
Ljubiša Stanković Digital Signal Processing 29
∗
where An /2 and − Bn /2 are real and imaginary part of Xn . Since Xn = X− n
holds for real-valued signals, the values of An and Bn are
T/2
"
2 2πnt
A n = Xn + X− n = x (t) cos( )dt,
T T
− T/2
T/2
"
Xn − X− n 2 2πnt
Bn = = x (t) sin( )dt. (1.21)
−j T T
− T/2
An = Hn + H−n
Bn = Hn − H−n .
The coefficients calculated by (1.24) are the Hartley series coefficients. For a
real-valued and even signal x (t) = x (−t) this transform reduces to
T/2
"
1 2πnt
Cn = x (t) cos( )dt
T T
− T/2
0.6 0.6
0.4 0.4
x1(t)
x (t)
2
0.2 0.2
0 0
(a) (b)
-1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1
t t
0.6 0.6
0.4 0.4
x (t)
x6(t)
30
0.2 0.2
0 0
(c) (d)
-1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1
t t
Figure 1.4 Reconstruction of a signal using the Fourier series. Reconstructed signal is denoted
by x M (t), where M indicates the number of coefficients used in reconstruction.
with X0 = 1/8. Note that the relation between the Fourier coefficients in (a)
(b) ( a)
and (b) is 2X2n = Xn . The reconstruction is presented in Fig.1.5.
(c) For the signal xc (t) extended with its reversed version follows
1/2
" "1
Cn = te− j2πnt dt + (1 − t)e− j2πnt dt
0 1/2
1/2
"
(−1)n − 1
=2 t cos(2πnt)dt =
2π 2 n2
0
0.6 0.6
0.4 0.4
x1(t)
x2(t)
0.2 0.2
0 0
(a) (b)
-1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1
t t
0.6 0.6
0.4 0.4
x30(t)
x6(t)
0.2 0.2
0 0
(c) (d)
-1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1
t t
Figure 1.5 Reconstruction of a periodic signal, with a zero interval extension before using the
Fourier series.
0.6 0.6
0.4 0.4
x1(t)
x2(t)
0.2 0.2
0 0
(a) (b)
-1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1
t t
0.6 0.6
0.4 0.4
x30(t)
x6(t)
0.2 0.2
0 0
(c) (d)
-1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1
t t
Figure 1.6 Reconstruction of a periodic signal after an even extension before using the Fourier
series (cosine Fourier series).
Ljubiša Stanković Digital Signal Processing 33
A system transforms one signal (input signal) into another signal (output
signal). Assume that x (t) is the input signal. The system transformation will
be denoted by an operator T {◦}. The output signal can be written as
A system is linear if, for any two signals x1 (t) and x2 (t) and arbitrary
constants a1 and a2 , it holds
for any t0 .
Linear time-invariant (LTI) systems are fully described by their re-
sponse to the impulse signal. If we know the impulse response of these
systems,
h(t) = T {δ(t)},
then for arbitrary signal x (t) at the input, the output can be calculated, by
using (1.3), as
⎧ ⎫
⎨ "∞ ⎬
y(t) = T { x (t)} = T x (τ )δ(t − τ )dτ
⎩ ⎭
−∞
"∞ "∞
Linearity Time−invariance
= x (τ )T {δ(t − τ )}dτ = x (τ )h(t − τ )dτ.
−∞ −∞
x ( t ) ∗ t h ( t ) = h ( t ) ∗ t x ( t ). (1.29)
Example 1.7. Find a convolution of signals x (t) = u(t + 1) − u(t − 1) and h(t) =
e − t u ( t ).
⋆By using the convolution definition, we get
"∞ "1
y(t) = x (τ )h(t − τ )dτ = 1 · e−(t−τ ) u(t − τ )dτ
−∞ −1
⎧ & t +1
⎪
⎪ e−λ dλ = e−t (e − 1/e), for t ≥ 1
t −1
" ⎨ t −1
⎪
& t +1 − λ
=− e−λ u(λ)dλ = −(t+1) ,
⎪ 0 e dλ = 1 − e
⎪
for − 1 ≤ t < 1
t +1 ⎪
⎩
0 for t < −1.
"∞
|h(τ )|dτ < ∞ (1.30)
−∞
since
"∞ "∞
|y(t)| = | x (t − τ )h(τ )dτ | ≤ | x (t − τ )h(τ )|dτ
−∞ −∞
"∞ "∞
= | x (t − τ )|| h(τ )|dτ ≤ Mx |h(τ )|dτ < ∞,
−∞ −∞
if (1.30) holds.
It can be shown that the absolute value integrability of the impulse
response is the necessary condition for a linear time-invariant system to be
stable as well.
Ljubiša Stanković Digital Signal Processing 35
The Fourier series has been introduced and presented for periodic signals,
with a period T. Assume now that we extend the period to infinity, while not
changing the signal. This case corresponds to the analysis of an aperiodic
signal x (t). Its transform, the Fourier series coefficients normalized by the
period, is given by
T/2
" "∞
− j2πnt/T
lim Xn T = lim x (t)e dt = x (t)e− jΩt dt (1.31)
T →∞ T →∞
− T/2 −∞
"∞
X (Ω) = x (t)e− jΩt dt, (1.32)
−∞
is called the Fourier transform (FT) of a signal x (t). For the Fourier trans-
form existence it is sufficient that a signal is absolutely integrable. There
are some signals that do not satisfy this condition, whose Fourier transform
exists in a form of generalized functions, such as delta function.
The inverse Fourier transform (IFT) can be obtained by multiplying
both sides of (1.32) by e jΩτ and integrating over Ω,
"∞
1
x (t) = X (Ω)e jΩt dΩ. (1.33)
2π
−∞
36 Continuous-Time Signals and Systems
Example 1.8. Calculate the Fourier transform of x (t) = Ae−at u(t), a > 0.
⋆According to the Fourier transform definition we get
"∞
A
X (Ω) = Ae− at e− jΩt dt = .
( a + jΩ)
0
where
"0 "∞
2Ω
Xa (Ω) = −e at e− jΩt dt + e− at e− jΩt dt = . (1.35)
ja2 + jΩ2
−∞ 0
It results in
2
X (Ω) = . (1.36)
jΩ
Example 1.10. Find the Fourier transform of δ(t), x (t) = 1 and u(t).
⋆The Fourier transform of δ(t) is
"∞
FT{δ(t)} = δ(t)e− jΩt dt = 1. (1.38)
−∞
FT{1} = 2πδ(Ω).
Finally,
! 6
sign(t) + 1 1
FT{u(t)} = FT = + πδ(Ω). (1.39)
2 jΩ
"∞
y(t) = x (t) ∗t h(t) = Ae j(Ω0 (t−τ )+ ϕ) h(τ )dτ
−∞
"∞
= Ae j(Ω0 t+ ϕ) h(τ )e− jΩ0 τ dτ = H (Ω0 ) x (t), (1.40)
−∞
where
"∞
H (Ω) = h(t)e− jΩt dt (1.41)
−∞
is the Fourier transform of h(t). The linear time-invariant system does not
change the form of an input complex harmonic signal Ae j(Ω0 t+ ϕ) . It remains
complex harmonic signal after passing through the system, with the same
frequency Ω0 . The amplitude of the input signal x (t) is changed for | H (Ω0 )|
and the phase is changed for arg { H (Ω0 )}.
1. Linearity
where X1 (Ω) and X2 (Ω) are the Fourier transforms of signals x1 (t) and
x2 (t), separately.
2. Realness
The Fourier transform of a signal is real (i.e., X ∗ (Ω) = X (Ω)), if
x ∗ (−t) = x (t),
since
"∞ "∞
t→−t
X ∗ (Ω) = x ∗ (t)e jΩt dt = x ∗ (−t)e− jΩt dt = X (Ω), (1.43)
−∞ −∞
if x ∗ (−t) = x (t).
3. Modulation
"∞
FT{ x (t)e jΩ0 t } = x (t)e jΩ0 t e− jΩt dt = X (Ω − Ω0 ) (1.44)
−∞
FT{2x (t) cos(Ω0 t)} = X (Ω − Ω0 ) + X (Ω + Ω0 ).
4. Shift in time
"∞
FT{ x (t − t0 )} = x (t − t0 )e− jΩt dt = X (Ω)e− jt0 Ω . (1.45)
−∞
5. Time-scaling
"∞
1 Ω
FT{ x ( at)} = x ( at)e− jΩt dt = X ( ). (1.46)
| a| a
−∞
6. Convolution
"∞ "∞
FT{ x (t) ∗t h(t)} = x (τ )h(t − τ )e− jΩt dτdt (1.47)
−∞ −∞
"∞ "∞
t−τ →u
= x (τ )h(u)e− jΩ(τ +u) dτdu = X (Ω) H (Ω).
−∞ −∞
Ljubiša Stanković Digital Signal Processing 39
7. Multiplication
"∞ "∞
1
FT{ x (t)h(t)} = x (t) H (θ )e jθt dθe− jΩt dt (1.48)
2π
−∞ −∞
"∞
1
= H (θ ) X (Ω − θ )dθ = X (Ω) ∗Ω H (Ω) = H (Ω) ∗Ω X (Ω).
2π
−∞
9. Differentiation
⎧ ⎛ ⎞⎫
! 6 ⎨d "∞ ⎬
dx (t) ⎝ 1
FT = FT X (Ω)e jΩt dΩ⎠ = jΩX (Ω). (1.50)
dt ⎩ dt 2π ⎭
−∞
10. Integration
The Fourier transform of
"t
x (τ )dτ
−∞
Then,
⎧ ⎫
⎨ "t ⎬
FT = FT{ x (t)}FT{u(t)} =
x (τ )dτ (1.51)
⎩ ⎭
−∞
* +
1 1
X (Ω) + πδ(Ω) = X (Ω) + πX (0)δ(Ω).
jΩ jΩ
40 Continuous-Time Signals and Systems
It can be written as
where Xh (Ω) is the Fourier transform of the Hilbert transform of the signal
x (t). From Example 1.9 with the signal x (t) = sign(t) and the duality prop-
erty of the Fourier transform pair, obviously the inverse Fourier transform
of sign(Ω) is j/(πt). Therefore, the analytic part of a signal, in the time
domain, reads as
"∞
j 1 x (τ )
x a (t) = x (t) + jxh (t) = x (t) + x (t) ∗t = x (t) + j dτ. (1.54)
πt π t−τ
−∞
p.v.
where p.v. stands for Cauchy principal value of the considered integral.
1.3.3 Relationship Between the Fourier Series and the Fourier Transform
Consider an aperiodic signal x (t), with the Fourier transform X (Ω). As-
sume that the signal is of a limited duration (i.e., x (t) = 0 for |t| > T0 /2).
Then,
T"0 /2
X (Ω) = x (t)e− jΩt dt. (1.55)
− T0 /2
The periodic signal x p (t) can be expanded into Fourier series with the
coefficients
T/2
"
1
Xn = x p (t)e− j2πnt/T dt. (1.56)
T
− T/2
Ljubiša Stanković Digital Signal Processing 41
T/2
" T"0 /2
− j2πnt/T
x p (t)e dt = x (t)e− jΩt dt|Ω=2πn/T
− T/2 − T0 /2
or
1
Xn =
X (Ω)|Ω=2πn/T . (1.57)
T
It means that the Fourier series coefficients are the samples of the Fourier
transform, divided by T. The only condition in the derivation of this relation
is that the signal duration is shorter than the period of periodic extension
(i.e., T > T0 ). The sampling interval in frequency is
2π 2π
∆Ω = , ∆Ω < .
T T0
It should be smaller than 2π/T0 , where T0 is the signal x (t) duration. This
is a form of the sampling theorem in the frequency domain. The sampling
theorem in the time domain will be discussed later.
In order to write the Fourier series coefficients in the Fourier transform
form, note that a periodic signal x p (t), formed by a periodic extension of
x (t) with period T, can be written as
∞ ∞
x p (t) = ∑ x (t + nT ) = x (t) ∗t ∑ δ(t + nT ). (1.58)
n=−∞ n=−∞
since
% ; "∞
∞ ∞
FT ∑ δ(t + nT ) = ∑ δ(t + nT )e− jΩt dt
n=−∞ n=−∞−∞
∞ * +
jΩnT 2π ∞ 2π
= ∑ e = ∑ δ Ω− T n .
T n=−
(1.60)
n=−∞ ∞
42 Continuous-Time Signals and Systems
When a signal
x (t) = A(t)e jφ(t) (1.61)
is not of a simple analytic form, it may be possible, in some cases, to obtain
an approximative expression for its Fourier transform by using the method
of stationary phase.
The method of stationary phase states that if the phase function φ(t) is
monotonous and the amplitude A(t) is sufficiently smooth function, then
<
"∞
2πj
A(t)e jφ(t) e− jΩt dt ≃ A(t0 )e jφ(t0 ) e− jΩt0 , (1.62)
|φ′′ (t0 )|
−∞
Ωi ( t ) = φ′ ( t )
Around this point the phase can be expanded into a Taylor series as
1
φ(t) − Ωt = [φ(t0 ) − Ωt0 ] + [φ′ (t0 ) − Ω] + φ′′ (t0 )t2 + ...
2
Ljubiša Stanković Digital Signal Processing 43
"∞ "∞
1 ′′ ( t 2
A(t)e j(φ(t)−Ωt) dt ∼
= A(t0 )e j(φ(t0 )−Ωt0 ) ej 2 φ 0 )t dt
−∞ −∞
where A(t) ∼
= A(t0 ) is also used. With
<
"∞
j 12 at2 2πj
e dt =
| a|
−∞
8πt0 + 10π = Ω
Ω − 10π
t0 =
8π
and
φ′′ (t0 ) = 8π. (1.63)
The amplitude of X (Ω) is
'< ' =
' 2π '' 2π
' 2 2
| X (Ω)| ≃ A(t0 ) ' ' = exp (−( t 0 − 1 ) t 0 )
' φ′′ (t0 ) ' 8π
( >* + ?* + )
1 Ω − 10π 2 Ω − 10π 2
= exp − −1 (1.64)
2 8π 8π
The signal, stationary phase approximation of the Fourier transform and the
numerical value of the Fourier transform amplitudes are shown in Fig.1.7
44 Continuous-Time Signals and Systems
1
x(t)
-1
-4 -3 -2 -1 0 1 2 3 4
t
1
0.5
0
-100 -80 -60 -40 -20 0 20 40 60 80 100
Ω
1
Numeric calculation
|X(Ω)|
0.5
0
-100 -80 -60 -40 -20 0 20 40 60 80 100
Ω
Figure 1.7 The signal (top), along with the stationary phase method approximation of its
Fourier transform and the Fourier transform obtained by a numeric calculation (bottom).
and * +
Ω (2N −2)/(2N −1)
φ′′ (t0 ) = 2N (2N − 1) a . (1.65)
2Na
The amplitude and phase of X (Ω), according to (1.62), are
' '
2 2
' 2π '
'
| X (Ω)| ≃ A (t0 ) ' ′′ ' (1.66)
φ ( t0 ) '
* +1/(2N −1) '' * +1/(2N −1) ''
Ω ' 2π Ω '
= A2 ( )' '
2Na ' (2N − 1)Ω 2aN '
* +1/(2N −1)
(1 − 2N ) Ω
arg { X (Ω)} ≃ φ(t0 ) − Ωt0 + π/4 = Ω + π/4
2N 2aN
the method of stationary phase states that if the Fourier transform phase
function θ (t) is monotonous and the amplitude B(t) is sufficiently smooth
function, then
<
"∞
1 jθ (Ω) jΩt jθ (Ω0 ) jΩ0 t j
x (t) = B(Ω)e e dΩ ≃ B(Ω0 )e e , (1.68)
2π 2π |θ ′′ (Ω0 )|
−∞
−θ ′ (Ω0 ) = t,
and
t g = −θ ′ (Ω)
is the group delay.
Example 1.13. Consider a system with transfer function
aΩ0 + b = t
t−b
Ω0 =
a
and
θ ′′ (Ω0 ) = − a.
The impulse response is
<
2 j
h(t)) ≃ exp(−Ω20 )e− jaΩ0 /2− jbΩ0 + jΩ0 t
2π |θ ′′ (Ω0 )|
* + =
t − b 2 j((t−b)2 /(2a)+π/4) 1
= exp(− )e .
a 2πa
Example 1.14. For a system with frequency response H (Ω) = | H (Ω)| e j0 the im-
pulse response is h(t). Find the impulse response of the systems with transfer
functions shown in Fig.1.8 with:
(a) Ha (Ω) = | H (Ω)| e− j4Ω ,
2
(b) Hb (Ω) = | H (Ω)|@e− j2πΩ , and A
3
(c) Hc (Ω) = | H (Ω)| 4 + 14 cos(2πΩ2 ) e j0 .
"∞
1
h a (t) = H (Ω)e− j4Ω e jΩt = h(t − 4).
2π
−∞
"∞
1 2
hb (t) = H (Ω)e− j2πΩ e jΩt dΩ.
2π
−∞
|Hb(jΩ)|
|Ha(jΩ)|
|H (jΩ)|
c
Ω Ω Ω
arg{Hb(jΩ)}
arg{Ha(jΩ)}
arg{H (jΩ)}
c
Ω Ω Ω
|ha(t)|=|h(t-4)|
|hb(t)|
|h (t)|
c
t t t
Figure 1.8 Frequency response of systems (amplitude, top row, and phase, middle row) with
corresponding impulse responses (amplitude - bottom row).
48 Continuous-Time Signals and Systems
"∞
X (s) = L{ x (t)} = x (t)e−st dt, (1.69)
−∞
if
lim e−(s+ a)t = 0
t→∞
or σ + a > 0, that is, σ > − a. Therefore, the region of convergence of this
Laplace transform is the region where σ > − a. The point s = − a is the pole of
Ljubiša Stanković Digital Signal Processing 49
"∞ "∞
−σt −σt − jΩt
FT{ x (t)e }= x (t)e e dt = x (t)e−st dt = X (s). (1.70)
−∞ −∞
γ"+ jT
1
x (t) = lim X (s)est ds
2πj T →∞
γ− jT
Since the Laplace transform will be used to describe linear systems de-
scribed by linear differential equations we will consider only the relation of
the signal derivatives with the corresponding forms in the Laplace domain.
In general the Laplace transform of the first derivative dx (t)/d(t) of a signal
x (t) is
"∞ "∞
dx (t) −st
e dt = s x (t)e−st dt = sX (s).
dt
−∞ −∞
This relation follows by integration in part of the first integral, with the
assumption that the values of x (t)e−st are zero as t → ±∞.
In many applications it has been assumed that the systems are causal
with corresponding causal signals used in calculations. In these cases x (t) =
0 for t < 0, i.e., x (t) = x (t)u(t). Then the so called one-sided Laplace trans-
form (unilateral Laplace transform) is used. Its definition is
"∞
X (s) = x (t)e−st dt.
0
d( x (t)u(t)) dx (t)
= u ( t ) + x (0 ) δ ( t ).
dt dt
"∞ "∞
dx (t) −st
e dt = x (t)e−st |0∞ + s x (t)e−st dt = sX (s) − x (0).
dt
0 0
"∞ n "∞
d x (t) −st
n
e dt = s n
x (t)e−st dt − sn−1 x (0) − sn−2 x ′ (0) − ... − x (n−1) (0)
dt
0 0
= s X (s) − sn−1 x (0) − sn−2 x ′ (0) + ... − x (n−1) (0).
n
"t
1
L{ x (τ )dτ } = L{u(t) ∗t x (t)} = X (s)},
s
0
&∞
since L{u(t)} = 0 e−st dt = 1/s.
The initial and final values of the signal are x (0) = lims→∞ sX (s) and
x (∞) = lims→0 sX (s), respectively.
After we have established the relation between the Laplace transform and
signals derivatives we may use it to analyze the systems described by
differential equations. Consider a causal system
with the initial conditions x (0) = x ′ (0) = x (n−1) (0) = 0. The Laplace trans-
form of both sides of this differential equation is
Y (s) b s M + ... + b1 s + b0
H (s) = = M N .
X (s) a N s + ... + a1 s + a0
|H(jΩ)|2
|H(jΩ)|
Ω Ω Ω
Figure 1.9 Squared amplitude of the frequency response of a Butterworth filter of order N.
1
| H ( jΩ)|2 = B C2N .
Ω
1+ Ωc
1
H ( jΩ) H (− jΩ) = B C2N
jΩ
1+ jΩc
1
H (s) H (−s) = B C2N for s = jΩ.
s
1+ jΩc
54 Continuous-Time Signals and Systems
Im{s}
Im{s}
Re{s} Re{s} Re{s}
s0 , s1 , ..., s N −1
within the left side of the s plane, where Re{s} < 0, i.e., π/2 < αk < 3π/2.
The symmetric poles with Re {s} > 0 are the poles of H (−s). They are not
used in the filter design.
Example 1.18. Design a lowpass Butterworth filter with:
(a) N = 3 with Ωc = 1,
(b) N = 4 with Ωc = 3.
⋆(a) Poles for N = 3 with Ωc = 1 have the phases
2πk + π π
αk = + , for k = 0, 1, 2.
6 2
Ljubiša Stanković Digital Signal Processing 55
with
c 1
H (s) = √ √ =
(s + 1
−j 3 1 3 (s2 + s + 1)(s + 1)
2 2 )( s + 2 +j 2 )( s + 1)
2πk + π π
αk = + , for k = 0, 1, 2, 3.
8 2
with
c
H (s) =
(s2 + 2.296s + 9)(s2 + 5.543s + 9)
9
= 2
(s + 2.296s + 9)(s2 + 5.543s + 9)
In practice we usually do not know the filter order, but its passband
frequency Ω p and stopband frequency Ωs , with a maximal attenuation in
the passband a p [dB] and a minimal attenuation in the stopband a p [dB], as
shown in Fig.1.11. Based on these values we can calculate the order N and
the critical frequency Ωc needed for a filter design.
56 Continuous-Time Signals and Systems
1
A
p
|H(jΩ)|2
A
s
Ω Ω
p s
Figure 1.11 Specification of a Butterworth filter parameters in the passband and stopband.
1 2
B C2N ≥ A p (1.72)
Ωp
1+ Ωc
1 2
B C2N ≤ As .
Ωs
1+ Ωc
Nearest greater integer is assumed for the filter order N. Then we can
use any of the relations in (1.72) with equality sign to calculate Ωc . If we
' '2
choose the first one then Ωc will satisfy ' H ( jΩ p )' = A2p , while if we use
the second relation the value of Ωc will satisfy | H ( jΩs )|2 = A2s . These two
values differ. However both of them are within the defined criteria for the
transfer function.
The relation
a = 20 log A
or A = 10a/20 should be used for the attenuation given in [dB] .
All other filter forms, like passband and highpass, may be obtained
from a lowpass filter with appropriate signal modulations. These modula-
tions will be discussed for discrete-time filter forms in Chapter V.
Part II
57
Chapter 2
Discrete-Time Signals and Transforms
T
HE
tion in time. A continuous-time signal is converted into a sequence
of numbers, defining the discrete-time signal. The basic definitions
of discrete-time signals and their transforms are presented in this chapter.
The key fact in the conversion from a continuous-time signal into a sequence
of numbers is that these two signal representations are equivalent under cer-
tain conditions. The discrete-time signal may contain the same information
as the original continuous-time signal. The sampling theorem is fundamen-
tal for this relation between two signal forms. It is presented in this chapter,
after basic definitions of discrete-time signals and systems are introduced.
59
60 Discrete-Time Signals and Transforms
Δt t n
Figure 2.1 Signal discretization: continuous-time signal (left) and corresponding discrete-
time signal (right).
context will always be clear, so that there is no doubt what kind of signal is
considered. Notation x [n] is sometimes used in literature for discrete-time
signals, instead of x (n).
as illustrated in Fig.2.3.
The discrete unit-step signal is defined by
!
1, for n ≥ 0
u(n) = . (2.5)
0, for n < 0
Ljubiša Stanković Digital Signal Processing 61
1 1
x(n)=u(n)
δ(n)
0 0
-1 -1
(a) (b)
-10 0 10 -10 0 10
t n
1 1
x(n)=b(n)
sin(nπ/4)
0 0
-1 -1
(c) (d)
-10 0 10 -10 0 10
n n
Figure 2.2 Illustration of discrete-time signals: (a) unit-step function, (b) discrete-time im-
pulse signal, (c) boxcar signal b(n) = u(n + 2) − u(n − 3), and (d) discrete-time sinusoid.
4 4
2 2
-2 δ (n+2)
x(n)
0 0
-2 -2
-4 -4
-5 0 5 -5 0 5
n n
4 4
2 2
- δ(n-1 )
3δ(n)
0 0
-2 -2
-4 -4
-5 0 5 -5 0 5
n n
δ ( n ) = u ( n ) − u ( n − 1)
n
u(n) = ∑ δ ( k ).
k=−∞
x ( n + N ) = x ( n ). (2.7)
x (n) = x (−n).
where xe (n) and xo (n) are its even and odd part, respectively.
Ljubiša Stanković Digital Signal Processing 63
⋆For a signal x (n) we can form its even and odd part as
x (n) + x (−n)
xe (n) =
2
x (n) − x (−n)
xo (n) = .
2
Summing these two parts, the signal x (n) is reconstructed. Note that xo (0) =
0.
A signal is Hermitian if
x (n) = x ∗ (−n).
1 N # $
PAV = lim ∑ | x (n)|2 = | x (n)|2 , (2.9)
N →∞ 2N + 1 n=− N
# $
where | x (n)|2 is used to denote an average over large number of signal
values, as N → ∞. The average power of signals with a finite energy (energy
signals) is PAV = 0. For power signals (when 0 < PAV < ∞) the energy is
infinite, Ex → ∞.
Example 2.3. The energy of signal x (n) is Ex = 10. The energy of its even part is
Exe = 3. Find the energy of its odd part.
⋆The energy of signal is
∞ ∞
Ex = ∑ | x (n)|2 = ∑ | xe (n) + xo (n)|2
n=−∞ n=−∞
∞
= ∑ [ xe (n) + xo (n)][ xe (n) + xo (n)]∗
n=−∞
∞ ∞ ∞
= ∑ | xe (n)|2 + ∑ | xo (n)|2 + ∑ [ xo (n) xe∗ (n) + xe (n) xo∗ (n)].
n=−∞ n=−∞ n=−∞
64 Discrete-Time Signals and Transforms
The terms xo (n) xe∗ (n) and xe (n) xo∗ (n) in the last sum correspond to odd
signals
For the signals xe (n) and xo (n), satisfying the previous relation, we say that
they are orthogonal.
Therefore, for the energies Ex , Exe , and Exo , holds
Ex = Ex e + Ex o .
A discrete system T {·} is linear if for any two signals x1 (n) and x2 (n) and
any two constants a1 and a2 holds
holds
T { x (n − n0 )} = y(n − n0 ),
for any t0 .
For any input signal x (n) the signal at the output of a linear time-
invariant discrete system can be calculated if we know the output to the
impulse signal. The output to the impulse signal, h(n) = T {δ(n)}, is the
impulse response.
Ljubiša Stanković Digital Signal Processing 65
3 3
2 2
1 1
h(n)
x(n)
0 0
-1 -1
-2 -2
-4 -2 0 2 4 6 8 -4 -2 0 2 4 6 8
n n
Figure 2.4 Input signal and impulse response.
x ( n ) ∗ n h ( n ) = h ( n ) ∗ n x ( n ). (2.15)
Example 2.4. Calculate discrete-time convolution of signals x (n) and h(n) shown
in Fig. 2.4.
⋆By definition, according to Fig. 2.5, we have
∞
y (0) = ∑ x (k)h(−k) = 1 − 1 + 2 = 2,
k=−∞
∞
y (1) = ∑ x (k)h(1 − k ) = −1 − 1 + 1 + 4 = 3.
k=−∞
3 3
2 2
h(-k )
1 1
x(k)
0 0
-1 -1
-2 -2
-4 -2 0 2 4 6 8 -4 -2 0 2 4 6 8
n n
3 3
2 2
h(1- k)
h(2- k)
1 1
0 0
-1 -1
-2 -2
-4 -2 0 2 4 6 8 -4 -2 0 2 4 6 8
n n
Figure 2.5 Signals for the output y(0), y(1), and y(2) calculation.
8
6
4
y(n)
2
0
-2
-4 -2 0 2 4 6 8
n
Example 2.5. Calculate the convolution of signals x (n) = n[u(n) − u(n − 10)] and
h ( n ) = u ( n ).
⋆The convolution is
∞
y(n) = ∑ k[u(k) − u(k − 10)]u(n − k ) =
k=−∞
⎧ n
⎪
⎪ n +1
⎨ ∑ k=n 2 for 0≤n≤9
k =0
= ∑ k=
⎪ 9
0≤k ≤9 and k≤n ⎪
⎩ ∑ k = 45 for n>9
k =0
n+1
=n [u(n) − u(n − 10)] + 45u(n − 10).
2
Ljubiša Stanković Digital Signal Processing 67
Therefore |y(n)| < ∞ if (2.16) holds. It can be shown that the absolute sum
convergence of the impulse response is the necessary condition for a linear
time-invariant discrete system to be stable as well.
x (n∆t)∆t −→ x (n)
Ω∆t −→ ω, (2.19)
x (n) = Ae−α|n|
2.2.1 Properties
With respect to the signal shift and modulation the Fourier transform
of discrete-time signals behaves in the same way as the Fourier transform
of continuous-time signals,
and
FT{ x (n)e jω0 n } = X (e j(ω −ω0 ) ). (2.27)
Example 2.9. The Fourier transform of a discrete-time signal x (n) is X (e jω ).
Find the Fourier transform of y(n) = x (2n).
⋆For y(n) = x (2n) the Fourier transform is
∞
FT{ x (2n)} = ∑ x (2n)e− jωn
n=−∞
∞
x (n) + (−1)n x (n)
= ∑ e− jωn/2
n=−∞ 2
∞
1
= ∑ [x(n) + e− jnπ x(n)]e− jωn/2
2 n=− ∞
1 1
= [ X (e jω/2 ) + X (e j(ω/2+π ) )] = [ X (e jω/2 ) + X (e j(ω +2π )/2 )]. (2.28)
2 2
The period of this Fourier transform is 2π. Period of X (e jω/2 ) is 4π.
Example 2.10. Calculate the Fourier transform of the discrete-time signal (rectan-
gular window),
w R ( n ) = u ( N + n ) − u ( n − N − 1 ). (2.29)
Write the Fourier transform of a Hann(ing) window
1
w H (n) = [1 + cos(nπ/N )] [u( N + n) − u(n − N − 1)] .
2
⋆By definition
N
1 − e− jω (2N +1) sin(ω 2N2+1 )
WR (e jω ) = ∑ e− jωn = e jωN = . (2.30)
n=− N 1 − e− jω sin(ω/2)
1
N=4
W (ejω)
w (n)
0.5
R
R
0
-10 0 10 -π ω π
n
1
N=8
W (ejω)
w (n)
0.5
R
R
0
-10 0 10 -π ω π
n
1
N=8
W (ejω)
w (n)
0.5
H
-10 0 10 -π ω π
n
Figure 2.7 Discrete-time signal in a form of rectangular window of the widths 2 N + 1 = 9 and
2N + 1 = 17 samples (top and middle), and a Hann(ing) window with 2N + 1 = 17 (bottom).
The time domain values are on the left while the Fourier transforms of these discrete-time
signals are on the right.
As the window width increases in the time domain the main lobe width in
the Fourier domain is narrowing. The first zero value of the Fourier transform
of a rectangular window is at ω (2N + 1)/2 = π, i.e., at ω = 2π/(2N + 1)
where 2N + 1 is the signal duration. In the case of a Hann(ing) window
the main lobe is wider as compared to the rectangular window of the same
width, but its convergence is much faster with very reduced oscillations in
the Fourier transform, Fig.2.7.
72 Discrete-Time Signals and Transforms
∞ ∞
FT{ x (n) ∗n h(n)} = ∑ ∑ x (k )h(n − k )e− jnω = X (e jω ) H (e jω ), (2.33)
n=−∞ k=−∞
∞
H (e jω ) = ∑ h(n)e− jωn
n=−∞
Example 2.11. Find the output of a discrete linear time-invariant system with
frequency response H (e jω ) if the input signals are:
(a) x (n) = Ae jω0 n and (b) x (n) = A cos(ω0 n + ϕ). What is the output if
the impulse response h(n) is real-valued?
∞ ∞
y(n) = ∑ h(k ) x (n − k ) = ∑ h(k) Ae jω0 (n−k)
k=−∞ k=−∞
∞
= Ae jω0 n ∑ h(k) Ae− jω0 k = Ae jω0 n H (e jω0 )
k =−∞
' '
' ' jω
= A 'H (e jω0 )' e j(ω0 n+arg{ H (e 0 }) .
A j ( ω0 n + ϕ ) A − j ( ω0 n + ϕ )
x (n) = A cos(ω0 n + ϕ) = e + e .
2 2
A '' '
' jω0
y(n) = 'H (e jω0 )' e j(ω0 n+ ϕ)+ j arg{ H (e )}
2
A '' '
' − jω
+ 'H (e− jω0 )' e− j(ω0 n+ ϕ)+ j arg{ H (e 0 )} .
2
H (e jω ) = H ∗ (e− jω )
Ljubiša Stanković Digital Signal Processing 73
and
∞ ∞
H (e jω ) = ∑ h(n) cos(ωn) + j ∑ h(n) sin(ωn)
n=−∞ n=−∞
' '2 ' '2
' ' ' '
'H (e jω )' = 'H (e− jω )'
! ∞ 6
∑n=−∞ h(n) sin(ωn)
jω
arg{ H (e ) = arctan = − arg{ H (e− jω ).
∑∞ n=−∞ h (n ) cos( ωn )
∞ ∞ "π
1
∑ ∑ 2π X (e jω )e jωn y∗ (n)dω
x (n)y∗ (n) = (2.35)
n=−∞ n=−∞ −π
"π
( ) "π
∞
1 1
=
2π
jω
X (e ) ∑ (e y(n)) dω = 2π X (e jω )Y∗ (e jω )dω.
− jωn ∗
−π n=−∞ −π
1 ' '2
' '
Pxx (e jω ) = lim 'X N (e jω )' , (2.37)
N →∞ 2N + 1
N N
1
Pxx (e jω ) = lim ∑ ∑ x (n) x ∗ (m)e− jω (n−m) . (2.38)
N →∞ 2N + 1 n=− N m=− N
2N
1
Pxx (e jω ) = lim ∑ (2N + 1 − |k |)r (k )e− jωk
N →∞ 2N + 1
k=−2N
with a period 2Ω0 . It is very important to note that X p (Ω) = X (Ω) for
|Ω| < Ω0 if
Ω0 > Ω m .
In this case, it is possible to make transformation from X (Ω) to X p (Ω) and
back without losing any information.
Of course, that would not be the case if Ω0 > Ωm did not hold. By
periodic extension of X (Ω), in that case, overlapping (aliasing) would have
occurred in X p (Ω). It would not be reversible. Then it would not be possible
to recover X (Ω) from X p (Ω). The periodic extension is illustrated in Fig. 2.8.
The periodic function X p (Ω) can be expanded into Fourier series with
coefficients
"Ω0 "∞
1 jπΩn/Ω0 1
X− n = X p (Ω)e dΩ = X (Ω)e jπΩn/Ω0 dΩ.
2Ω0 2Ω0
− Ω0 −∞
The integration limits are extended to the infinity since X (Ω) = X p (Ω)
within the basic period interval and X (Ω) = 0 outside this interval.
76 Discrete-Time Signals and Transforms
X(Ω)
-Ω Ω Ω
m m
X (Ω) = X(Ω)
p
-Ω <Ω<Ω
0 0
- Ω0 - Ωm Ωm Ω0 Ω
Figure 2.8 The Fourier transform of a signal, with X (Ω) = 0 for |Ω| > Ωm (top) and its
periodically extended version, with period 2Ω0 > 2Ωm (bottom).
"∞
1
x (t) = X (Ω)e jΩt dΩ. (2.41)
2π
−∞
with ∆t = π/Ω0 .
Ljubiša Stanković Digital Signal Processing 77
"∞ "Ω0
1 jΩt 1
x (t) = X (Ω)e dΩ = X p (Ω)e jΩt dΩ (2.44)
2π 2π
−∞ − Ω0
"Ω0
( )
∞
1
= ∑ Xn e jπnΩ/Ω0 e jΩt dΩ
2π n=−∞
− Ω0
"Ω0
( )
∞
1 jπnΩ/Ω0
= ∑ x (−n∆t)∆te e jΩt dΩ
2π n=−∞
− Ω0
as
∞ π
sin( ∆t (t − n∆t))
x (t) = ∑ x (n∆t) π . (2.45)
n=−∞ ∆t (t − n∆t )
The signal x (t), for any t, is expressed in terms of its samples x (n∆t).
Example 2.13. The last relation can be used to prove that X (Ω) = X (e jω ) with
Ω∆t = ω and |ω | < π for the signals sampled at the rate satisfying the
sampling theorem.
⋆Starting from
"∞
X (Ω) = x (t)e− jΩt dt
−∞
the signal x (t), satisfying the sampling theorem, can by written in terms of
samples, according to the third row of (2.44), as
"Ω0
( )
∞
1 − j∆tnθ
x (t) = ∑ x (n∆t)∆te e jθt dθ.
2π n=−∞
− Ω0
It follows
"∞ "Ω0
( )
∞
1 − j∆tnθ
X (Ω) = ∑ x (n∆t)∆te e jθt dθe− jΩt dt
2π n=−∞
−∞ − Ω0
∞ "Ω0
= ∑ x (n∆t)∆t δ(θ − Ω)e− j∆tnθ dθ
n=−∞
− Ω0
∞
= ∑ x (n∆t)∆te− j∆tnΩ for |Ω| < Ω0 (2.46)
n=−∞
78 Discrete-Time Signals and Transforms
resulting in
∞
X (Ω) = ∑ x (n)e− jωn for |ω | < π
n=−∞
Example 2.14. If the highest frequency in a signal x (t) is Ωm1 and the highest
frequency in a signal y(t) is Ωm2 what should be the sampling interval for the
signal x (t)y(t) and for the signal x (t − t1 )y∗ (t − t2 )? The highest frequency
Ωm in a signal is used in the sense that the Fourier transform of the signal is
zero for |Ω| > Ωm .
2 2 2
X p (Ω) = ... + + + + ...
1 + (Ω + 20π )2 1 + Ω2 1 + (Ω − 20π )2
Ljubiša Stanković Digital Signal Processing 79
Thus, the value of X p (Ω) at the period ending points ±10π will approxi-
mately be X p (±10π ) = 2/(1 + 100π 2 ) ∼ = 0.002. Comparing with the maxi-
mum value X p (0) = 2, it means that the expected error due to the discretiza-
tion of this signal (since it does not strictly satisfy the sampling theorem) will
be of a 0.1% order.
(b) The discrete-time signal obtained by sampling x (t) = exp(− |t|)
with ∆t = 0.1 is x (n) = 0.1e−0.1|n| . Its Fourier transform is already calculated
with A = 0.1 and α = 0.1, eq.(2.22). The result is
1 − e−0.2
X (e jω ) = 0.1 . (2.48)
1 − 2e−0.1 cos(ω ) + e−0.2
1
y(n) = cos(nπ/4 + π/4)∆t
2
corresponding to the continuous-time signal
1 π 1
y(t) = cos(n + π/4) = cos(25πt + π/4).
2 4∆t 2
2.4 PROBLEMS
Problem 2.1. Check the periodicity and find the period of signals:
(a) x (n) = sin(2πn/32), (b) x (n) = cos(9πn/82), (c) x (n) = e jn/32 , and (d)
x (n) = sin(πn/5) + cos(5πn/6) − sin(πn/4).
Problem 2.2. Check the linearity and time-invariance of the discrete system
described by equation
y(n) = x (n) + 2.
h(n) = h1(n)*h2(n)*h2(n)
12 3
10 2
8
1
x(n)
6
4 0
2 -1
0
-2
-2 0 2 4 6 8 -2 0 2 4 6 8
n n
Figure 2.9 Problem 2.7, impulse response h(n) (left) and Problem 2.14, discrete signal x (n)
(right).
Problem 2.5. Find the convolution of signals x (n) = e−|n| and h(n) = u(n +
5) − u ( n − 6).
Problem 2.6. A discrete system consists of systems with impulse responses
h1 (n) = e− an u(n), h2 (n) = e−bn u(n), and h3 (n) = u(n). Find the impulse
response of the resulting system for:
(a) Systems h1 (n), h2 (n), and h3 (n) connected in parallel,
(b) System h1 (n) connected in parallel with a cascade of systems h2 (n)
and h3 (n).
Problem 2.7. Consider three causal linear time-invariant systems in cas-
cade. Impulse responses of these systems are h1 (n), h2 (n), and h2 (n), re-
spectively. The impulse response of the second and the third system is
h2 (n) = u(n) − u(n − 2), while the impulse response of the whole system,
h ( n ) = h 1 ( n ) ∗ n h 2 ( n ) ∗ n h2 ( n ),
Find
∞
S= ∑ ne−n/2 .
n =0
82 Discrete-Time Signals and Transforms
H (e jω ) ∼
= jω for small ω,
i.e., ' '
dH (e jω ) '' d2 H (e jω ) ''
= j and = 0.
dω 'ω =0 dω 2 'ω =0
Problem 2.11. Find the Fourier transform of the following discrete-time
signal (triangular window)
* +
|n|
wT (n) = 1 − [u(n + N ) − u(n − N − 1)].
N+1
with N being an even number.
Problem 2.12. Find the value of integral
"π
1 sin2 (( N + 1)ω/2)
I= dω.
2π sin2 (ω/2)
−π
w(n) = w H (n + N ) + w H (n) + w H (n − N )
Problem 2.14. A discrete-time signal x (n) is given in Fig. 2.9 (right). Without
calculating its Fourier transform X (e jω ) find
Using this Fourier transform find the center of gravity of signal x (n) =
e−n/4 u(n) defined by
∞
∑ nx (n)
n=−∞
ng = ∞ .
∑ x (n)
n=−∞
sin(nπ/3)
(a) h(n) = , with h(0) = 1/3,
nπ
sin2 (nπ/3)
(b) h(n) = ,
(nπ )2
sin((n − 2)π/4)
(c) h(n) = .
( n − 2) π
Show that the frequency response of the system with h(n) = sin(nπ/3)/nπ
is H (e jω ) = 1 for |ω | ≤ π/3 and H (e jω ) = 0 for π/3 < |ω | < π. Find the
frequency responses in other two cases. Find the systems output to the input
signal x (n) = sin(nπ/6).
where xh (n) is the Hilbert transform of x (n). Find the impulse response of
the system that transforms a signal x (n) into its Hilbert transform (Hilbert
transformer).
Problem 2.20. For a signal whose Fourier transform is zero for frequencies
Ω ≥ Ωm = 2π f m = π/∆t show that
"∞
sin(π (t − τ )/∆t)
x (t) = x (τ ) dτ.
π (t − τ )
−∞
Problem 2.21. Sampling of a signal is done twice, with the sampling interval
∆t = 2π/Ωm that is twice larger than the sampling interval required by the
sampling theorem (∆t = π/Ωm is required). After first sampling process,
the discrete-time signal x1 (n) = ∆tx (n∆t) is formed, while after the second
sampling process signal x2 (n) = ∆tx (n∆t + a) is formed. Show that we can
reconstruct continuous-time signal x (t) based on x1 (n) and x2 (n) if a ̸= k∆t,
that is, if samples x1 (n) and x2 (n) do not overlap in continuous-time.
Problem 2.23. Show that the relation among the amplitudes of a signal
x (n) and its even and odd parts xe (n) = [ x (n) + x (−n)]/2 and xo (n) =
[ x (n) − x (−n)]/2 is
√
As (n) ≤ | xe (n)| + | xo (n)| ≤ 2As (n)
@ A
with As (n) > 0 defined by A2s (n) = | x (n)|2 + | x (−n)|2 /2.
2.5 SOLUTIONS
T { x ( n − N ) = x ( n − N ) + 2 = y ( n − N ).
h(n) = T {δ(n)}.
86 Discrete-Time Signals and Transforms
It can be written as
h(n) = T {u(n) − u(n − 1)}.
∞ ∞
2−1
∑ |h(n)| = 1 + ∑ 2− n = 1 + = 2.
n=−∞ n =1 1 − 2−1
The system is stable since the sum of absolute values of impulse response is
finite.
∞
y (0) = ∑ x (k ) x (−k ) = x (0) x (0) = 1
k=−∞
∞
y (1) = ∑ x ( k ) x (1 − k ) = x (0 ) x (1 ) + x (1 ) x (0 ) = 2
k=−∞
∞
y(−1) = ∑ x (k ) x (−1 − k ) = 0
k=−∞
∞
y (2) = ∑ x ( k ) x (2 − k ) = 3
k=−∞
...
1.5 1.5
1 1
x(-k )
x(k)
0.5 0.5
0 0
-0.5 -0.5
-15 -10 -5 0 5 10 15 -15 -10 -5 0 5 10 15
k k
1.5 1.5
1 1
x(1-k )
x(2-k )
0.5 0.5
0 0
-0.5 -0.5
-15 -10 -5 0 5 10 15 -15 -10 -5 0 5 10 15
k k
1.5 6
x(n)* x(n)
1
x(-1-k )
4
0.5
2
0
-0.5 0
-15 -10 -5 0 5 10 15 -15 -10 -5 0 5 10 15
k n
with !
1, for k ≤ n + 5
u((n − k ) + 5) =
0, for k > n + 5
and !
1, for k ≤ n − 6
u((n − k ) − 6) =
0, for k > n − 6
we get
!
1, for n − 6 < k ≤ n + 5
(u((n − k) + 5) − u((n − k) − 6)) =
0, elsewhere.
88 Discrete-Time Signals and Transforms
Since !
k, for k ≥ 0
|k | = ,
−k, for k < 0
we have three cases:
1) For n + 5 ≤ 0, i.e., n ≤ −5, we have k ≤ 0 for all terms. Therefore |k | = −k,
n +5
1 − e11 e −5 − e 6
y(n) = ∑ e k = e n −5 = en
k = n −5
1−e 1−e
e0.5 e−5.5 − e5.5 sinh 5.5
= en 0.5 −0.5 = en
e e − e0.5 sinh 0.5
2) For n − 5 ≥ 0, the lowest k = n − 5 is greater than 0. Then k ≥ 0 for all
terms and |k | = k with
n +5
1 − e−11 5
−n e − e
−6
y(n) = ∑ e − k = e − n +5 = e
k = n −5
1 − e −1 1 − e −1
e−0.5 e5.5 − e−5.5 sinh 5.5
= e−n −0.5 0.5 = e−n .
e e − e−0.5 sinh 0.5
3) For −5 < n < 5, index k can assume positive and negative values. The
convolution is split into two sums as
n +5 −1 n +5 5− n n +5
y(n) = ∑ e−|k| = ∑ ek + ∑ e−k = ∑ e−k + ∑ e−k
k = n −5 k = n −5 k =0 k =1 k =0
1−e −( 5 − n ) 1 − e n +6)
−(
= e −1 + =
1 − e −1 1 − e −1
1 − e n −5 1/2 1 − e
−(n+6)
= e−1/2 1/2 + e
e − e−1/2 e1/2 − e−1/2
1
= 0.5 (e−0.5 − en−5.5 + e0.5 − e−n−5.5 ) =
e − e−0.5
−e−5.5 (en + e−n ) + e−0.5 + e0.5
=
e0.5 − e−0.5
cosh 0.5 − e−5.5 cosh(n)
= .
sinh 0.5
Ljubiša Stanković Digital Signal Processing 89
where
∞
h23 (n) = ∑ h 3 ( m ) h2 ( n − m ) = h 2 ( n ) ∗ h 3 ( n ).
m=−∞
The impulse response of the whole system is
h(n) = h1 (n) + h23 (n) = h1 (n) + h2 (n) ∗ h3 (n),
with
∞
h2 ( n ) ∗ h3 ( n ) = ∑ e−b(n−m) u (n − m )u (m )
m=−∞
n
1 − e b ( n +1) e−bn − eb
= u(n) ∑ e−b(n−m) = e−bn 1 − eb
u ( n ) =
1 − eb
u ( n ).
m =0
90 Discrete-Time Signals and Transforms
It follows
∞
e−(1/2+ jω )
H (e jω ) = ∑ ne−n/2 e− jωn = .
n =0 (1 − e−(1/2+ jω ) )2
The output for a real-valued h(n) is
' '
' '
y(n) = 5 'H (e jπ/10 )' sin(πn/5 + arg{ H (e jπ/10 })
' ' ' '
' ' ' '
− 3 'H (e jπ/6 )' cos(πn/3 + π/6 + 'H (e jπ/6 )')
=14.1587 sin(πn/5 − 1.1481)
− 5.7339 cos(πn/3 + π/6 − 1.6605).
Ljubiša Stanković Digital Signal Processing 91
Solution 2.10. For the impulse response h(n) the frequency response is
a + 2b = 1/2
a + 4b = 0
1
h(n) = δ(n + 1) − δ(n − 1) − (δ(n + 2) − δ(n − 2)).
4
Solution 2.11. Note that
1
wT (n) = w R (n) ∗n w R (n)
N+1
1 1 sin2 (ω N2+1 )
WT (e jω ) = WR (e jω )WR (e jω ) = .
N+1 N + 1 sin2 (ω/2)
Ljubiša Stanković Digital Signal Processing 93
sin(ω N2+1 )
X (e jω ) = .
sin(ω/2)
This signal is the rectangular window, x (n) = u(n + N/2) − u(n − N/2 − 1).
Its energy is
1
w H (n) = [1 + cos(nπ/N )] [u( N + n) − u(n − N − 1)] .
2
w(n) = w H (n) + w H (n − N )
1 1
= [1 + cos(nπ/N )] + [1 + cos((n − N )π/N )]
2 2
1 1
= 1 + cos(nπ/N ) + cos(nπ/N − π ) = 1.
2 2
The same holds for − N ≤ n ≤ −1 when
w(n) = w H (n + N ) + w H (n) = 1.
window WH (e jω ), is
For
K
w(n) = ∑ w H (n + kN )
k=−K
we get
⎧
⎪
⎪ 0D E for n < −(K + 1) N
⎪
⎪ π
⎨ 12 1 + cos((n + KN ) N ) for −(K + 1) N + 1 ≤ n ≤ −KN + 1
w(n) = 1D E for −KN ≤ n ≤ KN − 1
⎪
⎪ 1 π
⎪
⎪ 2 1 + cos((n − KN ) N ) for KN ≤ n ≤ (K + 1) N − 1
⎩
0 for n > ( K + 1) N − 1
with
K
1 − e− jω (2K +1) N
W (e jω ) = WH (e jω ) ∑ e− jωkN = e jωKN
k=−K 1 − e− jωN
sin(ω (2K + 1) N/2)
= WH (e jω ) .
sin(ωN/2)
Similar results hold for the Hamming and triangular window. The results
can be generalized for shifts of N/2, N/4,...
For very large K the second term variations in W (e jω ) are much faster
than the variations of WH (e jω ). Thus, for large K the Fourier transform
W (e jω ) approaches to the Fourier transform of a rectangular window of the
width (2K + 1) N.
Solution 2.14. Based on the definition of the Fourier transform of discrete-
time signals,
∞
X (e j0 ) = ∑ x (n) = 7,
n=−∞
∞
X (e jπ ) = ∑ x (n)(−1)n = 1,
n=−∞
"π
X (e jω )dω = 2πx (0) = 4π,
−π
Ljubiša Stanković Digital Signal Processing 95
1B C
Re{ X (e jω )} = X (e jω ) + X ∗ (e jω ) .
2
The inverse Fourier transform of Re { X (e jω )} is
1
y(n) = ( x (n) + x ∗ (−n)).
2
Solution 2.15. The Fourier transform of y(n) is
> ?
∞ ∞
d
Y (e jω ) = ∑ ne−n/4 u(n)e− jωn = j ∑ e−n/4 e− jωn
n=−∞ dω n =0
d 1 e − 1/4 − jω
=j = .
dω 1 − e 1/4− jω
− (1 − e−1/4− jω )2
is
π/3
" 'π/3
1 e jωn '' sin(πn/3)
h(n) = e jωn dω = = .
2π 2jπn '−π/3 πn
−π/3
The value of frequency response at the input signal frequency ω = ±π/6 is
H (e± jπ/6 ) = 1. The output signal is, y(n) = sin(nπ/6).
(b) The frequency response, in this case, is H (e jω ) ∗ω H (e jω ), resulting in
y(n) = 0.25 sin(nπ/6).
(c) Output signal in this case is y(n) = sin((n − 2)π/6) = sin(nπ/6 − π/3).
96 Discrete-Time Signals and Transforms
π ∞
X (e jω ) = ∑ [δ(ω − 0.2π + 2kπ )e jπ/4 + δ(ω + 0.2π + 2kπ )e− jπ/4 ]
100 k=− ∞
∞
π
+ ∑
j100 k=−∞
[δ(ω − 0.9π + 2kπ ) − δ(ω + 0.9π + 2kπ )].
π ∞
X (e jω ) = ∑ [δ(ω − 0.4π + 2kπ )e jπ/4 + δ(ω + 0.4π + 2kπ )e− jπ/4 ]
50 k=− ∞
π ∞
+ ∑ [δ(ω − 1.8π + 2kπ ) − δ(ω + 1.8π + 2kπ )].
j50 k=− ∞
Ljubiša Stanković Digital Signal Processing 97
H(ejω), X(ejω)
1
H(ejω), X(ejω)
1
jω jω
H(e ), X(e )
1
Figure 2.11 Illustration of the system output with various sampling intervals (a)-(c).
π
X (e jω ) = [δ(ω − 0.4π )e jπ/4 + δ(ω + 0.4π )e− jπ/4 ]
50
π
+ [δ(ω − 1.8π + 2π ) − δ(ω + 1.8π − 2π )]
j50
π
= [δ(ω − 0.4π )e jπ/4 + δ(ω + 0.4π )e− jπ/4 ]
50
π
+ [δ(ω + 0.2π ) − δ(ω − 0.2π )].
j50
x(n)
Figure 2.12 Illustration of the aliasing caused frequency change, from signal sin (90πt) to
signal − sin(10πt).
H(ejω) h(n)
1 2/π
-2 π -π 0 π 2π ω
0 n
Figure 2.13 Frequency and impulse response of the discrete-time Hilbert transformer.
1.5 1.5
1 1
X (Ω)
X(Ω)
p
0.5 0.5
0 0
-3 -2 -1 0 1 2 3 4 5 6 7 -3 -2 -1 0 1 2 3 4 5 6 7
Ω/Ω1 Ω/Ω1
Figure 2.14 Problem 2.19: illustration of the Fourier transform periodic extension.
∞
sin(π (t − n∆t)/∆t)
x (t) = e j4Ω1 t ∑ x (n∆t) with ∆t = π/Ω1 .
n=−∞ π (t − n∆t)/∆t
100 Discrete-Time Signals and Transforms
Solution 2.20. For signal whose Fourier transform is zero for frequencies
Ω ≥ Ωm = 2π f m = π/∆t hods
where !
1 for |Ω| < π/∆t
H (Ω) = .
0 for |Ω| ≥ π/∆t
The impulse response of H (Ω) is
1 π/∆t
& jΩt sin(πt/∆t)
h(t) = e dΩ = .
2π −π/∆t πt
"∞ "∞
sin(π (t − τ )/∆t)
x (t) = x (τ )h(t − τ )dτ = x (τ ) dτ.
π (t − τ )
−∞ −∞
Relation (1.60)
* + % ; % ;
2π ∞ 2π ∞ ∞
∑ δ Ω − ∆t k = FT ∑ δ(t + n∆t) = FT ∑ δ(t − n∆t)
∆t k=− ∞ n=−∞ n=−∞
is used.
Ljubiša Stanković Digital Signal Processing 101
"∞ ∞
x (t) = x p (t) ∗t h(t) = x (τ ) ∑ δ(τ − n∆t)h(t − τ )∆tdτ
−∞ n=−∞
∞ ∞ π
sin( ∆t (t − n∆t))
= ∑ x (n∆t)h(t − n∆t)∆t = ∑ x (n∆t) π . (2.52)
n=−∞ n=−∞ ∆t ( t − n∆t )
∆Ωm
with a reduction of the sampling interval to ∆t = π/(Ωm + 2 ) with
respect to ∆t = π/Ωm .
Solution 2.21. The Fourier transforms of discrete-time signals, in continu-
ous frequency notation, are periodically extended versions of X (Ω) with the
period 2π/∆t,
∞
X1 ( Ω ) = ∑ X (Ω + 2πn/∆t),
n−−∞
∞
X2 ( Ω ) = ∑ X (Ω + 2πn/∆t)e j(Ω+2πn/∆t)a .
n−−∞
H(Ω)
X(Ω)
- Ωm Ωm Ω
X (Ω) = X(Ω)
p
H(Ω) -Ω <Ω<Ω
0 0
X(Ω)
-Ω -Ω Ω Ω Ω
0 m m 0
Xp(Ω) = X(Ω)
- Ω0 < Ω < Ω0
X(Ω)
-Ω -Ω Ω Ω Ω
0 m m 0
Figure 2.15 Smoothed filter in the sampling theorem illustration (first two graphs) versus
original sampling theorem relation within filtering framework.
Similarly for negative frequencies, within the basic period −Ωm < Ω < 0,
follows
X (Ω)e j2πa/∆t − X2 (Ω)e− jΩa
X (Ω) = 1 for a ̸= k∆t.
e j2πa/∆t − 1
Ljubiša Stanković Digital Signal Processing 103
with * +
1 x (t0 + ∆t) + x (t0 − ∆t)
Ω0 = arccos .
∆t 2x (t0 )
The condition for a unique solution is that the argument of cosine is 0 ≤
Ω0 ∆t ≤ π, limiting the approach to small values of ∆t.
In addition, here we will discuss the discrete complex-valued signal.
For a complex sinusoid x (n) = A exp( j2πk0 n/N + φ0 ), with available two
samples x (n1 ) = A exp( jϕ(n1 )) and x (n2 ) = A exp( jϕ(n2 )), from
x ( n1 )
= exp( j2πk0 (n1 − n2 )/N )
x ( n2 )
follows
2πk0 (n1 − n2 )/N = ϕ(n1 ) − ϕ(n2 ) + 2kπ,
where k is an arbitrary integer. Then
ϕ ( n1 ) − ϕ ( n2 ) k
k0 = N+ N. (2.53)
2π (n1 − n2 ) n1 − n2
k0 = 5 + 16k/4,
| x (n)|2 + | x (−n)|2
| xe (n)|2 + | xo (n)|2 = = A2s (n).
2
Obviously
F | xe (n)|2 ≤ A2s (n) and | xo (n)|2 ≤ A2s (n). Replacing | xo (n)| =
A2s (n) − | xe (n)|2 into | xe (n)| + | xo (n)| we get
F
| xe (n)| + | xo (n)| = | xe (n)| + A2s (n) − | xe (n)|2 .
2.6 EXERCISE
Exercise 2.1. Calculate the convolution of signals x (n) = n[u(n) − u(n − 3)]
and h(n) = δ(n + 1) + 2δ(n) − δ(n − 2).
Exercise 2.2. Find the convolution of signals x (n) = e−|n| and h(n) = u(3 −
n ) u (3 + n ).
Exercise 2.3. The output of a linear time-invariant discrete system to the
input signal x (n) = u(n) is y(n) = ( 31n + n)u(n). Find the impulse response
h(n). Is the system stable?
Exercise&2.4. For signal x (n&) = nu(5 − n)u(n + 5) find the values of X (e j0 ),
π π
X (e jπ ), −π X (e jω )dω, and −π | X (e jω )|2 dω without the Fourier transform
calculation. Check the results by calculating the Fourier transform.
Exercise 2.5. For a signal x (n) at an instant m a signal y(n) = x (m −
n) x ∗ (m + n) is formed. Show that the Fourier transform of y(n) is real-
valued. What is the Fourier transform of y(n) if x (n) = A exp( jan2 /4 +
j2ω0 n)? Find the Fourier transform of z(m) = x (m − n) x ∗ (m + n) for a given
n.
Note: The Fourier transform of y(n) is the Wigner distribution of x (n)
for a given m, while the Fourier transform of z(m) is the Ambiguity function
of x (n) for a given n.
Exercise 2.6. For a signal x (n) with Fourier transform X (e jω ) find the
Fourier transform of x (2n). Find the Fourier transform of y1 (2n) = x (2n)
and y1 (2n + 1) = 0. What is the Fourier transform of x (2n + 1) and what is
the Fourier transform of y2 (2n) = 0 and y2 (2n + 1) = x (2n + 1). Check the
result by showing that Y1 (e jω ) + Y2 (e jω ) = X (e jω ).
Exercise 2.7. For a real-valued signal find the relation between the Fourier
transform of signal X (e jω ) and the Hartley transform
∞
H (e jω ) = ∑ x (n)[cos(ωn) + sin(ωn)].
n=−∞
Write this relation if the signal is real-valued and even, x (n) = x (−n).
Exercise 2.8. Systems with impulse responses h1 (n), h2 (n) and h3 (n) are
connected in cascade. If the impulse responses h2 (n) = h3 (n) = u(n) −
u(n − 2) and the resulting impulse response is h(n) = δ(n) + 5δ(n − 1) +
10δ(n − 2) + 11δ(n − 3) + 8δ(n − 4) + 4δ(n − 5) + δ(n − 6). Find the impulse
response h1 (n).
106 Discrete-Time Signals and Transforms
for k = 0, 1, 2, ..., N − 1.
In order to establish the relation between the DFT with the Fourier
transform of discrete-time signals, consider a discrete-time signal x (n) of
limited duration. Assume that nonzero samples of x (n) are within 0 ≤ n ≤
N0 − 1. Its Fourier transform is
N0 −1
X (e jω ) = ∑ x (n)e− jωn .
n =0
The DFT values can be considered as the frequency domain samples of the
Fourier transform of discrete-time signals, taken at ∆ω = 2π/N. There are
N frequency samples within the period −π ≤ ω < π,
'
'
X (k ) = X (e j2πk/N ) = X (e jω )' . (3.2)
ω =k∆ω =2πk/N
107
108 Discrete Fourier Transform
x(n)
0 N0 n
xp(n)
0 N n
"T
1
Xk = x p (t)e− j2πkt/T dt.
T
0
Assuming that the sampling theorem is satisfied, the integral can be re-
placed by a sum (in the sense of Example 2.13)
N −1
1
Xk = ∑ x (n∆t)e− j2πkn∆t/T ∆t
T n =0
Ljubiša Stanković Digital Signal Processing 109
with x p (t) = x (t) within 0 ≤ t < T. Using T/∆t = N, x (n∆t)∆t = x (n) and
X (k ) = TXk this sum can be written as
N −1
X (k ) = ∑ x (n)e− j2πkn/N . (3.3)
n =0
Therefore, the relation between the DFT and the Fourier series coeffi-
cients is
X (k ) = TXk . (3.4)
Sampling the Fourier transform of a discrete-time signal corresponds to
the periodical extension of the original discrete-time signal in time by the
period N. The period N in time is equal to the number of samples of the
Fourier transform within one period in frequency. We can conclude that this
periodic extension in time (discretization in frequency) will not influence
the possibility to recover the original signal if the original discrete-time
signal duration was not longer than N (the number of samples in the Fourier
transform of discrete-time signal).
The inverse DFT is obtained by multiplying both sides of the DFT
definition (3.1) by e j2πkm/N and summing over k
N −1 N −1 N −1
∑ X (k )e j2πmk/N = ∑ x (n) ∑ e j2πk(m−n)/N
k =0 n =0 k =0
with
N −1
1 − e j2π (m−n)
∑ e j2πk(m−n)/N = = Nδ(m − n),
k =0 1 − e j2π (m−n)/N
for 0 ≤ m, n ≤ N − 1. The inverse discrete Fourier transform (IDFT) of signal
x (n) is
1 N −1
X (k )e j2πnk/N .
N k∑
x (n) = (3.5)
=0
for 0 ≤ n ≤ N − 1.
The signal calculated by using the IDFT is, by definition, periodic with
the period N since
N −1
1
x (n + N ) = ∑ X (k )e j2π (n+ N )k/N = x (n).
N k =0
Therefore the DFT of a signal x (n) calculated using the signal samples
within 0 ≤ n ≤ N − 1 assumes that the signal x (n) is periodically extended
110 Discrete Fourier Transform
with period N as
∞
IDFT{DFT{ x (n)}} = ∑ x (n + mN )
m=−∞
∞
with ∑ x (n + mN ) = x (n) for 0 ≤ n ≤ N − 1.
m=−∞
The values of this periodical extension within the basic period are equal to
x (n). This is a circular extension of signal x (n). The following notations are
also used for this kind of the signal x (n) extension
assuming that the initial DFT was calculated for signal samples x (n) within
0 ≤ n ≤ N − 1.
In literature it is quite common to use the same notation for both x (n)
and IDFT{DFT{ x (n)}} having in mind that any DFT calculation with N
signal samples implicitly assumes a periodic extension of the original signal
x (n) with period N. Thus, we will use this kind of notation, except in the
cases when we want to emphasize a difference in the results when the
inherent periodicity in the signal (when the DFT is used) is not properly
taken into account.
Example 3.1. For the signals x (n) = 2 cos(2πn/8) for 0 ≤ n ≤ 7 and x (n) =
2 cos(2πn/16) for 0 ≤ n ≤ 7 plot the periodic signals IDFT {DFT{ x (n)}} with
N = 8 without calculating the DFTs.
∞
IDFT{DFT{ x (n)}} = ∑ x (n + 8m)
m=−∞
Example 3.3. For a signal x (n) whose values are x (0) = 1, x (1) = 1/2, x (2) = −1,
and x (3) = 1/2 find the DFT with N = 4. What is the IDFT for n = −2?
Ljubiša Stanković Digital Signal Processing 111
x(n) x(n)
0 N=8 n 0 N=8 n
...x(n-N)+x(n)+x(n+N)+.. ...x(n-N)+x(n)+x(n+N)+..
0 N=8 n 0 N=8 n
Figure 3.2 Signals x (n) = 2 cos(2πn/8) for 0 ≤ n ≤ 7 (left) and x (n) = 2 cos(2πn/16) for
0 ≤ n ≤ 7 (right) along with their periodic extensions IDFT {DFT{ x (n)}} with N = 8.
The IDFT is
1 3
[1 + cos(2πk/4) + (−1)k+1 ]e j2πnk/4 ,
4 k∑
x (n) =
=0
for 0 ≤ n ≤ 3. The DFT and IDFT inherently assume the signal and its Fourier
transform periodicity. Thus the result for n = −2 is
1 3 k 1 3 k
x (−2) = ∑ X (k )e j2π (−2) 4 = ∑ X (k )e j2π (4−2) 4 = x (4 − 2) = x (2) = −1.
4 k =0 4 k =0
Example 3.4. Assume that there is a routine to calculate the DFT of x (n) for
0 ≤ n ≤ N − 1 as X (k ) = DFT{ x (n)} = R{ x (n)}. How to use it to calculate the
DFT of a signal x (n) whose values are given within − N/2 ≤ n ≤ N/2 − 1?
⋆A periodic extension of the signal x (n) is assumed when the DFT
is calculated. It means that in the DFT calculation the signal x (n), defined
within − N/2 ≤ n ≤ N/2 − 1, will be extended with the period N. Here, we
112 Discrete Fourier Transform
Here, we have used the property that for a signal y(n) periodic with a period
N holds ∑nN=−01 y(n) = ∑nM=+MN −1 y(n) for any M (Generalize the result for the
DFT calculation and inversion for a signal x (n) defined within M ≤ n ≤
M + N − 1, using the given routine R{ x (n)}).
where
k
WN = e− j2πk/N
is used to simplify the notation, especially in graphical illustrations.
The number of additions to calculate a DFT is N − 1 for each X (k )
in (3.1). Since there are N DFT coefficients the total number of additions is
Ljubiša Stanković Digital Signal Processing 113
N ( N − 1). From the matrix from (3.6) we can see that the multiplications
are not needed for calculation of X (0). There is non need for multiplication
in the first term of each coefficient calculation as well. If we neglect the fact
that some other terms in matrix (3.6) may also assume values 1, −1, j, or − j
then the number of multiplications is ( N − 1)2 . The order of the number of
multiplications and the number of additions for the DFT calculation is N 2 .
The inverse DFT in a matrix form is
x = W−1 X, (3.9)
Most of the DFT properties can be derived in the same way as in the Fourier
transform and the Fourier transform of discrete-time signals.
N −1
1 2π
IDFT{ X (k )e− j2πkn0 /N } = ∑ X (k )e− j2πkn0 /N e j N kn
N k =0
N −1
1 2π
= ∑ X ( k ) e j N k ( n − n0 ) = x ( n − n 0 ). (3.10)
N k =0
x ∗ ( n ) = x ( N − n ).
X (k ) = X ∗ (k )
or
N −1 N −1 N −1
∑ x (n)e− j2πnk/N = ∑ x ∗ (n)e j2πnk/N = ∑ x ∗ ( N − n)e j2π ( N −n)k/N ,
n =0 n =0 n =0
X ∗ ( k ) = X ( N − k ).
Calculate the convolution x (n) ∗ x (n). Extend signals with period N = 7 and
calculate the circular convolution (corresponding to the DFT based convolu-
tion calculation with N = 7, which is longer than the signal duration). Com-
pare the results. What value of N should be used for the period so that the
direct convolution corresponds to one period of the circular convolution?
⋆Signal x (n) and its reversed version x (−n), along with the
shifted signal used in the convolution calculation, are presented in Fig.3.3.
In the circular (DFT) calculation, for example, at n = 0, the con-
volution value is
6
x p (n) ∗ x p (n) = ∑ x p (m) x p (0 − m) = 1 + 1 + 1 = 3.
m =0
In addition to the term x (0) x (0) = 1 which exists in the aperiodic convolution,
two terms for m = 3 and m = 4 appeared due to the periodic extension of
116 Discrete Fourier Transform
1.5 6
x(n)* x(n)
1 4
x(n)
0.5
2
0
-0.5 0
-15 -10 -5 0 5 10 15 -15 -10 -5 0 5 10 15
n n
1.5 1.5
xp(-m+1 )
1 1
x (m)
0.5 0.5
p
0 0
-0.5 -0.5
-15 -10 -5 0 5 10 15 -15 -10 -5 0 5 10 15
n n
1.5 1.5
xp(-m+3 )
1 1
x (-m )
0.5 0.5
p
0 0
-0.5 -0.5
-15 -10 -5 0 5 10 15 -15 -10 -5 0 5 10 15
n n
1.5 6
xp(n)* xp(n)
x (-m+5 )
1 4
0.5
2
0
p
-0.5 0
-15 -10 -5 0 5 10 15 -15 -10 -5 0 5 10 15
n n
Figure 3.3 Illustration of the discrete-time signal convolution and circular convolution for
signals whose length is 5 and the circular convolution is calculated with N = 7.
the signal. They made that the circular convolution value differs from the
convolution of original aperiodic signals. The same situation occurred for
n = 1 and n = 2. For n = 3, 4, and 5 the correct result for aperiodic convolution
is obtained using circular convolution. It could be concluded that if the signal
in circular convolution were separated by at least two more zero values (if the
period N were N ≥ 9) this difference would not occur, Fig.3.4 for N = 9. Then
one period of circular convolution 0 ≤ n ≤ N − 1 would correspond to the
original aperiodic convolution.
1.5
1
xp(-m )
xp(m)
0.5
0
-0.5
-15 -10 -5 0 5 10 15 n
n
1.5 6
xp(n)* xp(n)
xp(-m+8 )
1 4
0.5
2
0
-0.5 0
-15 -10 -5 0 5 10 15 -15 -10 -5 0 5 10 15
n n
Figure 3.4 Illustration of the discrete-time signal circular convolution for signals whose
length is 5 and the circular convolution is calculated with N = 9.
Duration of the input signal x (n) may be much longer that the dura-
tion of the impulse response h(n). For example, an input signal may have
tens of thousands of samples, while the impulse response of a discrete sys-
tem duration is, for example, tens of samples, M ≫ L. A direct convolution
would be calculated (after first L − 1 output samples) as
n
y(n) = ∑ x ( m ) h ( n − m ).
m = n − L +1
For each output sample, L multiplications would be used. For a direct DFT
application in the convolution calculation we should wait until the end of
the signal and then zero-pad both the input signal and the impulse response
118 Discrete Fourier Transform
For the convolutions yk (n) = xk (n) ∗n h(n) calculation the signals xk (n)
and h(n) should be of duration N + L − 1 only. These convolutions can be
calculated after each N ≪ M input signal samples. The output sequence
yk (n) duration is N + L − 1. Since yk (n), k = 0, 1, . . . , K − 1, are calculated
with step N in time, they overlap, although the input signals xk (n) are
nonoverlapping. For two successive yk (n) and yk+1 (n) and L ≤ N, L − 1
samples within kN + N ≤ n < kN + N + L − 1 overlap. This should be taken
into account, by summing the overlapped output samples in y(n), after the
individual convolutions yk (n) = xk (n) ∗n h(n) are calculated using the DFTs,
Fig.3.5.
2π 2π
ω= k or Ω = k, for 0 ≤ k ≤ N/2 − 1, (3.15)
N N∆t
and the other part being a shifted version of the negative frequencies (in the
original aperiodic signal)
2π 2π
ω= (k − N ) or Ω = (k − N ), for N/2 ≤ k ≤ N − 1. (3.16)
N N∆t
Illustration of the frequency correspondence to the frequency index in the
DFT is given in Fig.3.6
Ljubiša Stanković Digital Signal Processing 119
x(n)
0 n
h(n)
0 n
x (n)
1
0 n
x (n)
2
0 n
x (n)
3
0 n
y1(n)
0 n
y2(n)
0 n
y (n)
3
0 n
y(n)
0 n
Figure 3.5 Illustration of the convolution calculation when the input signal duration is much
longer then the duration of the system impulse response.
120 Discrete Fourier Transform
X(Ω)|Ω=2πk/(NΔt)
-N/2 0 N/2-1 k
Ω=2πk/(NΔt)
X(k)
0 N k
Figure 3.6 Relation between the frequency in continuous-time and the DFT frequency index.
holds in the case when the sampling theorem is satisfied, then we see that by
increasing N in the DFT calculation, the density of sampling (interpolation)
in the Fourier transform of the original signal increases. The DFT interpo-
lation by zero padding the signal in the time domain is illustrated in Fig.
3.7.
The same holds for the frequency domain. If we calculate DFT with N
samples and then add, for example, N zeros after the region corresponding
to the highest frequencies, then by the IDFT of this 2N point DFT, we will
interpolate the original signal in time. All zero values in the frequency
domain should be inserted between two parts (regions) of the original DFT
corresponding to positive and negative frequencies.
Ljubiša Stanković Digital Signal Processing 121
x(n)
n
X(k)
k
x(n)
n
X(k)
k
x(n)
n
X(k)
Figure 3.7 Discrete-time signal and its DFT (top two subplots). Discrete-time signal zero-
padded and its DFT interpolated (two subplots in the middle). Zero-padding (interpolation)
factor was 2. Discrete-time signal zero-padded and its DFT interpolated (two bottom subplots).
Zero-padding (interpolation) factor was 4. According to the duality property, the same holds if
X (k) were signal in the discrete-time and x (−n) was its Fourier transform.
122 Discrete Fourier Transform
Example 3.6. The Hann(ing) window for a signal within − N/2 ≤ n ≤ N/2 − 1, is
1 2πn
w(n) = [1 + cos( )], for − N/2 ≤ n ≤ N/2 − 1. (3.18)
2 N
If the original signal values are within 0 ≤ n ≤ N − 1 then the Hann(ing)
window form is
1 2πn
w(n) = [1 − cos( )], for 0 ≤ n ≤ N − 1. (3.19)
2 N
Present the zero-padded forms of Hann(ing) windows with 2N samples.
⋆The zero-padded form of the Hann(ing) windows used for window-
ing data within the intervals − N/2 ≤ n ≤ N/2 − 1 and 0 ≤ n ≤ N − 1 are
shown in Fig.3.8. The DFTs of windows (3.18) and (3.19) are W (k ) = N [δ(k) +
δ(k − 1)/2 + δ(k + 1)/2]/2 and W (k ) = N [δ(k) − δ(k − 1)/2 − δ(k + 1)/2]/2,
respectively. After the presented zero-padding the window DFT realness
property w pz (n) = w pz (n − 2N ) is preserved (for an even N in the case
− N/2 ≤ n ≤ N/2 − 1 and for an odd N for data within 0 ≤ n ≤ N − 1).
"∞ "∞
x (t) = 1
2π X (Ω)e jΩt dΩ, X (Ω) = x (t)e− jΩt dt.
−∞ −∞
∞
x p (t) = ∑ x (t + mT )
m=−∞
∞ T/2
"
x p (t) = ∑ Xn e j2πnt/T
, Xn = 1
T x (t)e− j2πnt/T dt,
n=−∞
− T/2
1
Xn = X (Ω)|Ω=2πn/T .
T
Ljubiša Stanković Digital Signal Processing 123
w(n)
-N/2 0 N/2-1 n
wp(n)
0 N n
w (n)
p
0 N 2N n
w(n)
0 N n
wp(n)
0 N n
w (n)
p
0 2N n
Figure 3.8 Zero-padding of the Hann(ing) windows used to window data within − N/2 ≤
n ≤ N/2 − 1 and 0 ≤ n ≤ N − 1.
124 Discrete Fourier Transform
∞
2π
X (e jω ) = ∑ X (Ω + m ) .
m=−∞ ∆t |Ω=ω/∆t
The Fourier transform of the discrete-time signal is a periodic exten-
sion X (e jω ), ω = Ω∆t, of the Fourier transform X (Ω) of a continuous-
time signal. There is no overlapping (aliasing) if the width of the
Fourier transform of the original continuous-time signal is shorter
than the extension period 2π/∆t.
4. Discrete-time periodic signal (discrete Fourier transform)
∞
x p (n) = ∑ x (n + mN ) = x p (t)|t=n∆t ,
m=−∞
N −1 N −1
x p (n) = 1
N ∑ X (k )e j2πnk/N , X (k ) = ∑ x (n)e− j2πnk/N ,
k =0 n =0
x(t) X(Ω)
X(Ω)
x(t)
t Ω
)
jω
x(n)
X(e
-π π
n ω
Xn
p
- T/2 T/2
t n
j2πk/T
xp(n) = x(n) X(k) = X(e ) = TX
k
- N/2 ≤ n < N/2 - N/2 ≤ k < N/2
x (n)
X(k)
p
n k
Figure 3.9 Aperiodic continuous-time signal and its Fourier transform (first row). Discrete-
time signal and its Fourier transform (second row). Periodic continuous-time signal and its
Fourier series coefficients (third row). Periodic discrete-time signal and its discrete Fourier
transform (DFT), (fourth row).
126 Discrete Fourier Transform
differences between the DFT and inverse DFT calculation are in the sign of
the exponent and the division of the final result by N.
Here we will present an algorithm based on splitting the signal x (n),
with N samples, into two signals x (n) for 0 ≤ n ≤ N/2 − 1 and x (n) for
N/2 ≤ n ≤ N − 1, whose duration is N/2. It is assumed that N is an even
number. By definition, a DFT of a signal with N samples is
N −1
DFT N { x (n))} = X (k ) = ∑ x (n)e− j2πnk/N
n =0
N/2−1 N −1
= ∑ x (n)e− j2πnk/N + ∑ x (n)e− j2πnk/N
n =0 n= N/2
N/2−1 @ A
= ∑ x (n) + x (n + N/2)(−1)k e− j2πnk/N
n =0
with
g(n) = x (n) + x (n + N/2).
For an odd number k = 2r + 1, follows
N/2−1
DFT N/2 {h(n)} = X (2r + 1) = ∑ h(n)e− j2πnr/( N/2)
n =0
where
h(n) = ( x (n) − x (n + N/2))e− j2πn/N .
In this way, we split one DFT of N elements into two DFTs of N/2
elements. Having in mind that the direct calculation of a DFT with N
elements requires an order of N 2 operations, it means that we will reduce
the calculation complexity, since N 2 > ( N/2)2 + ( N/2)2 . An illustration of
this calculation, with N = 8, is shown in Fig. 3.10. We can continue and
split N/2 DFTs into N/4 DFTs, and so on. A complete calculation scheme is
shown in Fig. 3.11. We can conclude that in the FFT algorithms an order of
N log2 N of operations is required. Here it is assumed that log 2 N = p is an
integer, i.e., N = 2 p . This a decimation-in-frequency algorithm.
Ljubiša Stanković Digital Signal Processing 127
x(0) X(0)
x(1) X(2)
DFT
4
x(2) X(4)
x(3) X(6)
x(4) 0 X(1)
-1 W8
x(5) 1 X(3)
-1 W8 DFT
4
x(6) 2 X(5)
-1 W8
x(7) X(7)
-1 W3
8
x(0) X(0)
x(1) 0 X(4)
-1 W8
x(2) 0 X(2)
-1 W8
x(3) 2 0 X(6)
-1 W8 -1 W8
x(4) 0 X(1)
-1 W8
x(5) 1 0 X(5)
-1 W8 -1 W8
x(6) 2 0 X(3)
-1 W8 -1 W8
x(7) X(7)
-1 W3 -1 W8
2 -1 W8
0
8
Nadditions = N log2 N.
For the number of multiplications we can see that in the first stage there are
( N/2 − 1) multiplications. In the second stage there are 2 ( N/4 − 1) mul-
B be 4C( N/8 − 1) multiplications. Finally
tiplications. In the next stage would
in the last stage would be 2 p−1 2Np − 1 = N2 ( N N − 1) = 0 multiplications
p
(N = 2 or p = log2 N). The total number of multiplications, in this algo-
rithm, is
* + * + * + * +
N N N p −1 N
Nmultiplicat. = −1 +2 −1 +4 − 1 + ... + 2 −1
2 4 8 2p
N N N N N
= − 1 + − 2 + − 4 + ... + −
2 2 2 2 2
N N 1 − 2p
= p − (1 + 2 + 22 + ... + 2 p−1 ) = p −
2 2 1−2
N N
= log2 N − ( N − 1) = [log2 N − 2] + 1.
2 2
If the multiplications by j and − j were excluded the number of multiplica-
tions would be additionally reduced.
Example 3.7. Consider a signal x (n) within 0 ≤ n ≤ N − 1. Assume that N is an
even number. Show that the DFT of x (n) can be calculated as two DFTs, one
using the even samples of x (n) and the other using odd samples of x (n).
⋆By definition
N −1
X (k ) = ∑ x (n)e− j2πkn/N
n =0
N/2−1 N/2−1
= ∑ x (2m)e− j2πk2m/N + ∑ x (2m + 1)e− j2πk(2m+1)/N
m =0 m =0
N/2−1 N/2−1
= ∑ xe (m)e− j2πkm/( N/2) + e− j2πk/N ∑ xo (m)e− j2πkm/( N/2) , (3.20)
m =0 m =0
where xe (m) = x (2m) and xo (m) = x (2m + 1) are even and odd samples of the
signal, respectively. Thus, a DFT of N elements is split into two DFTs of N/2
elements. Two DFTs of N/2 elements require an order of 2 ( N/2)2 = N 2 /2
operations. It is less than N 2 . In this way, if N/2 is an even number, we can
continue and split two DFTs of N/2 elements into four DFTs of N/4 elements,
and so on. This is a decimation-in-time algorithm, Fig.3.12.
Ljubiša Stanković Digital Signal Processing 129
x(0) X(0)
x(4) X(1)
W0 -1
8
x(2) X(2)
W0 -1
8
x(6) X(3)
W0 -1 W2 -1
8 8
x(1) X(4)
W0 -1
8
x(5) X(5)
W0 -1 W1 -1
8 8
x(3) X(6)
W0 -1 W2 -1
8 8
x(7) X(7)
W0 -1 W2 -1 W3 -1
8 8 8
M −1
X (3k) = ∑ g(n)e− j2πmk/M
m =0
M −1
X (3k + 1) = ∑ r (n)e− j2πmk/M
m =0
130 Discrete Fourier Transform
D E
with r (n) = x (m) + ax (m + M ) + a2 x (m + 3M ) e− j2πm/(3M) ,
M −1
X (3k + 2) = ∑ p(n)e− j2πmk/M
m =0
D E
with p(n) = x (m) + a2 x (m + M ) + ax (m + 3M ) e− j2π2m/(3M) , where a =
e− j2π/3 . Thus, a DFT of N = 3M elements is split into three DFTs of N/3 = M
elements. Three DFTs of N/3 elements require an order of 3 ( N/3)2 = N 2 /3
operations. If, for example, M = N/3 is an even number, we can continue
and split three DFTs of N/3 elements into six DFTs of N/6 elements, and so
on.
A periodic signal x (t), with a period T, can be reconstructed if its Fourier se-
ries is with limited number of nonzero coefficients so that Xk = 0 for k > k m
corresponding to frequencies greater than Ωm = 2πk m /T. The periodic sig-
nal can be reconstructed from the samples taken at ∆t < π/Ωm = 1/(2 f m ).
The number of samples within the period is N = T/∆t.
The reconstructed signal is
N −1 t
sin[(n − ∆t )π ]
x (t) = ∑ x (n∆t)
N sin[(n − t
n =0 ∆t )π/N ]
N −1 t
sin[(n − ∆t )π ]
x (t) = ∑ x (n∆t)e j(n−t/∆t)π/N t
n =0 N sin[(n − ∆t ) π/N ]
for an even N.
Example 3.9. Samples of a signal x (t) are taken with step ∆t = 1. Obtained discrete-
time values are x (n) = [0, 2.8284, − 2, 2.8284, 0, − 2.8284, 2, − 2.8284] for 0 ≤
n ≤ N − 1 with N = 8. Assuming that the signal satisfies the sampling
theorem find its value at t = 1.5. Check the accuracy if the original signal
values were known, x (t) = 3 sin(3πt/4) + sin(πt/4).
⋆Using the reconstruction formula for an even N we get
7
sin[(n − 1.5)π ]
x (1.5) = ∑ x(n)e j(n−1.5)π/8 8 sin[(n − 1.5)π/8] = −0.2242.
n =0
Ljubiša Stanković Digital Signal Processing 131
-2
-4
0 2 4 6 8
time
This result is equal to the original signal value. Calculation is repeated with
0 ≤ t ≤ 8, with step 0.01. The reconstructed values of x (t) are presented in
Fig.3.13.
km
x (t) = ∑ Xk e j2πkt/T . (3.21)
k=−k m
km
x (n∆t) = ∑ Xk e j2πkn/N .
k =−k m
( N −1)/2
T km T
x (n∆t)∆t = ∑ Xk e j2πkn/N = N
N k=− ∑ Xk e j2πkn/N .
k m k=−( N −1)/2
132 Discrete Fourier Transform
With x (n∆t)∆t = x (n) and TXk = X (k ) this form reduces to the DFT and the
inverse DFT
( N −1)/2 N −1
1
X (k )e j2πkn/N , x (n)− j2πkn/N .
N k=−(∑ ∑
x (n) = X (k ) =
N −1)/2 n =0
N −1 N −1
2 N −1 N −1 2
1 n j2πk t
− j2πk N 1 t n
x (t) = ∑ ∑ x (n)e e T = ∑ ∑ x (n∆t)e j2πk( T − N )
T N
k=− N2−1 n=0 n=0 k=− N −1
2
N −1 n )( N −1)/2 1 − e
j2π ( Tt − N
n )N
1 − j2π ( Tt − N
=
N ∑ x (n∆t)e t n
n =0 1 − e j2π ( T − N )
N −1 π
sin[ ∆t (t − n∆t)]
= ∑ x (n∆t) π
N sin[ N∆t (t − n∆t)]
.
n =0
This is the reconstruction formula that can be used to calculate x (t) for any
t based on the signal values at x (n∆t) with ∆t < π/Ωm = 1/(2 f m ).
In a similar way the reconstruction formula for an even number of
samples N can be obtained.
The sampling theorem reconstruction formula of aperiodic signals
follows as a special case as N → ∞, since for a small argument
π π
sin[ (t − n∆t)] → (t − n∆t)
N∆t N∆t
and
∞ π
sin[ ∆t (t − n∆t)]
x (t) → ∑ x (n∆t) π .
n=−∞ ∆t (t − n∆t )
Example 3.10. For a signal x (t) whose period is T it is known that the signal has
components corresponding to the nonzero Fourier series coefficients at k1 , k2 ,
..., k K . What is the minimal number of signal samples needed to reconstruct
the signal? What condition the sampling instants and the frequencies should
satisfy for the reconstruction?
⋆The signal x (t) can be reconstructed by using the Fourier series
(1.11). In calculations, a finite number of K nonzero terms will be used,
K
x (t) = ∑ Xkm e j2πkm t/T .
m =1
Ljubiša Stanković Digital Signal Processing 133
Since there are K unknown values Xk1 , Xk2 ,...,XkK the minimal number of
equations to calculate their values is K. The equations are written for K time
instants
K
∑ Xkm e j2πkm ti /T = x (ti ), for i = 1, 2, ..., K
m =1
or
ΦX= y, X = Φ −1 y
where
⎡ ⎤
e j2πk1 t1 /T e j2πk2 t1 /T ... e j2πkK t1 /T
⎢ e j2πk1 t2 /T e j2πk2 t2 /T ... e j2πkK t2 /T ⎥
Φ=⎢
⎣ ...
⎥
⎦
... ... ...
e j2πk1 tK /T e j2πk2 tK /T ... e j2πkK tK /T
The reconstruction condition is det ∥Φ∥ ̸= 0 for selected time instants ti and
given frequency indices k i .
Assume that the signal x (t) is sampled with ∆t. The discrete-time form
of this signal is
x (n) = Ae jω0 n ∆t,
with ω0 = Ω0 ∆t.
In order to compute the DFT of this signal, we will assume a value of
N and calculate
N −1
X (k ) = ∑ Ae jω0 n e− j2πnk/N ∆t.
n =0
N −1
1 − e jω0 N e− j2πk
X (k ) = A ∑ e jω0 n e− j2πnk/N ∆t = A ∆t (3.24)
n =0 1 − e jω0 e− j2πk/N
sin( N (ω0 − 2πk/N )/2)
= Ae j(( N −1)(ω0 −2πk/N )/2) ∆t (3.25)
sin((ω0 − 2πk/N )/2)
' '
' sin( N (ω0 − 2πk/N )/2) '
| X (k)| = | A| '' ' ∆t. (3.26)
sin((ω0 − 2πk/N )/2) '
ω0 = 2πk0 /N
N −1
X (k ) = A ∑ e j2πk0 n/N e− j2πnk/N ∆t = N Aδ(k − k0 )∆t. (3.27)
n =0
X(k)
x(n)
n k
X(k)
x(n)
n k
Figure 3.14 Sinusoid x (n) = cos(8πn/64) and its DFT with N = 64 (top row) and sinusoid
x (n) = cos(8.8πn/64) and its DFT absolute value, with N = 64 (bottom row).
N −1 B 2πk C
X (k ) = ∑ w(n) Ae jω0 n e− j2πnk/N ∆t = W e j( N −ω0 ) ∆t,
n =0
1
w(n) = [1 − cos(2nπ/N )] [u(n) − u(n − N − 1)] .
2
N −1
1
X H (k ) = ∑ [1 − cos(2nπ/N )] Ae jω0 n e− j2πnk/N ∆t
n =0 2
- .
A N −1 1 j2nπ/N 1 − j2nπ/N
Ae jω0 n e− j2πnk/N ∆t
2 n∑
= 1− e − e
=0 2 2
- .
1 1 1
= X R ( k ) − X R ( k − 1) − X R ( k + 1) ,
2 2 2
1
A= [ X (k0 ) + X (k0 + 1) + X (k0 − 1)].
N∆t
Ljubiša Stanković Digital Signal Processing 137
x(n)w(n)
XH(k)
n k
x(n)w(n)
XH(k)
n k
Figure 3.15 Sinusoid x (n) = cos(8πn/64) multiplied by a Hann(ing) window and its DFT
with N = 64 (top row) and sinusoid x (n) = cos(8.8πn/64) multiplied by a Hann(ing) window
and its DFT absolute value, with N = 64 (bottom row).
3.7.2 Displacement
A relation of the maximum DFT value with the few surrounding values of
the windowed DFT is used to calculate correction, the displacement bin of
the estimated frequency. If we apply a window function w(n) in the DFT
calculation, we get
B 2πk C
X (k ) = W e j( N −ω0 ) ∆t.
! 6
k̂0 = arg max | X (k)| ,
0≤ k ≤ N −1
138 Discrete Fourier Transform
by
X0 = | X (k̂0 )|
and two neighboring samples by
∂X (k̂0 + d)/∂d = 0,
2π
Ω0 = (k̂0 + d) (3.32)
N∆t
2π
for 0 ≤ k̂0 ≤ N/2 − 1 and Ω0 = N∆t (( k̂ 0 + d) − N ) for N/2 ≤ k̂0 ≤ N − 1.
Example 3.11. A sinusoidal signal x (t) = A exp( jΩ0 t) is sampled with a sampling
interval ∆t = 1/128 and N0 = 64 samples are considered. Prior to the DFT
calculation, the signal is zero padded four times. The DFT maximum is
detected at k̂0 = 95. The maximum DFT value is X (95) = 0.9. Neighboring
values are X (96) = 0.7 and X (94) = 0.3. Calculate the displacement bin d and
estimate the value of Ω0 .
Ljubiša Stanković Digital Signal Processing 139
X(0) X(1)
X(Ω), X(k)
X(Ω), X(k)
X(-1)
Ω, k Ω, k
Figure 3.16 Illustration of the displacement bin correction for a true maximum position
calculation based on three neighboring values (full range – left and zoomed graph – right)
.
The DFT of signal satisfies many desirable properties. Also its calculation is
simple and efficient using the FFT algorithm. With the DFT calculation the
signal periodic extension is assumed and embedded in the discrete trans-
form. However, in the DFT case the periodic signal extension will, in gen-
eral, introduce significant signal change (corresponding to discontinuities in
continuous time) at the period ending points Fig.3.17 (first and second row).
This change will significantly worsen the DFT coefficients convergence and
increases number of coefficients in signal reconstruction. In order to reduce
this effect and to improve convergence of the signal transform coefficients
the signal could be extended in an appropriate way.
The discrete cosine transforms (DCT) and discrete sine transforms
(DST) are used to analyze real-valued discrete signals, periodically extended
to produce even or odd signal forms, respectively. However, this extension
is not straightforward for discrete-time signals. Consider a discrete-time
signal of duration N, when x (n) assumes nonzero values for 0 ≤ n ≤ N − 1.
If we try with a direct extension (using all signal values) and form a periodic
signal y(n), whose basic period is of duration 2N, as
!
x (n) for 0 ≤ n ≤ N − 1
y(n) =
x (2N − n − 1) for N ≤ n ≤ 2N − 1
the obtained signal is not even, Fig.3.17(third row). It is obvious that y(n)
does not satisfy the condition y(n) = y(−n) = y(2N − n), required for a
real-valued DFT.
The same holds for an odd extension, Fig.3.17(fourth row),
!
x (n) for 0 ≤ n ≤ N − 1
y(n) = .
− x (2N − n − 1) for N ≤ n ≤ 2N − 1
Thus we have not achieved one of our goals to have a real-valued transform
after a real-valued signal periodic extension. However from Fig.3.17(third
and fourth row) we can see that the signals y(n) are even (or odd) with
respect to the vertical line at n = −1/2. Thus, if we add zeros between each
sample of y(n) and assume that the position which was at n = −1/2 in the
initial signal is the new coordinate origin n = 0 in the new signal z(n), then
these signals will be even and odd, respectively, Fig.3.17(last two rows).
This is just one of possible extensions to make the original discrete-
time signal even (or odd). Several forms of the DCT and DST are defined
based on other ways of getting an even (odd) signal extension.
The most commonly used is the so called DCT-II or just DCT. It will
be presented here. Signal extension for this transform corresponds to the
Ljubiša Stanković Digital Signal Processing 141
x(n)
z(n)
z(n)
Figure 3.17 Illustration of a signal x (n), its periodic extension corresponding to the DFT, an
even and odd discrete-time signal extension corresponding to the DCT and DST of type II.
There are two main advantages of this transform over the standard DFT
calculation. The DCT coefficients are real-valued for a real-valued signal.
142 Discrete Fourier Transform
This transform can produce a better energy concentration than the DFT. In
order to understand why a better energy concentration can be obtained we
will compare the DCT to the standard DFT
N −1
X (k ) = ∑ x (n)e− j2πnk/N ), 0≤k≤ N−1
n =0
Only N terms of the transform are used and the DCT values are obtained.
Since the basis functions are orthogonal the inverse DCT is obtained
2π (2m+1)k
by multiplying both sides of the DCT by cos ( 4N ) and summing over
0 ≤ k ≤ N − 1,
N −1 N −1
2π (2n + 1)k 2π (2m + 1)k
∑ 2x(n) ∑ wk cos(
4N
) cos(
4N
)
n =0 k =0
N −1
2π (2m + 1)k
= ∑ wk C (k ) cos(
4N
),
k =0
Ljubiša Stanković Digital Signal Processing 143
we get
N −1
1 2π (2n + 1)k
x (n) =
N ∑ wk C (k ) cos(
4N
). (3.34)
k =0
A symmetric relation, with the same coefficients in the time and fre-
quency domain, is
N −1
2π (2n + 1)k
C (k ) = vk ∑ x (n) cos(
4N
)
n =0
N −1
2π (2n + 1)k
x (n) = ∑ vk C (k ) cos(
4N
),
k =0
√ √
where v0 = 1/N and vn = 2/N for n ̸= 1.
In a similar way the discrete sine transforms are defined. The most
common form is the DST of type II (DST-II), whose definition is
N −1
2π (2n + 1)
S(k ) = ∑ 2x(n) sin( 2N
(k + 1))
n =0
z(2n + 1) = y(n)
z(2n) = 0.
X(k)
x(n)
n k
[x(n) x(n)]
X (k)
2
n k
C(k)
y(n)
n k
with N terms of the transform being used. The DST is the imaginary part of
this DFT.
Example 3.12. Consider a signal
x (n) = cos(2π (2n + 1)/64) + 0.75 cos(7π (2n + 1)/64).
Calculate its DFT with N = 32. Plot the periodic extension of this signal. Plot
the even extension y(n) of x (n). Calculate the DFT (the DCT) of such a signal
and discuss the results.
⋆Signal x (n), along with its extended versions and corresponding
transforms, is presented in Fig.3.18. Better energy concentration in the DCT
is due to the introduced symmetry in y(n). The artificial discontinuity in the
DFT, which causes its slow convergence, is eliminated in the DCT.
Here we will present two discrete signal transform that can be calculated
without using multiplications. One of them will be used to explain the basic
principle of the wavelet transform calculation as well.
Let us consider a two-sample signal x (n), with N = 2. The correspond-
ing two-sample DFT is
1
X (k) = ∑ x(n)e− j2πnk/2 = x(0) + (−1)k x(1).
n =0
It can be calculated without using multiplications, X (0) = x (0) + x (1) and
X (1) = x (0) − x (1). Now we can show that it is possible to define basis
functions for any signal duration in such a way that the multiplications
are not used in the signal transformation. These transform values will be
denoted by H (k ). For two-sample signal case
X(k)
x(n)
n k
[x(n) x(n)]
X2(k)
n k
y (n)
C(k)
c
n k
y (n)
S(k)
s
n k
Figure 3.19 Signal and its periodic extensions, corresponding to: the DFT (second row), the
cosine transform (third row), and the sine transform (fourth row). Positive frequencies for the
DFT are shown.
Example 3.14. For the signal shown in Fig. 3.20 calculate the two-sample The DFT
for each pair of signal samples
x(n)
n
y (n)
L
n
y (n)
H
Figure 3.20 Original signal x (n) and its two-sample lowpas part y L (n) and highpass part
y H ( n ).
for 0 ≤ n ≤ N/2 − 1.
In some cases the smoothed version y L (n), with a half of the samples of
the original signal, (3.20), is quite good representative of the original signal,
so there is no need to use corrections. Note that for many instants correction
is zero as well.
There are two possibilities to continue and apply the two-point DFT
scheme to a signal with N samples. One of them is in further splitting of
both y L (n) and y H (n) into their low and highpass parts. It leads to a discrete
Walsh-Hadamard transform Fig.3.21. In the other case the splitting is done
for the lowpass part only, while the highpass correction is kept as it is.
It leads to the Haar wavelet transform, Fig.3.22. These two forms will be
explained in details next.
148 Discrete Fourier Transform
x(n)
yH(n)
y (n)
L
n n
yHL(n)
y (n)
LL
n n
yHH(n)
y (n)
LH
n n
Figure 3.21 Illustration of the procedure leading to the Walsh-Hadamard transform calcula-
tion.
Let us continue the idea of splitting both (lowpass and highpass) parts of the
signal and define a transformation of a four-sample signal. For this signal
form two auxiliary two-sample signals y L (n) and y H (n) as
y L (0 ) = x (0 ) + x (1 ), y L (1 ) = x (2 ) + x (3 ) (3.37)
y H (0 ) = x (0 ) − x (1 ), y H (1 ) = x (2 ) − x (3 ). (3.38)
yH(n)
x(n)
n n
yLH(n)
y (n)
L
n n
y (n)
LL
Figure 3.22 Illustration of the procedure leading to the Haar wavelet transform calculation.
H (0 ) = y L (0 ) + y L (1 ) = x (0 ) + x (1 ) + x (2 ) + x (3 ).
H (1 ) = y L (0 ) − y L (1 ) = x (0 ) + x (1 ) − x (2 ) − x (3 ).
H (3 ) = y H (0 ) + y H (1 ) = x (0 ) − x (1 ) + x (2 ) − x (3 ).
H (4 ) = y H (0 ) − y H (1 ) = x (0 ) − x (1 ) − x (2 ) + x (3 ).
150 Discrete Fourier Transform
By replacing the values of y L (n) and y H (n) with signal values x (n),
we get the transformation equation
⎡ ⎤ ⎡ ⎤⎡ ⎤ ⎡ ⎤
H (0) 1 1 1 1 x (0 ) x (0 )
⎢ H (1) ⎥ ⎢ 1 −1 −1 ⎥ ⎢ ⎥ ⎢ ⎥
⎢ ⎥=⎢ 1 ⎥ ⎢ x (1 ) ⎥ = T 4 ⎢ x (1 ) ⎥ , (3.40)
⎣ H (2) ⎦ ⎣ 1 −1 1 −1 ⎦ ⎣ x (2 ) ⎦ ⎣ x (2 ) ⎦
H (3) 1 −1 −1 1 x (3 ) x (3 )
It would correspond to [ H (0), H (4), H (2), H (6), H (1), H (5), H (3), H (7)]T
order of coefficients in the Walsh transform with dyadic ordering (3.41).
Recursive construction of a Hadamard transform matrix H2N is easy
using the Kronecker product of T2 defined by (3.36) and HN ,
- .
HN HN
H2N = T2 ⊗ HN = .
HN −HN
Order [ H (0), H (1), H (3), H (2), H (6), H (7), H (5), H (4)] T in (3.41) would
correspond to a Walsh transform with sequency ordering.
Calculation of the Walsh-Hadamard transforms requires only addi-
tions. For an N-order transform the number of additions is ( N − 1) N.
152 Discrete Fourier Transform
Consider again two pairs of signal samples, x (0), x (1) and x (2), x (3). The
high frequency parts of these pairs are calculated as y H (n) = x (2n) − x (2n +
1), for n = 0, 1. They are used in the Haar transform without any further
modification. Since they represent highpass Haar transform coefficients
they will be denoted, in this case, by W (2) = y H (0) = x (0) − x (1) and
W (3) = y H (1) = x (2) − x (3). The lowpass coefficients of these pairs are
y L (0) = x (0) + x (1) and y L (1) = x (2) + x (3). The highpass and lowpass
parts of these signals are calculated as y LH (0) = [ x (0) + x (1)] − [ x (2) + x (3)]
and y LL (0) = [ x (0) + x (1)] + [ x (2) + x (3)]. For a four-sample signal the
transformation ends here with W (1) = y LH (0) and W (0) = y LL (0). Note
that the order of coefficients is such that the lowest frequency coefficient
corresponds to the transform index k = 0. Matrix form for a four-sample
signal is
⎡ ⎤ ⎡ ⎤⎡ ⎤
W (0) 1 1 1 1 x (0)
⎢ W (1) ⎥ ⎢ 1 1 −1 −1 ⎥ ⎢ ⎥
⎢ ⎥ ⎢ ⎥ ⎢ x (1) ⎥ .
⎣ W (2) ⎦ = ⎣ 1 −1 0 0 ⎦ ⎣ x (2) ⎦
W (3) 0 0 1 −1 x (3)
In specific, for the highest N/2 coefficients the Haar transform does
only one addition (of two signal values) for each coefficient. For next N/4
coefficients the Haar wavelet uses 4 signal values with 3 additions and so
on. The total number of additions is for a Haar transform is
N N N N
Nadditions = (2 − 1) + (4 − 1) + (8 − 1) + ... + ( N − 1).
2 4 8 N
For N of the form N = 2m we can write
1 1 1 1
Nadditions = N log2 N − N ( + 2 + 3 + ... + m )
2 2 2 2
1 1 − 21m
= N log2 N − N = N log2 N − ( N − 1) = N [log2 N − 1] + 1.
2 1 − 12
Calculate its Haar and Walsh-Hadamard transform with N = 16. Discuss the
results.
⋆Signal x (n) is presented in Fig.3.23. In full analogy with (3.43) a
Haar transformation matrix of order N = 16 is formed. For example, higher
coefficients are just two-sample signal transforms,
Although there are some short duration pulses (x (2), x (3), x (13)), the Haar
transform coefficients W (2), W (3), ..., W (8), W (10), W (11), W (12), W (13),
W (15) are zero-valued, Fig.3.23. This is the result of its property to decom-
pose the high frequency signal region into short duration (two-sample) basis
functions. Then a short duration pulse is contained in high frequency part
of only one Haar coefficient. That is not the case in the Fourier transform
(or Walsh-Hadamard transform) where a single delta pulse will cause that
all coefficients are nonzero, Fig.3.24. Transformation matrix T16 is obtained
from T8 using (3.42).
Property that high-frequency coefficients are well localized in time and
they represent a short duration signal components is used in image compres-
sion where adding high frequency coefficients adds details into an image,
with important property that one detail in the image corresponds to one (a
few) nonzero coefficient. Reconstruction with the Haar transform with dif-
ferent number of coefficients is presented in Fig.3.23. As explained it can be
considered as "a zooming" a signal toward the details when the higher fre-
quency coefficients are added. Since a half of the coefficients are zero-valued
154 Discrete Fourier Transform
W(k)
x(n)
n k
x0-1 (n)
x (n)
0
n n
x0-1,9,14 (n)
(n)
0-1,9
x
n k
Figure 3.23 Signal x (n) and its discrete Haar transform H (k ). Reconstructed signals: using
H (0) presented by x0 (n), using two coefficients H (0) and H (1) denoted by x0−1 (n), using
H (0), H (1), and H (9) denoted by x0−1,9 (n), and using H (0), H (1), H (9), and H (14) denoted
by x0−1,9,14 (n). Vertical axes scales for the signal and transform are different.
H(k)
x(n)
n k
x(n)
x(n)
n n
W(k)
W(k)
k k
H(k)
H(k)
k k
Figure 3.25 The Haar wavelet transform (second row) and the Walsh-Hadamard transform
(third row) for high frequency long duration signals (first row). Vertical axes scales for the
signal and transform are different.
3.10 PROBLEMS
Problem 3.1. Calculate the DFT of signals using the smallest possible value
of N: a) x (n) = δ(n), b) x (n) = δ(n) + δ(n − 1) − 2jδ(n − 2) + 2jδ(n − 3) +
δ(n − 4), and c) x (n) = an (u(n) − u(n − 10)).
Problem 3.2. If the signals g(n) and f (n) are real-valued show that their
DFTs, G (k ) and F (k ), can be obtained from the DFT Y (k ) of the signal
y(n) = g(n) + jh(n).
Problem 3.3. The relationship between the DFT index and the continuous
signal frequency is given by
%
2πk/( N∆t) for 0 ≤ k ≤ N/2 − 1
Ω=
2π (k − N )/( N∆t) for N/2 ≤ k ≤ N − 1.
Problem 3.7. Find the signal whose DFT is Y (k ) = | X (k )|2 and X (k ) is the
DFT of x (n) = u(n) − u(n − 3) with period N = 10.
Problem 3.8. What is the relation between the discrete Hartley transform
(DHT) of real-valued signals
N −1 * +
2πnk 2πnk
H (k ) = ∑ x (n) cos + sin
n =0
N N
Ljubiša Stanković Digital Signal Processing 157
and the DFT? Express the DHT in terms of the DFT and the DFT in terms of
the DHT.
Problem 3.9. Show that the DCT of a signal x (n) with N samples, defined
by
N −1
2πk 1
C (k ) = ∑ 2x(n) cos( 2N (n + 2 ))
n =0
as
πk
N −1 2πk n πk
C (k ) = Re{e− j 2N ∑ y(n)e− j N } = Re{e− j 2N DFT{y(n)}}.
n =0
z(2n + 1) = y(n)
z(2n) = 0.
(a) What are the real and imaginary parts of Z (k ) = DFT{z(n)}? How
they are related to the DCT and DST of x (n)? (b) The signal x (n) is applied
as an input to a system with impulse response h(n) such that h(n) is of
duration shorter than N, defined within 0 ≤ n ≤ N − 1, and x (n) ∗n h(n) is
also within 0 ≤ n ≤ N − 1. The DCT of the output signal is calculated. How
it is related to the DCT and DST of x (n)?
N −1
yk (n + ( N − 1)) = ∑ x (n + m)e− j2πmk/N
m =0
158 Discrete Fourier Transform
so that its value yk ( N − 1) at the last instant of the signal duration is equal
to the DFT of signal, for a given k,
N −1
y k ( N − 1) = ∑ x (m)e− j2πmk/N = DFT{ x (n)} = X (k ).
m =0
Note that the system is causal since yk (n) uses only x (n) at instant n and
previous instants.
Show that the output signal yk (n) is related to previous output value
yk (n − 1) by the equation
3.11 SOLUTIONS
Solution 3.1. The DFT assumes that the signals are periodic. In order to
calculate the DFT we have to assume a period of signals first. Period N
should be greater or equal to the duration of signal, so that the signal values
do not overlap. Larger values of N will increase the density of the frequency
domain samples, but will also increase the computation time.
a) For this signal any N ≥ 1 is acceptable, producing
X (k ) = 1, k = 0, 1, ..., N − 1,
with period N.
b) We may use any N ≥ 5. Using N = 5 we get:
5−1
X (k ) = ∑ x(n)e− j2πnk/5 = 1 + e− j2πk/5 − 2je− j4πk/5 + j2e− j6πk/5 + e− j8πk/5
n =0
= 1 + 2 cos(2πk/5) − 4 sin(4πk/5).
Ljubiša Stanković Digital Signal Processing 159
c) For a period N ≥ 10
9
1 − a10 e− j2πk(10/N )
X (k ) = ∑ (ae− j2πk/N )n = 1 − ae− j2πk/N
.
n =0
Solution 3.2. From y(n) = g(n) + j f (n) the real and imaginary parts g(n)
and f (n) can be obtained as
DFT{y∗ (n)} = Y ∗ ( N − k ).
Y (k ) + Y ∗ ( N − k ) Y (k ) − Y ∗ ( N − k )
G (k) = and F (k ) = .
2 2j
N −1
X1 ( k ) = ∑ x (n)(−1)n e− j2πnk/N .
n =0
For 0 ≤ k ≤ N/2 − 1
N −1 N −1
N
X1 ( k ) = ∑ x (n)e− jπn e− j2πnk/N = ∑ x (n)e− j2πn(k+ N/2)/N = X (k + ).
n =0 n =0 2
For N/2 ≤ k ≤ N − 1
N −1 N −1
N
X1 ( k ) = ∑ x (n)e jπn e− j2πnk/N = ∑ x (n)e− j2πn(k− N/2)/N = X (k − ).
n =0 n =0
2
160 Discrete Fourier Transform
N −1 N −1
Y (k ) = ∑ y(n)e− j2πnk/N = ∑ [x(n) + (−1)n x(n)]e− j2πnk/N
n =0 n =0
N −1
N
= ∑ [x(n) + x(n)e− jπnN/N ]e− j2πnk/N = X (k) + X (k + 2
)
n =0
N −1 N −1
Z (k ) = ∑ z(n)e− j2πnk/N = ∑ [x(n) − (−1)n x(n)]e− j2πnk/N
n =0 n =0
N
= X ( k ) − X ( k + ).
2
Obviously Y (k ) + Z (k ) = X (k ).
Solution 3.5. For the convolution calculation, using the DFT, the minimal
number N is N = K + L − 1 = 4, where K = 2 is the duration of x (n) and
L = 3 is the duration of h(n). With N = 4 follows
X (k ) = 1 − e− j2πk/4
H (k ) = 2 − e− j2πk/4 + 2e− j4πk/4
Y (k ) = X (k ) H (k ) = (1 − e− j2πk/4 )(2 − e− j2πk/4 + 2e− j4πk/4 )
= 2 − 3e− j2πk/4 + 3e− j4πk/4 − 2e− j6πk/4 .
The signal is
Solution 3.6. The circular convolution of y(n) = x (n) ∗ h(n) has the DFT
Y (k ) = X (k) H (k ) with
N −1
1 1
X (k ) = ∑ [e j4πn/N + 2j e j2πn/N − 2j e− j2πn/N ]e− j2πnk/N
n =0
N N
= Nδ(k − 2) + δ ( k − 1) − δ ( k + 1)
2j 2j
Ljubiša Stanković Digital Signal Processing 161
and
N −1
1 1
H (k) = ∑ [ 2 e j4πn/N + 2 e− j4πn/N + e j2πn/N ]e− j2πnk/N
n =0
N N
= δ(k − 2) + δ(k + 2) + Nδ(k − 1).
2 2
The value of Y (k ) is
N2 N2
Y (k ) = δ ( k − 2) + δ ( k − 1 ).
2 2j
Since
( )∗
N −1 N −1
IDFT{ X ∗ (k )} = ∑ X ∗ (k )e j2πnk/N = ∑ X (k )e− j2πnk/N
k =0 k =0
( )∗
N −1
= ∑ X (k )e j2πk( N −n)/N = x∗ ( N − n)
k =0
we get
Thus,
N −1
2πnk X (k ) + X ( N − k ) H (k ) + H ( N − k )
∑ x (n) cos
N
=
2
=
2
n =0
N −1
2πnk X ( N − k ) − X (k ) H (k ) − H ( N − k )
∑ x (n) sin
N
=
2j
=
2
.
n =0
2H (k ) = X (k ) + X ( N − k ) − j[ X ( N − k ) − X (k )].
2X (k ) = H (k ) + H ( N − k ) − j[ H (k ) − H ( N − k )].
Solution 3.9. We can split the DCT sum into an even and odd part
N −1
2πk 1
C (k ) = ∑ 2x(n) cos( 2N (n + 2 )) =
n =0
N/2−1 N/2−1
2πk 1 2πk 1
∑ 2x (2n) cos(
2N
(2n + )) + ∑ 2x (2n + 1) cos(
2 2N
(2n + 1 + )).
2
n =0 n =0
N/2−1
2πk 1
∑ 2x (2n + 1) cos(
2N
(2n + 1 + ))
2
n =0
N/2−1
2πk 1
= ∑ 2x ( N − 2m − 1) cos(
2N
( N − 2m − 1 + )).
2
m =0
Shifting now the summation index in this sum for N/2 + m = n follows
N/2−1
2πk 1
∑ 2x ( N − 2m − 1) cos(
2N
( N − 2m − 1 + ))
2
m =0
N −1
2πk 1
= ∑ 2x (2N − 2n − 1) cos(
2N
(2N − 2n − )).
2
n= N/2
Ljubiša Stanković Digital Signal Processing 163
Now we can go back to the DCT and to replace the second sum, to get
N/2−1
2πk 1
C (k ) = ∑ 2x (2n) cos(
2N
(2n + ))
2
n =0
N −1 N −1
2πk 1 2πk 1
+ ∑ 2x (2N − 2n − 1) cos(
2N
(2n + )) = ∑ y(n) cos(
2 2N
(2n + ))
2
n= N/2 n =0
and
N −1
Z (k ) = DFT{z(n)} = e− j2πk/(4N ) ∑ 2x(n)e− j2πnk/(2N )
n =0
jπk/(2N )
Z (k )e = Y (k) = 2X (k/2).
b) If the signal x (n) is input to a system then the DCT is calculated for
It has been assumed that all x (n), h(n), and x (n) ∗n h(n) are zero-valued
outside 0 ≤ n ≤ N − 1 (it means that the duration of x (n) and h(n) should
be such that their convolution is within 0 ≤ n ≤ N − 1) . Then for a signal
zh (n) related to xh (n) = x (n) ∗n h(n) in the same way as z(n) to x (n) in a)
we can write
k k k k
DFT{zh (n)}e jπk/(2N ) = 2Xh ( ) = 2X ( ) H ( ) = Y (k ) H ( ).
2 2 2 2
Then
k
Ch (k ) = DCT{ xh (n)} = Re{Y (k ) H ( )e− jπk/(2N ) }
2
− jπk/(4N ) k k
= Re{Y (k)e } Re{ H ( )} − Im{Y (k)e− jπk/(4N ) } Im{ H ( )}
2 2
k k
= C (k) Re{ H ( )} + S(k) Im{ H ( )}.
2 2
The system output is x (n) ∗n h(n) = xh (n) = IDCT{Ch )k )}, (3.34). Transform
H (k/2) is the DFT of zero-padded h(n) with factor 2. Only first half of the
DFT samples are then used.
Solution 3.11. For the signal yk (n) we may write
N −1
yk (n) = ∑ x (n − N + 1 + m)e− j2πmk/N .
m =0
For 0 ≤ n ≤ N − 1
Therefore H (2r ) is
N/2−1 - .
2πrn 2πrn
H (2r ) = ∑ g(n) cos
N/2
+ sin
N/2
n =0
where g(n) = x (n) + x (n + N/2). This is a DHT of g(n) with N/2 samples.
Note: For odd frequency indices k = 2r + 1 we can write
N −1 - .
2π (2r + 1)n 2π (2r + 1)n
H (2r + 1) = ∑ x (n) cos + sin .
n =0
N N
N/2−1 - .
2πnr 2πnr
H (2r + 1) = ∑ f (n) cos
N/2
+ sin
N/2
n =0
where
N 2πn N 2πn
f (n) = [ x (n) − x (n + )] cos + [ x ( − n) − x ( N − n)] sin .
2 N 2 N
This is again a DHT of a signal f (n) with N/2 samples.
166 Discrete Fourier Transform
where |X| is the vector whose elements are the DFT values | X (k )|, k =
0, 1, ..., 15. Maximal value is at k = 3, meaning that the frequency estimation
without displacement √ bin would be (2π · 3)/16 = 1.1781, while the true
frequency is (2π · 2 3)/16 = 1.3603. The error is 13.4%.
For the zero-padded signal (interpolated DFT), with a factor of 4,
15 √ 15 √
3n/16 − j2πnk/64
| X (k)| = | ∑ e j4π e | =| ∑ e j2π (8 3−k)n/64
|
n =0 n =0
' B C'
' sin π (8√3 − k )/4 '
' '
= '' B √ C '' .
' sin π (8 3 − k )/64 '
@ √ A @ √ A
Maximal value is obtained for k = 8 3 = 14, where 8 3 denotes the
nearest integer value. Then
' B C'
' sin π (8√3 − 14)/4 '
' '
| X (14)| = '' B √ C '' = 15.9662,
' sin π (8 3 − 14)/64 '
' B C'
' sin π (8√3 − 15)/4 '
' '
| X (15)| = '' B √ C '' = 13.9412
' sin π (8 3 − 15)/64 '
' B C'
' sin π (8√3 − 13)/4 '
' '
| X (13)| = '' B √ C '' = 14.8249,
' sin π (8 3 − 13)/64 '
| X (15)| − | X (13)|
d = 0.5 = −0.1395.
2 | X (14)| − | X (15)| − | X (15)|
√
The true frequency index would be 8 3 = 13.8564, with the true frequency
2π · 13.8564/64 = 1.3603. The correct value of frequency index is shifted
Ljubiša Stanković Digital Signal Processing 167
4
3
2
x( n)
1
0
-1
-2
-3
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
n
from the nearest integer k = 14 (on the frequency grid) for 14 − 13.8564 =
−0.1436, when the interpolation is done. Thus, the obtained displacement
bin value −0.1395 is close to the true shift value −0.1436. The estimated
frequency, using the displacement bin, is 1.3608. As compared to the true
frequency the error is 0.03%.
If the displacement formula is applied on √ the DFT values, without
interpolation, we would get d = 0.3356, while 2 3 = 3.4641 is displaced
from the nearest integer for 0.4641.
3.12 EXERCISE
Exercise 3.1. Find the DFT of x (n) = δ(n) − δ(n − 3) with N = 4 and N = 8.
Exercise 3.2. Calculate the DFT of signal x (n) = sin(nπ/4) for 0 ≤ n < N
with N = 8 and N = 16.
Exercise 3.3. For a real-valued signal the DFT is calculated with N = 8 and
the following DFT values are known: X (0) = 1, X (2) = 2 − j, X (5) = j,
X (7) = 3. Find the remaining values. What are the values of x (0) and
∑7n=0 x (n)?
Exercise 3.4. Signal x (n) is presented in Fig. 3.26. Find X (0), X (4), and X (8),
where X (k ) is the DFT of x (n) calculated with N = 16.
Exercise 3.5. Prove that for an arbitrary real-valued signal x (n), defined for
0 ≤ n < N, where N is an even integer, the DFT value X ( N/2) is real-valued.
168 Discrete Fourier Transform
4
3
2
X( k)
1
0
-1
-2
-3
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
k
5. Find x (0).
6. Calculate ∑15
n =0 x ( n ).
7. Calculate ∑15 2
n=0 | x ( n )| .
8. Calculate ∑15 n
n=0 (−1) x (n ).
Fourier transform of discrete signals and the DFT are used for
T
HE
direct signal processing and calculations. A transform that gener-
alizes these transforms, in the same way as the Laplace transform
generalizes the Fourier transform of continuous signals, is the z-transform.
This transform provides an efficient tool for qualitative analysis and design
of the discrete systems.
∞
X (z) = ∑ x (n )z−n , (4.1)
n=−∞
x ( n ) = a n u ( n ) + b n u ( n ),
where a and b are complex numbers, | a| < |b|. Find the z-transform of this
signal and its region of convergence.
169
170 z-Transform
a a
Im{z}
Im{z}
Im{z}
b b
Re{z} Re{z} Re{z}
a a
Im{z}
Im{z}
Im{z}
b b
Re{z} Re{z} Re{z}
4.2.1 Linearity
∞
Z{ ax (n) + by(n)} = ∑ [ ax (n) + by(n)]z−n = aX (z) + bY (z)
n=−∞
with the region of convergence being at least the intersection of the regions
of convergence of X (z) and Y (z). In special cases the region can be larger
than the intersection of the regions of convergence of X (z) and Y (z) if
some poles, defining the region of convergence, cancel out in the linear
combination of transforms.
4.2.2 Time-Shift
∞ ∞
Z{ x (n − n0 )} = ∑ x ( n − n0 ) z − n = ∑ x (n)z−(n+n0 ) = X (z)z−n0 .
n=−∞ n=−∞
Example 4.3. For a causal signal x (n) = x (n)u(n) find the z-transform of x (n +
n0 )u(n), for n0 ≥ 0.
172 z-Transform
n =0 n =0
> ?
∞
= z n0 ∑ x (n)z−n − x (0) − x (1)z−1 − ... − x (n0 − 1)z−n0 +1
n =0
@ A
= zn0 X (z) − x (0) − x (1)z−1 − ... − x (n0 − 1)z−n0 +1 .
For n0 = 1 follows Z{ x (n + 1)u(n)} = zX (z) − x (0). Note that for this signal
x ( n + n0 ) u ( n ) ̸ = x ( n + n0 ) u ( n + n0 ).
4.2.4 Differentiation
d N X (z)
Z{n(n − 1)...(n − N − 1) x (n)u(n)} = (−1) N z N .
dz N
Ljubiša Stanković Digital Signal Processing 173
∞
Z{ x (n) ∗ y(n)} = Z{ ∑ x (m)y(n − m)}
m=−∞
∞ ∞ ∞ ∞
= ∑ ∑ x (m )y (n − m )z−n = ∑ ∑ x ( m ) y ( l ) z − m − l = X ( z )Y ( z )
n=−∞ m=−∞ l =−∞ m=−∞
with the region of convergence being at least the intersection of the re-
gions of convergence of X (z) and Y (z). In the case of a product of two z-
transforms it may happen that some poles are canceled out causing that the
resulting region of convergence is larger than the intersection of the individ-
ual regions of convergence.
= lim [ x ( N + 1) − x (0)].
N →∞
Thus,
lim [ x ( N + 1) − x (0)] = zX (z) − x (0) − X (z),
N →∞
identify possible regions of convergence and find the inverse z-transform for
each of them.
⋆Obviously the z-transform has the poles z1 = 1/2 and z2 = 1/3. Since
there are no poles in the region of convergence there are three possibilities
to define the region of convergence: 1) |z| > 1/2, 2) 1/3 < |z| < 1/2, and 3)
|z| < 1/3. The signals are obtained by using power series expansion for each
case.
1) For the region of convergence |z| > 1/2 the z-transform should be
written in the form
1 1
X (z) = 1
+ 1
.
1− 2z −3z(1 − 3z )
Both of these sums converge for |z| > 1/2. The resulting power series expan-
sion of X (z) is
∞
1 −n 1 ∞ 1 −n
X (z) = ∑ 2n
z −
3z n∑ n
z
n =0 =0 3
∞ ∞
1 1
= ∑ n z−n − ∑ n z−n .
n =0
2 n =1
3
2) For 1/3 < |z| < 1/2 the z-transform should be written in the form
−2z 1
X (z) = + 1
.
1 − 2z −3z(1 − 3z )
The corresponding geometric series are
∞ 0
1 1
= ∑ (2z)n = ∑ 2−n z−n for |2z| < 1 or |z| <
1 − 2z n=0 n=−∞ 2
∞ * + ' '
1 1 n ∞
1 −n '1' 1
1
= ∑ = ∑ n z for '' '' < 1 or |z| > .
1 − 3z n =0
3z n =0
3 3z 3
They converge for 1/3 < |z| < 1/2. The resulting power series expansion is
0
1 ∞ 1 −n
∑ 2− n z − n −
3z n∑
X (z) = −2z n
z
n=−∞ =0 3
−1 1 −n ∞
1
=− ∑ n
z − ∑ n z−n .
n=−∞ 2 n =1
3
The inverse z-transform for this region of convergence is
1 1
x (n) = − u(−n − 1) − n u(n − 1).
2n 3
3) For |z| < 1/3 we can write
−2z 1
X (z) = + .
1 − 2z 1 − 3z
The corresponding geometric series are
∞ 0
1 1
= ∑ (2z)n = ∑ 2−n z−n for |2z| < 1 or |z| <
1 − 2z n=0 n=−∞ 2
∞ 0
1 1
= ∑ (3z)n = ∑ 3−n z−n for |3z| < 1 or |z| < .
1 − 3z n=0 n=−∞ 3
1 1
X (z) = e a/z = 1 + ( a/z) + ( a/z)2 + ( a/z)3 + ...
2! 3!
follows
1 2 1
x (n) = δ(n) + aδ(n − 1) + a δ(n − 2) + a3 δ(n − 3)+
2! 3!
1
= an u ( n ).
n!
The series converges for any z except z = 0.
z2 + 1
X (z) =
(z − 1/2)(z2 − 3z/4 + 1/8)
find the signal x (n) if the region of convergence is |z| > 1/2.
⋆ The denominator of X (z) will be rewritten in the form
z2 + 1 z2 + 1
X (z) = =
(z − 1/2)(z − z1 )(z − z2 ) (z − 1/2)2 (z − 1/4)
where z1 = 1/2 and z2 = 1/4. Writing X (z) in the form of partial fractions
A B C
X (z) = + +
(z − 12 )2 z − 1
2 z− 1
4
or from
1 1 1 1
(z2 + 1) = A(z − ) + B(z − )(z − ) + C (z − )2 . (4.4)
4 2 4 2
For z = 1/4 we get 17/16 = C/16 or C = 17. Value of z = 1/2 gives
1 1 1
( + 1) = A ( − )
4 2 4
178 z-Transform
For the region of convergence |z| > 1/2 and a parameter | a| ≤ 1/2 holds
∞
1 1 −1 −1
= a = z (1 + az + a 2 z −2 + . . . ) = ∑ a n −1 z − n .
z−a z (1 − z ) n =1
In general the inversion is calculated by using the Cauchy relation from the
complex analysis
O
1
zm−1 dz = δ(m),
2πj
C
Ljubiša Stanković Digital Signal Processing 179
where C is any closed contour line within the region of convergence. The
complex plane origin is within the contour. By multiplying both sides of
X (z) by zm−1 , after integration along the closed contour within the region
of convergence we get
O ∞ O
1 1
zm−1 X (z)dz = ∑ 2πj zm−1 x (n)z−n dz = x (m).
2πj n=−∞
C C
where zi are the poles of zn−1 X (z) within the integration contour C that is
in the region of convergence and k is the pole order. If the signal is causal,
n ≥ 0, and all poles of zn−1 X (z) within contour C are simple (first-order
poles with k = 1) then, for a given instant n,
M N
x (n) = ∑ [zn−1 X (z)(z − zi )]|z=zi . (4.5)
zi
Y ( z ) = X ( z ) H ( z ).
It means that the z-transform exists at |z| = 1, i.e., that the circle
|z| = 1
Ljubiša Stanković Digital Signal Processing 181
a a a
Im{z}
Im{z}
Im{z}
1 1 1
c
b b
Re{z} Re{z} Re{z}
2 4
60 h (n) 1.5 h2(n) 3 h3(n)
1
40 1 2
20 0.5 1
0 0 0
Figure 4.3 Regions of convergence (gray) with corresponding signals. Poles are denoted by
"x".
| H (e jω )| = | H (z)||z=e jω .
where z0i are zeros and z pi are poles of the transfer function. For the
amplitude of frequency response we my write
' '
' B ' TO1 TO2 ...TO M
| H (e jω )| = '' 0 ''
A0 TP1 TP2 ...TPN
Example 4.9. Plot the frequency response of the causal notch filter with the transfer
function
z − e jπ/3
H (z) =
z − 0.95e jπ/3
TO1
| H (e jω )| =
TP1
O1
1.5
T P
1
ω
π/3
|H(ejω)|
Im{z}
0.5
0
Re{z} -2 0 π/3 2 ω
Figure 4.4 Poles and zeros of a first-order notch filter (left). The frequency response of this
notch filter (right).
B0 + B1 z−1 + ... + B M z− M
Y (z) = X ( z ).
1 + A1 z−1 + ... + A N z− N
184 z-Transform
5 1
y ( n ) − y ( n − 1) + y ( n − 2) = x ( n ). (4.6)
6 6
If the input signal is x (n) = 1/4n u(n) find the output signal.
1
Y (z) = X ( z ).
1 − 12 z−1 + 16 z−2
The z-transform of the input signal is X (z) = 1/(1 − 14 z−1 ) for |z| > 1/4. The
output signal z-transform is
z3
Y (z) = .
(z − 12 )(z − 13 )(z − 14 )
For a causal system the region of convergence is |z| > 1/2. The output signal
is the inverse z-transform of Y (z). For n > 0 it is
M N
y(n) = ∑ [zn−1 Y (z)(z − zi )]|z=zi
zi =1/2,1/3,1/4
z n +2 z n +2 z n +2
= 1 1
+ 1 1
+
(z − 3 )( z − 4 ) |z=1/2 (z − 2 )( z − 4 ) |z=1/3 (z − 12 )(z − 13 ) |z=1/4
1 8 3
=6 − n + n.
2n 3 4
For n = 0 there is no pole at z = 0. Thus, the above expressions hold for n = 0
as well. The output signal is
- .
6 8 3
y(n) = − + u ( n ).
2n 3n 4n
Note: This kind of solution assumes the initial values from the system causal-
ity and x (n) as y(0) = x (0) = 1 and y(1) − 5y(0)/6 = x (1), i.e., y(1) =
13/12.
Find its impulse response and discuss its behavior in terms of the system
coefficients.
Ljubiša Stanković Digital Signal Processing 185
⋆For the impulse response calculation the input signal is x (n) = δ(n)
with X (z) = 1. Then we have
The pole of this system is z = − A1 . The are two possibilities for the region
of convergence |z| > | A1 | and |z| < | A1 |. For a causal system the region of
convergence is |z| > | A1 |. Thus, the z-transform Y (z) can be expanded into a
geometric series with q = A1 z−1 = ( A1 /z) < 1
B CB C
Y (z) = B0 + B1 z−1 1 − A1 z−1 + A21 z−2 − A31 z−3 + ... + (− A1 z−1 )n + ...
∞ ∞
= B0 + B0 ∑ (− A1 )n z−n + B1 ∑ (− A1 )(n−1) z−n
n =1 n =1
with
y(n) = B0 δ(n) + (− A1 )n−1 (− A1 B0 + B1 )u(n − 1).
We can conclude that, in general, the impulse response has an infinite dura-
tion for any A1 ̸= 0. It is a result of the recursive relation between the output
y(n) and its previous value(s) y(n − 1). This kind of systems are referred to
as infinite impulse response (IIR) systems or recursive systems. If the value
of coefficient A1 is A1 = 0 then there is no recursion and
Then we have a system with a finite impulse response (FIR). This kind of
system produces an output to a signal x (n) as
They are called moving average (MA) systems. Systems without recursion
are always stable since a finite sum of finite signal values is always finite.
Systems that would contain only x (n) and the output recursions, in this
case,
y(n) + A1 y(n − 1) = B0 x (n)
are auto-regressive (AR) systems or all pole systems. This kind of systems
could be unstable, due to recursion. In our case the system is obviously
unstable if | A1 | > 1. Systems (4.7) are in general auto-regressive moving
average (ARMA) systems.
186 z-Transform
If the region of convergence were |z| < | A1 | then the function Y (z)
would be expanded into series with q = z/A1 < 1 as
* + ∞
B0 + B1 z−1 B0 B1
(− A1−1 z)n
A1 n∑
Y (z) = = z +
A1 z−1 (z/A1 + 1) A1 =0
0
B1 0
= B0 ∑ (− A1 )n−1 z−(n−1) + ∑ (− A1 )n z−n
n=−∞ A1 n=− ∞
−1 B1 0
= B0 ∑ (− A1 )n z−n + ∑ (− A1 )n z−n
n=−∞ A1 n=− ∞
with
B1
y(n) = B0 (− A1 )n u(−n − 1) + (− A1 )n u(−n).
A1
This system would be stable if |1/A1 | < 1 and unstable if |1/A1 | > 1, having
in mind that y(n) is nonzero for n < 0. This is an anticausal system since it
has impulse response satisfying h(n) = 0 for n ≥ 1.
Here, we have just introduced the notions. These systems will be
considered in Chapter 5 in details.
yi (n) = Ci λin .
y(0) = C1 + C2 + 1 = 1
C1 C
y (1) = + 2 +4=5
2 3
the constants C1 = 6 and C2 = −6 follow. The final solution is
- .
6 6
y(n) = n − n + 3n + 1 u(n).
2 3
Note: The z-transform based solution would assume y(0) = x (0) =
11/6 and y(1) = 5y(0)/6 + x (1) = 157/36. The solution with the initial
conditions y(0) = 1 and y(1) = 5 could be obtained from this solution with
appropriate changes of the first two samples of the input signal in order to
take into account the previous system state and to produce the given initial
conditions y(0) = 1 and y(1) = 5 .
If multiple polynomial roots are obtained, for example λi = λi+1 , then
yi (n) = λin and yi+1 (n) = nλin .
Consider a periodic signal x (n) with a period N and its DFT values X (k ),
1 N −1
x (n) = ∑ X (k )e j2πnk/N . (4.12)
N k =0
1 1 − e j2πk0 z− N '
'
= z N −1 X ( k 0 ) '
N 1 − e j2πk0 /N z−1 ' j2πk0 /N
z=e
1 z N − e j2πk0
= X (k0 ) lim j2πk0 /N
= X ( k 0 ).
N z→e j2πk0 /N z − e
190 z-Transform
z=ej2πk/N z=e
j2πk/N
k0=k
Im{z}
Im{z}
Im{z}
z=ej2πk0/N
k0≠ k
Re{z} Re{z} Re{z}
Figure 4.5 Zeros and the pole in Z{ xk (n)} (left), the pole in 1/ (1 − e j2πk0 n/N z−1 ) for k ̸= k0
(middle), and the pole in 1/ (1 − e j2πk0 n/N z−1 ) for k = k0 (right). Illustration is for N = 16.
y k ( N − 1) = X ( k ) δ ( k − k 0 ).
It requires just one additional complex multiplication for the last instant
and for one frequency. The total number of multiplications is 2N + 4. It
is reduced with respect to the previously needed 4N real multiplications.
The total number of additions is 4N + 2. It is increased. However the time
needed for a multiplication is much longer than the time needed for an
addition. Thus, the overall efficiency is improved. The efficiency is even more
improved having in mind that (4.18) is the same for calculation of X (k0 ) and
X (−k0 ) = X ( N − k0 ).
"∞ ∞ ∞
X (s) = x (t)e−st dt ∼
= ∑ x (n∆t)e−sn∆t ∆t = ∑ x (n)e−sn∆t
−∞ n=−∞ n=−∞
with x (n) = x (n∆t)∆t. Comparing this relation with the z-transform defini-
tion we can conclude that the Laplace transform of x (t) corresponds to the
z-transform of its samples with
z = exp(s∆t),
that is,
X (s) ↔ X (z)|z=exp(s∆t) . (4.19)
A point s = σ + jΩ from the Laplace domain maps into the point
z = re jω with r = eσ∆t and ω = Ω∆t. Points from the left half-plane in the
s domain, σ < 0, map to the interior of unit circle in the z domain, r < 1.
192 z-Transform
Example 4.14. A causal discrete-time signal x (n) has the Fourier transform X (e jω ).
Write its z-transform in terms of the Fourier transform of the discrete-time
signal, i.e., write the z-transform value based on its values on the unit circle.
⋆The signal can be expressed in term of its Fourier transform as
"π
1
x (n) = X (e jω )e jωn dω
2π
−π
∞ "π ∞
−n 1
X (z) = ∑ x (n)z =
2π
X (e jω ) ∑ e jωn z−n dω
n =0 −π n =0
"π
1 X (e jω )
= dω,
2π 1 − e jω z−1
−π
N −1
X (k ) = X (e jω )|ω =2πk/N = X (z)|z=e j2πk/N = ∑ x (n)z−n j2πk/N .
n =0
|z=e
Example 4.15. Consider a discrete-time signal with N samples different from zero
within 0 ≤ n ≤ N − 1. Show that all values of X (z), for any z, can be calculated
based on its N samples on the unit circle in the z-plane.
⋆If the signal has N nonzero samples, then it can be expressed in term
of its DFT as
N −1
1 N −1
∑ x (n)e− j2πnk/N and x (n) = X (k )e j2πnk/N .
N k∑
X (k ) =
n =0 =0
Thus, the z-transform of x (n), using only the values of the IDFT where the
original signal is nonzero, 0 ≤ n ≤ N − 1,
1 N −1 N −1 1 N −1 1 − z− N e j2πk
∑ ∑ X (k )e j2πnk/N z−n =
N k∑
X (z) = −1 j2πk/N
X (k)
N k =0 n =0 =0 1 − z e
Ljubiša Stanković Digital Signal Processing 193
N=16
jω j2π k/16
z=e z=e
π/Δt
Im{s}=Ω
Im{z}
Im{z}
0
- π/Δt
Figure 4.6 Illustration of the z-transform relation with the Laplace transform (left), the Fourier
transform of discrete signals (middle), and the DFT (right).
4.7 PROBLEMS
Problem 4.1. Find the z-transform and the region of convergence for the
following signals:
(a) x (n) = δ(n − 2),
(b) x (n) = a|n| u(n),
(c) x (n) = 21n u(n) + 31n u(n)
Problem 4.2. Find the z-transform and the region of convergence for the
following signals:
(a) x (n) = δ(n + 1) + δ(n) + δ(n − 1),
(b) x (n) = 21n [u(n) − u(n − 10)].
Problem 4.3. Using the z-transfrom property that
dX (z)
Y (z) = −z
dz
corresponds to
y(n) = nx (n)u(n)
194 z-Transform
in the discrete-time domain, with the same region of convergence for X (z)
and Y (z), find a causal signal whose z-transform is
(a) X (z) = e a/z , |z| > 0.
(b) X (z) = ln(1 + az−1 ), |z| > | a|.
Problem 4.4. (a) How the z-transform of x (−n) is related to the z-transform
of x (n)?
(b) If the signal x (n) is real-valued show that its z-transfrom satisfies
X ( z ) = X ∗ ( z ∗ ).
Problem 4.5. If X (z) is the z-transform of a signal x (n) find the z-transform
of
∞
y(n) = ∑ x ( k ) x ( n + k ).
k =−∞
z+1
X (z) = .
(2z − 1)(3z + 2)
3 − 56 z−1
H (z) = .
(1 − 14 z−1 )(1 − 13 z−1 )
Problem 4.10. Find the impulse response of a causal system whose transfer
function is
z+2
H (z) = .
( z − 2) z2
Problem 4.11. Find the inverse z-transform of
z2
X (z) = .
z2 + 1
5 1 5 3
y ( n ) − y ( n − 1) + y(n − 2) − y(n − 3) = 3x (n) − x (n − 1) + x (n − 2).
16 16 4 16
3 1
y ( n ) = x ( n ) − x ( n − 1) + x ( n − 2)
4 8
has a finite output duration for an infinite duration input x (n) = 1/4n u(n) .
Using the z-transform find the output to the input signal x (n) = u(n) −
u ( n − 6) .
11 1 3
y(n) − y(n − 1) + y(n − 2) = 2x (n) − x (n − 1)
6 2 2
x (n + 2) + 3x (n + 1) + 2x (n) = 0
with the initial condition x (0) = 0 and x (1) = 1. Signal x (n) is causal.
196 z-Transform
x ( n + 1) = x ( n ) + a n
to the input signal x (n) = 31n u(n) by a direct solution of the differential
equation in the discrete-time domain and by using the z-transform. The
initial conditions are y(n) = 0 for n < 0.
Problem 4.19. The first backward difference is defined as
∇ x ( n ) = x ( n ) − x ( n − 1 ),
∇ m x ( n ) = ∇ m −1 x ( n ) − ∇ m −1 x ( n − 1 ).
∆x (n) = x (n + 1) − x (n),
∆ m x ( n ) = ∆ m −1 x ( n + 1 ) − ∆ m −1 x ( n ).
sampled at ∆t = 1/60.
Ljubiša Stanković Digital Signal Processing 197
Problem 4.21. Plot the frequency response of the discrete system (comb
filter)
1 − z− N
H (z) =
1 − rz− N
with r = 0.9999 and r 1/N ∼
= 1. Show that this system has the same transfer
function as
4.8 SOLUTIONS
for any z ̸= 0.
(b) For this signal
∞ −1 ∞
(1 − a 2 ) z
X (z) = ∑ a|n| z−n = ∑ a−n z−n + ∑ an z−n = (1 − az)(z − a)
n=−∞ n=−∞ n =0
for |z| < 1/a and |z| > a. If | a| < 1 then the region of convergence is
a < |z| < 1/a.
(c) In this case
∞ ∞
1 −n 1 1 1
X (z) = ∑ n
z + ∑ n z−n = 1 −1
+ 1 −1
n =0 2 n =0 3 1 − 2 z 1 − 3z
2 − 56 z−1 z(2z − 56 )
X (z) = =
(1 − 12 z−1 )(1 − 13 z−1 ) (z − 12 )(z − 13 )
for |z| > 1/2 and |z| > 1/3. The region of convergence is |z| > 1/2.
Solution 4.2. (a) The z-transform is
∞
X (z) = ∑ (δ(n + 1) + δ(n) + δ(n − 1)) z−n =
n=−∞
1
= z + 1 + z −1 = z + 1 + .
z
198 z-Transform
j2π/10
z=e /2
Im{z}
z=1/2
Re{z}
The z-transform is
∞ 9 9
1 −n −n 1 − (2z)−10
X (z) = ∑ x (n )z−n = ∑ z = ∑ ( 2z ) = =
n=−∞ n =0 2
n
n =0 1 − (2z)−1
z−10 z10 − ( 12 )10 z10 − ( 12 )10
= =
z −1 z − 12 z9 (z − 12 )
The expression for X (z) is written in this way in order to find the region of
convergence, observing the zero-pole locations in the z-plane, Fig.4.7. Poles
are at z p1 = 0 and z p2 = 1/2. Zeros are z0i = e j2iπ/10 /2, Fig.4.7. Since the z-
transform has a zero at z0 = 1/2, it will cancel out the pole z p2 = 1/2. The
resulting region of convergence will include the whole z plane, except the
point at z = 0.
Solution 4.3. (a) For X (z) = e a/z holds
dX (z) a a
−z = z 2 e a/z = X (z)
dz z z
Ljubiša Stanković Digital Signal Processing 199
nx (n)u(n) = ax (n − 1)u(n)
dX (z)
since Z [nx (n)] = −z dz and z−1 X (z) = Z [ x (n − 1)]. It means that
a
x (n) = x ( n − 1)
n
It means that
a2 a3
x (1) = a, x (2) = , x (3 ) = ,...
2 2·3
or
an
x (n) = u ( n ).
n!
(b) For X (z) = ln(1 + az−1 )
Therefore
dX (z) az−1
Z [nx (n)] = −z =
dz 1 + az−1
nx (n) = a(− a)n−1 u(n − 1),
producing
−(− a)n
x (n) = u ( n − 1 ).
n
Solution 4.4. (a) The z-transform of signal x (−n) is
∞
X1 ( z ) = ∑ x (−n)z−n .
n=−∞
1
Y ( z ) = X ( z ) X ( ).
z
Solution 4.6. A direct expansion of the given transform into power series,
within the region of convergence, will be used. In order to find the signal
x (n) whose z-transform is X (z) = 2−13z , it should be written in a form of
' '
power series with respect to z−1 . Since the condition ' 3z '
2 < 1 does not
correspond to the region of convergence given in the problem formulation
we have to rewrite X (z) as
1 1
X (z) = − 2
.
3z 1 − 3z
'2'
Now the condition ' 3z ' < 1, that is |z| > 2 , corresponds to the problem for-
3
mulation region of convergence. In order to obtain the inverse z-transform,
write
1 1 1
X (z) = − 2
= − X1 ( z ) ,
3z 1 − 3z 3z
where
1
X1 ( z ) = 2
.
1 − 3z
Ljubiša Stanković Digital Signal Processing 201
∞ * +n ∞ * +n
2 2
X1 ( z ) = ∑ = ∑ z−n .
n =0 3z n =0 3
∞
X (z) = ∑ x (n ) z−n (4.21)
n=−∞
or
* + n −1
1 2
x (n) = − u ( n − 1 ).
3 3
Solution 4.7. Since the signal is causal the region of convergence is outside
the pole with the largest radius (outside the circle passing through this pole).
202 z-Transform
and
* + * +
B 1 B ∞ 2 n −n B ∞ 2 n − n −1 2
3z n∑
2
= − z = ∑ − z , |z| > .
3z 1 + 3z =0 3 3 n =0 3 3
3 − 56 z−1 A B
H (z) = 1 −1 1 −1
= 1 −1
+
(1 − 4 z )(1 − 3 z ) 1− 4z 1 − 13 z−1
with A = 1, B = 2.
(a) The region of convergence must contain |z| = 1, for a stable system. It is
|z| > 13 .
From
1 2
H (z) = + =
1 − 14 z−1 1 − 13 z−1
∞ * +n ∞ * +n
1 −n 1 1 1
= ∑ z +2 ∑ z−n , |z| > and |z| >
n =0
4 n =0
3 3 4
h ( n ) = ( 4− n + 2 × 3− n ) u ( n ).
(b) The region of convergence is 14 < |z| < 13 . The first term in H (z) is the
same as in (a), since it converges for |z| > 14 . It corresponds to the signal
4−n u(n). The second term must be rewritten in such a way that its geometric
series converges for |z| < 13 . Then
∞ −1
2 3z 1
1 −1
= −2 = −2 ∑ (3z)n = −2 ∑ (3z)−m with |z| < .
1− 3z
1 − 3z n =1
m=−n
m=−∞ 3
c) For an anticausal system the region of convergence is |z| < 14 . Now the
second term in H (z) is the same as in (b). For |z| < 14 the first term in H (z)
should be written as:
∞ −1
1 4z 1
=− = − ∑ (4z)n = − ∑ (4z)−m with |z| < .
1 − 14 z−1 1 − 4z n =1
m=−n
m=−∞ 4
204 z-Transform
The signal corresponding to this term is −4−n u(−n − 1). The impulse
response of the anticausal discrete system with given transfer function is
can be written as
1
H (z) = √ √
3 3
(1 − 4z)(z − 4 + j 14 )(z − 4 − j 14 )
√ √
with poles z1 = 1/4, z2 = 43 − j 14 , and z3 = 43 + j 14 . Since |z2 | = |z3 | = 1/2
possible regions of convergence are: 1) |z| < 1/4, 2) 1/4 < |z| < 1/2, and 3)
|z| > 1/2. In the first two cases the system is neither causal nor stable, while
in the third case the system is causal and stable since |z| = 1 and |z| → ∞
belong to the region of convergence.
The output to x (n) = 2 cos(nπ/2) = 1 + cos(nπ ) = 1 + (−1)n is y(n) =
H (e )|ω =0 × 1 + H (e jω )|ω =π × (−1)n = H (z)|z=1 + H (z)|z=−1 (−1)n =
jω
−0.8681 + 0.0945(−1)n .
Solution 4.10. The transfer function can be written as
z+2 A B C
H (z) = = + + 2.
z2 ( z − 2) z−2 z z
Az2 + Bz(z − 2) + C (z − 2) = z + 2
( A + B)z2 + (−2B + C ) − 2C = z + 2.
A+B=0
−2B + C = 1
−2C = 2,
z −1 1 1
H (z) = − 1
− 2− .
1 − 2z z z
Ljubiša Stanković Digital Signal Processing 205
The region of convergence for a causal system is |z| > 2. The inverse z-
transform for a causal system is the system impulse response
h ( n ) = 2n −1 u ( n − 1 ) − δ ( n − 2 ) − δ ( n − 1 ) = δ ( n − 2 ) + 2n −1 u ( n − 3 ).
For the region of convergence defined by |z| > 1 the signal is causal and
1 1
x (n) = [1 + (−1)n ] jn u(n) = [1 + (−1)n ]e jπn/2 u(n).
2 2
For n = 4k, where k ≥ 0 is an integer, x (n) = 1 , while for n = 4k + 2 the
signal values are x (n) = −1. For other n the signal is x (n) = 0.
For |z| < 1 the inverse z-transform is
1
x (n) = − [1 + (−1)n ] jn u(−n − 1).
2
Solution 4.12. The transfer function of this system is
3 − 54 z−1 + 3 −2
16 z 3 − 54 z−1 + 3 −2
16 z
H (z) = 5 −2 1 −3
=
1 − z −1 + 16 z − 32 z (1 − 12 z−1 + 1 −2 1 −1
16 z )(1 − 2 z )
1 1 1
= 1 −1
+B C2 + .
1 − 4z 1 − 14 z−1 (1 − 12 z−1 )
For a causal system the region of convergence is outside of the pole z = 1/2,
that is |z| > 1/2. Since
* +'
1 d z '
B C2 = '
da 1 − az − 1 '
1 − 14 z−1 a=1/4
' '
d ∞ n −(n−1) '' ∞ '
n−1 −(n−1) '
∞
1
= ∑ a z ' = ∑ na z ' = ∑ ( n + 1) n z − n ,
da n=0 ' n =0
' n =0
4
a=1/4 a=1/4
3 1
y ( n ) = x ( n ) − x ( n − 1) + x ( n − 2)
4 8
is
3 1
H ( z ) = 1 − z −1 + z −2 .
4 8
The z-transform of the input signal x (n) = 1/4n u(n) is
1
X (z) = ,
1 − 14 z−1
with the region of convergence |z| > 1/4. The output signal z-transform is
1
H (z) =
1 − 13 z−1
1 − z −6
X ( z ) = 1 + z −1 + z −2 + z −3 + z −4 + z −5 = .
1 − z −1
The z-transform of the output signal is
1 − z −6
Y (z) = = Y1 (z) − Y1 (z)z−6
(1 − z−1 )(1 − 1/3z−1 )
with
1 3/2 1/2
Y1 (z) = = − .
(1 − z−1 )(1 − 1/3z−1 ) 1 − z−1 1 − 13 z−1
Its inverse is - * +n .
3 1 1
y1 ( n ) = − u ( n ).
2 2 3
Ljubiša Stanković Digital Signal Processing 207
Im{z}
Im{z}
Im{z}
1/3 1/3
3/2 3/2 3/2
Figure 4.8 Poles and zeros of the system (left), input signal z-transform (middle), and the
z-transform of the output signal (right).
The output signal transform does not have a pole z = 3/2 since this pole is
canceled out. The output signal is
1 3 1
y(n) = u(n) − u ( n − 1).
3 n 2 3n −1
208 z-Transform
with
z 1 1
X (z) = = − .
z2 + 3z + 2 1 + z −1 1 + 2z−1
The inverse z-transform of X (z) is
Solution 4.17. The z-transforms of the left and right side of the equation are
z
zX (z) − zx (0) = X (z) +
z−a
- .
z 1 1 a
X (z) = = − .
(z − a)(z − 1) 1 − a z − 1 z − a
The inverse z-transform is
1 1 − an
x (n) = [u(n − 1) − an u(n − 1)] = u ( n − 1)
1−a 1−a
or
n −1
x (n) = ∑ ak , n > 0.
k =0
√ √
2 2
with λ1,2 = 4 ±j 4 . The homogenous solution is
√ √ √ √
2 2 n 2 2 n
yh (n) = C1 ( +j ) + C2 ( −j )
4 4 4 4
1 1
= C1 n e jnπ/4 + C2 n e− jnπ/4 .
2 2
A particular solution is of the input signal x (n) = 31n u(n) form. It is y p (n) =
A 31n u(n). The constant A is obtained by replacing this signal into (4.20)
√
1 2 1 1 1 1
A n− A + A n −2 = n
3 2 3n −1 4 3 3
√
3 2 9
A (1 − + ) = 1.
2 4
Its value is A = 0.886. The general solution is
1 jnπ/4 1 1
y(n) = yh (n) + y p (n) = C1 n
e + C2 n e− jnπ/4 + 0.886 n .
2 2 3
Since the system is causal with y(n) = 0 for n < 0 then the constants C1
and
√
C2 may be obtained from the initial condition following from
√
y(n) −
2 2
2 y(n − 1) + 14 y(n − 2) = x (n) as y(0) = x (0) = 1 and y(1) = 2 y (0 ) +
√
x (1) = 22 + 13 ,
C1 + C2 + 0.886 = 1 (4.23)
√ √ √ √ √
2 2 2 2 1 2 1
C1 ( +j )/2 + C2 ( −j )/2 + 0.886 = + ,
2 2 2 2 3 2 3
as C1 = 0.057 − j0.9967 = 0.9984 exp(− j1.5137) = C2∗ . The final solution is
1 1
y(n) = 2 × 0.9984 cos(nπ/4 − 1.5137) + 0.886 n .
2n 3
For the z-domain we write
√
2 1
Y (z) − Y ( z ) z −1 + Y ( z ) z −2 = X ( z )
2 4
with
1 1
Y (z) = √
1− 2 −1
+ 1 −2 1 − 13 z−1
2 z 4z
210 z-Transform
with
z3
Y (z) = √ √ √ √
2 2 2 2 1
(z − ( 4 +j 4 ))( z − ( 4 −j 4 ))( z − 3 )
Using, for example, the residual value based inversion of the z-transform,
M N
n −1
y (n) = ∑ [
√
z
√
Y ( z )( z − z )]
i | z = zi
2 2
z1,2,3 = 4 ± j 4 ,1/3
' '
' '
1 ' '1
= z n +2 √ √ '
' + z n +2 √ √ '
'
2− j 2 1 '√ √ 2+ j 2 1 ' √ √
(z − 4 )(z − 3 ) 2+ j 2 (z − 4 )(z − 3 ) z= 2− j 2
4 4
'
'
1 '
+ z n +2 √ √ √ √ '
'
(z − 2+4 j 2 )(z − 2−4 j 2 ) 'z=1/3
(√ √ ) n +2 (√ √ ) n +2
1 2+j 2 1 1 2−j 2 1
= √ √ √ − √ √ √
j 2 4 2+ j 2
− 1 j 2 4 2− j 2
− 1
2 4 3 2 4 3
1 1
+ √
3n +2 ( 1 − 1 2 + 1 )
9 3 2 4
√ √
1 −j 2 1 j 2 1
= e j(n+2)π/4 √ √ + e− j(n+2)π/4 √ √ + 0.886
2n +2 2+ j 2
− 13 2n +2 − 13 2− j 2 3n
4 4
√ √
1 2 1 2 1
= n e jnπ/4 √ √ 4
+ n e− jnπ/4 √ √ 4
+ 0.886 n
2 2+j 2− 3 2 2−j 2− 3 3
1 1
= 2 × 0.9984 n cos(nπ/4 − 1.5137) + 0.886 n ,
2 3
for n ≥ 1. For n = 0 there is no additional pole at z = 0 the previous result
holds for n ≥ 0.
Solution 4.19. The z-transform of the first backward difference is
Its z-transform is
Z [∇2 x (n)] = (1 − z−1 )2 X (z).
In the same way we get
m −1
Z [∆m x (n)] = ( z − 1)m X (z) − z ∑ ( z − 1 ) m − j −1 ∆ j x (0 ).
j =0
The values of TP1 and TO1 , and TP2 and TO2 , are almost the same for
any ω except ω = ±π/4 where the distance to the transfer function zero is
212 z-Transform
O1
1.5
T
P
1
|H(ejω)|
Im{z}
P
2
0.5
O
2
0
Re{z} -2 - π/4 0 π/4 2 ω
Figure 4.9 Location of zeros and poles for a second order system.
0, while the distance to the corresponding pole is small but finite. Based on
this analysis the amplitude of frequency response is presented in Fig.4.9.
The input discrete-time signal is
This system will filter out signal components at ω = ±π/4. The output
discrete-time signal is
z−
o
N
= 1 = e− j2πm
zom = e j2πm/N , m = 0, 1, ..., N − 1
Similarly, the poles are zmp = r1/N e j2πm/N , m = 0, 1, ..., N − 1. The frequency
response of the comb filter is
N −1 N −1
z − zom z − e j2πm/N
H (z) = ∏ = ∏ 1/N e j2πm/N
.
m=0 z − z pm m =0 z − r
Ljubiša Stanković Digital Signal Processing 213
| H (e jω )| ∼
= 1 for z ̸= e j2πm/N
| H (e jω )| = 0 for z = e j2πm/N .
4.9 EXERCISE
Exercise 4.1. Find the z-transform and the region of convergence for the
following signals:
(a) x (n) = δ(n − 3) − δ(n + 3),
(b) x (n) = u(n) − u(n − 20) + 3δ(n),
(c) x (n) = 1/3|n| + 1/2n u(n),
(d) x (n) = 3n u(−n) + 2−n u(n),
(e) x (n) = n(1/3)n u(n).
(f) x (n) = cos(n π2 ).
Exercise 4.2. Find the z-transform and the region of convergence for the
signals:
(a) x (n) = 3n u(n) − (−2)n u(n) + n2 u(n).
(b) x (n) = ∑nk=0 2k 3n−k ,
(c) x (n) = ∑nk=0 3k .
Exercise 4.3. Find the inverse z-transform of:
−8
(a) X (z) = 1z−z + 3, if X (z) is the z-transform of a causal signal x (n).
(b) X (z) = (zz−+22)z2 , if X (z) is the z-transform of a causal signal x (n).
2+3z−2 , if X (z ) is the z-transform of an unlimited-duration
(c) X (z) = 6z
6z2 −5z+1
signal x (n). Find ∑∞ n=−∞ x ( n ) in this case.
3 1
y ( n ) − y ( n − 1) + y ( n − 2 ) = x ( n ). (4.24)
4 8
to the input signal x (n) = nu(n) by:
(a) a direct solution in the time domain.
(b) using the z-transform.
The initial conditions are y(n) = 0 for n < 0, that is y(0) = x (0) = 0 and
y(1) = 3y(0)/4 + x (1) = 1.
Exercise 4.10. A causal discrete system is described by the difference equa-
tion
5 1
y ( n ) − y ( n − 1) + y ( n − 2 ) = x ( n ). (4.25)
6 6
If the input signal is x (n) = 1/4n u(n) find the output signal if the initial
value of the output was y(0) = 2.
Ljubiša Stanković Digital Signal Processing 215
Hint: Since y(0) does not follow from (4.25) obviously the system
output was "preloaded" before the input is applied. This fact can be taken
into account by changing the input signal at n = 0 to produce the initial
output. It is x (n) = 1/4n u(n) + δ(n). Now the initial conditions are y(0) = 2
and y(1) = 5/3 + 1/4 = 23/12 and we can apply the z-transform with this
new input signal.
1
x ( n + 2) − x ( n + 1) + x ( n ) = 0
2
with initial condition x (0) = 0 and x (1) = 1/2. The signal x (n) is causal.
and r = 0.9999 plot the amplitude of the frequency response and find the
output to the signal
T
RANSFORMATION
discrete-time systems is of high importance. Some discrete-time sys-
tems are designed and realized in order to replace or perform as
equivalents of continuous-time systems. It is quite common to design a
continuous-time system with desired properties, since the designing pro-
cedures in this domain are simpler and well developed. In the next step
the obtained continuous-time system is transformed into an appropriate
discrete-time system.
Consider an Nth order linear continuous-time system described by a
differential equation with constant coefficients
217
218 From Continuous to Discrete Systems
Δt t n
Figure 5.1 Sampling of the impulse response for the impulse invariance method.
h(n) = hc (n∆t)∆t.
Obviously this relation can be used only if the sampling theorem is satisfied
for the sampling interval ∆t. It means that the frequency response of the
continuous-time system satisfies the condition
and ∆t < π/Ωm . Otherwise the discrete-time version will not correspond to
the continuous-time version of the frequency response. Here, the discrete-
time system frequency response is related to a periodically extended form
of the continuous-time system frequency response H (Ω) as
∞
∑ H (Ω + 2kπ/∆t) = H (e jω ), Ω = ω/∆t.
k=−∞
a N s N + ... + a1 s + a0 k1 k2 kM
H (s) = = + + ··· + , (5.1)
b M s M + ... + b1 s + b0 s − s1 s − s2 s − sM
Ljubiša Stanković Digital Signal Processing 219
where only simple poles of the transfer function are assumed. The case of
multiple poles will be discussed later. The inverse Laplace transform of a
causal system, described by the previous transfer function, is
h c ( t ) = k 1 e s1 t u ( t ) + k 2 e s2 t u ( t ) + · · · + k M e s M t u ( t ).
h(n) = hc (n∆t)∆t = [k1 ∆tes1 n∆t u(n) + k2 ∆tes2 n∆t u(n) + ... + k M ∆tes M n∆t u(n)],
since u(n∆t) = u(n). The z-transform of the impulse response h(n) of the
discrete-time system is
k1 ∆t k2 ∆t k M ∆t
H (z) = −
+ −
+ ··· + . (5.2)
1−e s 1 ∆t z 1 1−e s 2 ∆t z 1 1 − es M ∆t z−1
By comparing (5.1) and (5.2) it can be concluded that the terms in the
transfer functions are transformed from the continuous-time to the discrete-
time case as
ki k i ∆t
→ . (5.3)
s − si 1 − esi ∆t z−1
If a multiple pole, of an (m + 1)th order, exists in the continuous-time
system transfer function then it holds
ki 1 dm k i
= .
( s − si ) m + 1 m! dsim s − si
si → esi ∆t .
s=jΩ
j2π/Δt
jπ/Δt jω
z=e
Im{z}
Im{s}
-j π/Δt
-j2 π/Δt
Re{s} Re{z}
forms assume that the discrete-time impulse response h(n) = hc (t)|t=+0 . Re-
mind that the theory of Fourier transforms in this case states that the inverse
Fourier transform IFT P { H ( jΩ)} = hc (t) where
Q the signal hc (t) is continuous
and IFT{ H ( jΩ)} = hc (t)|t=−0 + hc (t)|t=+0 /2 at the discontinuity points,
in this case at t = 0. The special case of discontinuity at t = 0 can be easily
detected by mapping H (s) into H (z) and by checking, for a causal system,
is the following relation satisfied
with
k1 = H (s)(s + 1)|s=−1 = −1,
'
1 ''
k2 = H (s)(s + )' = 2.
2 s=−1/2
Thus, we get
−1 2
H (s) = + .
s+1 s + 12
According to (5.3) the discrete-time system is
−1 2
H (z) = + .
1 − e −1 z −1 1 − e−1/2 z−1
Since limz→∞ H (z) = 1 obviously there is a discontinuity in the impulse
response and the resulting transfer function should be corrected as
−1 2
H (z) = + − 1/2.
1 − e −1 z −1 1 − e−1/2 z−1
Impulse and frequency responses of the systems with uncorrected and cor-
rected discontinuity effect are presented in Fig.5.3.
with
k1 = H (s)(s + 1/2)|s=−1/2 = −7,
k2 = 27/8,
222 From Continuous to Discrete Systems
0.5 0.5
0 0
-5 0 5 10 15 -5 0 5 10 15
4 jω
4 jω
|H(e )| |H(e )|
3 |H(jΩ)| 3 |H(jΩ)|
2 2
1 1
0 0
-2 0 2 -2 0 2
Figure 5.3 Impulse responses of systems in continuous and discrete-time domains (top). Am-
plitude of the frequency response of systems in continuous and discrete-time domains (bot-
tom). System without discontinuity correction (left) and system with discontinuity correction
(right).
'
'
k3 = H (s)(s + 1)2 ' = 5/4.
s=−1
The coefficient k4 follows, for example, from
H (0) = 1 = 2k1 + 3k2 + k3 + k4 ,
as
k4 = 29/8.
Thus, we get
−7 27/8 5/4 29/8
H (s) = 1
+ 1
+ 2
+ .
s+ 2 s+ 3 ( s + 1) s+1
According to (5.3) and (5.4) the discrete-time system is
−7 27/8
H (z) = +
1 − e−1/2 z−1 1 − e−1/3 z−1
'
d 5/4 ' 29/8
+ { } ' +
s
dsi 1 − e i z − 1 ' 1 − e −1 z −1
si =−1
−7z 27z/8 5e−1 z/4 29z/8
= + + + .
z−e − 1/2 z−e − 1/3 ( z − e −1 )2 z − e −1
Ljubiša Stanković Digital Signal Processing 223
s=jΩ
jω
z=e
Im{z}
Im{s}
1
-1 2/32/3 1.9894
Re{s} Re{z}
Figure 5.4 Pole-zero locations in the s-domain and the z-domain using the impulse invariance
method.
we can easily see that the poles are mapped according to s pi → es pi ∆t , Fig.5.4,
while there is no direct correspondence among zeros of the transfer functions.
Impulse responses of continuous-time system and discrete-time system are
presented in Fig.5.5.
"∞ ∞
X (s) = x (t)e−st dt ∼
= ∑ x (n)e−sn∆t = X (z)|z=es∆t .
−∞ n=−∞
This approximation leads to a relation between the Laplace domain and the
z-domain in the form of
z = es∆t .
224 From Continuous to Discrete Systems
0.3
h (t), h(n)
c
0.2
0.1
-0.1
0 5 10 15 20 25 30 35 40
jω
|H(jΩ)|, |H(e )|
1
0.5
0
-3 -2 -1 0 1 2 3
1
10
20log|H(jΩ)|
jω
0 20log|H(e )|
10
-1
10
-2
10
-3 -2 -1 0 1 2 3
Figure 5.5 Impulse responses of systems in continuous and discrete-time domains (top).
Amplitude of the frequency response of systems in continuous and discrete-time domains
(middle). Amplitude of the frequency response of systems in continuous and discrete-time
domains in logarithmic scale (bottom).
If we use this relation to map all zeros and poles of a continuous system
transfer function
z0i = es0i ∆t
z pi = es pi ∆t ,
Ljubiša Stanković Digital Signal Processing 225
s=jΩ
j2π/Δt
jπ/Δt jω
z=e
Im{z}
Im{s}
-j π/Δt
-j2 π/Δt
Re{s} Re{z}
Figure 5.6 Illustration of the zeros and poles mapping in the matched z−transform method.
dx (t)
y(t) =
dt
∼ n∆t) − x ((n − 1)∆t) .
y(n∆t) =
x (
∆t
The Laplace transform domain of the continuous-time first derivative is
1 − z −1
Y (z) = X ( z ). (5.6)
∆t
Based on (5.5) and (5.6) we can conclude that a mapping of the correspond-
ing differentiation operators from the continuous-time to the discrete-time
domain is
1 − z −1
s= . (5.7)
∆t
With a normalized discretization step ∆t = 1 this mapping is of the form
s = 1 − z −1 .
"t −∆t
t"
y(t) = x (t)dt ∼
= x (t)dt + x (n∆t)∆t.
−∞ −∞
y(n∆t) ∼
= y(n∆t − ∆t) + x (n∆t)∆t.
Ljubiša Stanković Digital Signal Processing 227
The Laplace and the z-transform domain forms of the previous integral
equations are
1
Y (s) = X (s)
s
∆t
Y (z) = X ( z ).
1 − z −1
1 − s → z −1 . (5.8)
Now we will consider the region that corresponds to the imaginary axis and
the left semi-plane of the s-domain (containing poles of a stable system),
Fig.5.7(left). The aim is to find the corresponding region in the z-domain.
If we start from the s-domain and the region in Fig.5.7(left), the first
mapping is to reverse the s-domain to −s and shift it for +1, as
1 − s → p.
1
Re{ p} = Re{ }
z
1
1 = Re{ }
x + jy
1 x − jy
1 = Re{ }
x + jy x − jy
228 From Continuous to Discrete Systems
s=0+jΩ p=1
z=ejω
Im{p}
Im{z}
Im{s}
-1
1-s → p p→ z
Figure 5.7 Illustration of the differentiation based mapping of the left s−semi-plane with the
imaginary axis (left), translated and reversed p−domain (middle), and the z−domain (right).
resulting in
x
1=
x2 + y2
or in * +2
1 1
( x − )2 + y2 = . (5.9)
2 2
Therefore, the imaginary axis in the s-plane is mapped onto a circle defined
by (5.9), Fig.5.7(right) in the z-plane. From the mapping relation 1 − s → z−1
it is easy to conclude that the origin s = 0 + j0 maps into z = 1 and that
s = 0 ± j∞ maps into z = ±0, according to 1/ (1 − s) → z.
Mapping of the imaginary axis into z-domain can also be analyzed
from
1 − (re jω )−1 1 − r −1 cos ω r −1
σ + jΩ → = +j sin ω.
∆t ∆t ∆t
For σ = 0 follows
1 − r −1 cos ω = 0 (5.10)
r = cos ω,
with
r −1 tan ω
Ω=sin ω = .
∆t ∆t
Obviously ω = 0 maps to Ω = 0 (with Ω ∼ = ω/∆t for small ω), and ω =
±π/2 maps into Ω → ±∞. Thus, the whole imaginary axis maps onto
−π/2 ≤ ω ≤ π/2. These values of ω could be used within the basic period.
Relation (5.10), with −π/2 ≤ ω ≤ π/2, is a circle defined by (5.9) if we
Ljubiša Stanković Digital Signal Processing 229
, ,
replace r = x2 + y2 and cos ω = x/ x2 + y2 with σ < 0 (semi-plane with
negative real values) being mapped into r < cos ω (interior of unit circle).
Example 5.4. A continuous-time system is described by a differential equation
3 1
y′′ (t) + y′ (t) + y(t) = x (t),
4 8
with zero initial conditions and the transfer function
1
H (s) = .
s + 4 s + 18
2 3
P ⋆A discrete-time
Q system transfer function is obtained by replacing
s = 1 − z−1 /∆t in H (s) as
1
H (z) = B C2
1 − z −1 3 1 − z −1 1
∆t + 4 ∆t + 8
(∆t)2
= 2
1 + 34 ∆t + 18 (∆t) − [2 + 34 ∆t]z−1 + z−2
with
y(n) = B0 x (n) + A1 y(n − 1) + A2 y(n − 2)
(∆t)2
B0 = = 0.1778
1 + 34 ∆t + 18 (∆t)2
[2 + 34 ∆t]
A1 = 3 1 2
= 1.6889
1+ 4 ∆t + 8 (∆t )
1
A2 = − = −0.7111,
1 + 34 ∆t + 18 (∆t)2
where ∆t = 1/2. For x (t) = u(t) in the continuous-time case
1
Y (s) = H (s) X (s) =
s(s2 + 34 s + 18 )
8 8 16
= + 1
− 1
s s+ 2 s+ 4
with
y(t) = [8 + 8e−t/2 − 16e−t/4 ]u(t).
The results of the difference equation for y(n) are compared with the exact
solution y(t) in Fig.5.8. The agreement is high. It could be additionally
improved by reducing the sampling interval, for example, to ∆t = 1/8.
230 From Continuous to Discrete Systems
10
y(t), y(n)
0 5 10 15
Figure 5.8 Exact solution of the difference equation y(t) in solid line and the discrete-time
system output y(n) in large dots for ∆t = 1/2 and in small dots for ∆t = 1/8..
"t −∆t
t"
x (n∆t) + x ((n − 1)∆t)
y(t) = x (t)dt ∼
= x (t)dt + ∆t
2
−∞ −∞
x ( n ) + x ( n − 1)
y ( n ) = y ( n − 1) + ∆t.
2
In the Laplace and the z-transform domain, these relations have the forms
1
Y (s) = X (s)
s
∆t 1 + z−1
Y (z) = X ( z ).
2 1 − z −1
Ljubiša Stanković Digital Signal Processing 231
2 1 − z −1
s→ . (5.11)
∆t 1 + z−1
y ( n ) = x ( n ) − x ( n − 1 ).
The same signal samples can used for the first-order forward derivative
approximation
y ( n − 1 ) = x ( n ) − x ( n − 1 ).
If we assume that the difference x (n) − x (n − 1) fits better to the mean
of y(n) and y(n − 1) than to any single one of them, then the derivative
approximation by using the difference equation
y ( n ) + y ( n − 1)
= x ( n ) − x ( n − 1 ),
2
produces the bilinear transform.
In order to prove that the imaginary axis in the s−domain corresponds
to the unit circle in the z−domain we may simply replace z = e jω into (5.11)
and obtain
1 − e− jω e jω/2 − e− jω/2 ω
2 −
= 2 jω/2 = 2j tan( ) → s∆t.
1+e jω e + e− jω/2 2
For s = σ + jΩ follows
σ=0
2 ω
Ω= tan( ).
∆t 2
2 ω ω
Ω= tan( ) ∼
= , for |ω | ≪ 1.
∆t 2 ∆t
232 From Continuous to Discrete Systems
From
s∆t
1+ 2
z= s∆t
1− 2
F
(1 + σ∆t 2
2 ) + ( Ω∆t
2 )
2
|z| = F
(1 − σ∆t 2
2 ) + ( Ω∆t
2 )
2
it may easily be concluded that σ < 0 maps into |z| < 1, since 1 + σ∆t
2 <
σ∆t
1 − 2 for σ < 0.
The bilinear transform mapping can be derived by using a series of
complex plane mappings. Since
s∆t
1+ 2 2
z= s∆t
= − 1,
1− 2 1 − s∆t
2
we can write
s∆t
1− → p1 ,
2
1
→ p2 ,
p1
2p2 − 1 → z.
s=0+jΩ p=1
Im{p1}
Im{s}
-1
p1 → p2
Re{s} Re{p1}
jω
z=e
Im{p }
Im{z}
1 1
2
Re{p2} 2p2-1 → z
Re{z}
Figure 5.9 Bilinear mapping illustration trough a series of elementary complex plane map-
pings.
and to stop all other possible signal components. The parameters are Q =
0.01, Ω1 = π/4, and Ω2 = 3π/5. The signal is sampled with ∆t = 1 and
the discrete-time signal x (n) is formed. Using the bilinear transform, design
the discrete system that corresponds to the continuous-time system with the
transfer function H (s).
⋆For the beginning just use the bilinear transform relation
1 − z −1
s→2 (5.12)
1 + z −1
and map H (s) to HB (z) without any pre-modification. The result is presented
in the first two subplots of Fig.5.10. The discrete frequencies are shifted since
the bilinear transform (5.12) made a nonlinear frequency mapping from the
234 From Continuous to Discrete Systems
0.016569(1 + z−1 )2
=
4.65327z−2 − 6.6272z−1 + 4.7195
0.0551(1 + z−1 )2
+
11.4677z 2 + 7.1556z−1 + 11.6879
−
0.003567(1 + z−1 )2
= −1
(z − 1.0071e j0.25π )(z−1 − 1.0071e− j0.25π )
0.0048(1 + z−1 )2
+
( z −1 − 1.0096e j0.6π )(z−1 − 1.0096e− j0.6π )
Ljubiša Stanković Digital Signal Processing 235
Ω1 Ω2
1
H(s)
s → 2(1-z
0.5
0
-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1
-1
)/(1+z )
1
HB(z)
-1
0.5
0
-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1
frequency ω/π or ΩΔt/π
s → 2(1-z
0.5
0
-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1
-1
)/(1+z )
1
H(z) ω =Ω Δt ω2=Ω2Δt
1 1
-1
0.5
0
-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1
frequency ω/π or ΩΔt/π
Figure 5.10 Amplitude of the continuous-time system with transfer function H (s) and the
amplitude of the transfer function HB (z) of the discrete-time system obtained by the bilinear
transform (first two subplots). A premodified system to take into account the frequency map-
ping nonlinearity in the bilinear transform Hd (s) and the amplitude of the transfer function
H (z) of the discrete-time system obtained by the bilinear transform of Hd (s) (last two subplots).
Sampling
Fourier transform
Method theorem
H (s)|s= jΩ → H (z)|z=e jω
condition
Impulse Invariance Yes, Ω = ω/∆t Yes
Matched z-transform No No
First-oder difference No No
tan(ω/2)
Bilinear transform Yes, Ω = ∆t/2 No
|H(ejω)|2
ω ω
-2 π -π - ωc ωc π 2π -2 π -π 0 π 2π
Figure 5.11 Lowpass filter frequency response: ideal case (left) and Butterworth type (right).
Ljubiša Stanković Digital Signal Processing 237
k i = H (s)(s − si )|s=si
k0 = (−0.3628 + j0.1503)/∆t,
k1 = (0.3628 − j0.8758)/∆t,
k2 = (0.3628 + j0.8758)/∆t,
k3 = (−0.3628 − j0.1503)/∆t.
Using the impulse invariance method we get the transfer function of the
discrete-time fourth order Butterworth filter
k0 ∆t k1 ∆t k2 ∆t k3 ∆t
H (z) = + + +
1 − es0 ∆t z−1 1 − es1 ∆t z−1 1 − es2 ∆t z−1 1 − es3 ∆t z−1
−0.3628 + j0.1503 0.3628 − j0.8758
= +
1 − eωc (−0.3827+ j0.9239) z−1 1 − eωc (−0.9239+ j0.3827) z−1
0.3628 + j0.8758 −0.3628 − j0.1503
+ + .
1 − eωc (−0.9239− j0.3827) z−1 1 − eωc (−0.3827− j0.9239) z−1
238 From Continuous to Discrete Systems
ωd = Ωd ∆t = 0.8284
H (z) =
ωd 4
z −1 2 − z −1 −1 −1
[4( 11−
+ z −1
) + 2ωd 0.7654 11+ z −1
−z )2 + 2ω 1.8478 1−z + ω 2 ]
+ ωd2 ][4( 11+ z −1 d 1 + z −1 d
0.4710
= z −1 2 − z −1 −1 −1
[4( 11−
+ z −1
) + 0.6863][4( 11−
+ 1.2626 11+ z −1
z
+ z −1
−z + 0.6863]
)2 + 3.0481 11+ z −1
P − 1
Q 4
0.4710 1 + z
=P QP Q
−
3.4237z − 6.6274z + 5.9484 1.6382z−2 − 6.6274z−1 + 7.7704
2 − 1
P Q4
0.084 1 + z−1
P
= −2 Q P Q
z − 1.9357z−1 + 1.7343 z−2 − 4.0455z−1 + 4.7433
0.084z−4 + 0.336z−3 + 0.504z−2 + 0.336z−1 + 0.084
= .
z−4 − 5.9810z−3 + 14.3z−2 − 16.1977z−1 + 8.2263
The transfer function (amplitude and phase) of the continuous-time
filter and the discrete-time filters obtained by using the impulse invariance
method and the bilinear transform are presented in Fig.5.12, within one
Ljubiša Stanković Digital Signal Processing 239
1.5
arg{H(e }
1
jω
|H(e |
jω
0
0.5
-π
0 ω ω
-π - π/2 0 π/2 π -π - π/2 0 π/2 π
Figure 5.12 Amplitude and phase of the fourth order Butterworth filter frequency response
obtained by using the impulse invariance method and bilinear transform.
2 1 − z −1
s=
∆t 1 + z−1
Ljubiša Stanković Digital Signal Processing 241
as
1.13273
H (z) = B C2
1 − z −1 − z −1 −z−1 + 1.13272 )
(2 1 + z −1
+ 1.1327)( 2 11+ z − 1 + 2.2653 11+ z −1
1.4533(1 + z−1 )3
=
(−0.8673z−1 + 3.1327)(3.0177z−2 − 5.434z−1 + 7.54)
−0.5553z−3 − 1.6658z−2 − 1.6658z−1 − 0.5553
=
z−3 − 5.4127z−2 + 9.0028z−1 − 9.0249
0.0615z3 + 0.1846z2 + 0.1846z + 0.0615
= .
z3 − 0.9975z2 + 0.5998z − 0.1108
The corresponding difference equation is
y(n) = 0.9975y(n − 1) − 0.5998y(n − 2) + 0.1108y(n − 3)
+ 0.0615x (n) + 0.1846x (n − 1) + 0.1846x (n − 2) + 0.0615x (n − 3).
Kπ ≤ Ω0 ∆t < (K + 1)π
will, after sampling, result into a component within the basic period of the
Fourier transform of discrete-time signal, corresponding to the continuous
K
signal at exp( j(Ω0 t − ∆t πt) This effect is known as aliasing. The most
obvious visual effect is when a wheel rotating with f 0 = 25 [Hz], Ω0 = 50π, is
sampled in a video sequence at ∆t = 1/50 [sec]. Then Ω0 ∆t = π corresponds
to exp( j(Ω0 t − 50πt)) = e j0 , i.e., the wheel looks as a static (nonmoving)
object.
HH (e jω ) = H (e j(ω −π ) ).
|H (ejω)|2
jω 2
1
|H(e )|
H
ω ω
-2 π -π -ω ω π 2π -2 π -π 0 π 2π
c c
x(n) y(n)
× h(n) ×
n n
(-1) (-1)
1.5 1.5
1 1
|HH(ejω|
|H(e |
jω
0.5 0.5
0 ω 0 ω
-π - π/2 0 π/2 π -π - π/2 0 π/2 π
Figure 5.15 Amplitude of frequency response of a lowpass Butterworth filter (left) and a filter
obtained from the lowpass Butterworth filter when z is replaced by −z (right).
jω 2
|H(ejω)|2
|HB(e )|
ω ω
-2 π -π - ωc ωc π 2π -2 π - π - ω0 0 ω0 π 2π
× h(n) ×
× h(n) ×
cos(ω0n) 2cos(ω0n)
Figure 5.17 Bandpass system realization using corresponding lowpass systems and signal
modulation.
∞ ∞
y(n) = h B (n) ∗ x (n) = ∑ h B (m) x (n − m) = 2 ∑ cos(ω0 m)h(m) x (n − m)
m=−∞ m=−∞
∞
=2 ∑ cos(ω0 n + ω0 m − ω0 n)h(m) x (n − m)
m=−∞
∞
=2 ∑ [cos(ω0 n) cos(ω0 m − ω0 n) − sin(ω0 n) sin(ω0 m − ω0 n)]h(m) x (n − m)
m=−∞
∞
= 2 cos(ω0 n) ∑ cos(ω0 (n − m)) x (n − m)h(m)
m=−∞
∞
+2 sin(ω0 n) ∑ sin(ω0 (n − m)) x (n − m)h(m).
m=−∞
z+2 z − 1a e jθ
Hs (z) = H (z) H A (z) = e− j2θ .
(z − 12 )(z − 13 )(z − 2) 1 − 1a e− jθ z
For a = 1/2 and θ = 0 we get
z+2 z−2
Hs (z) =
(z − 12 )(z − 13 )(z − 2) 1 − 2z
z+2
=− .
2(z − 12 )2 (z − 13 )
This system has the same frequency response amplitude as the initial system
' ' ' ' ' '
' ' ' ' ' '
'Hs (e jω )' = 'H (e jω ) H A (e jω )' = 'H (e jω )' .
Ljubiša Stanković Digital Signal Processing 247
where 0 < ai < 1 and θi , i = 1, 2, ..., N are arbitrary constants and phases. The
resulting frequency response amplitude is
' '
' '
'H A (e jω )' = 1.
This system can be used for multiple poles cancellation and phase correc-
tion.
1
Hi (z) = .
H (z)
It is obvious that
H (z) Hi (z) = 1
h ( n ) ∗ h i ( n ) = δ ( n ).
This kind of system can be used to reverse the signal distortion. For ex-
ample, assume that the Fourier transform of a signal x (n) is distorted dur-
ing transmission by a transfer function H (z), i.e., the received signal z-
transform is R(z) = H (z) X (z). In that case the distortion can be compen-
sated by processing the received signal using the inverse system. The output
signal is obtained as
1
Y (z) = R ( z ) = X ( z ).
H (z)
The system Hi (z) = 1/H (z) should be stable as well. It means that the poles
of the inverse system should be within the unit circle. The poles of the
inverse system are equal to the zeros of H (z).
The system H (z) whose both poles and zeros are within the unit circle
is called a minimum phase system.
248 From Continuous to Discrete Systems
If a system is the minimum phase system (with all poles and zeros
within |z| < 1) then this system has a minimum group delay out of all
systems with the same amplitude of the frequency response. Thus, any
nonminimum phase system will have a more negative phase compared
to the minimum phase system. The negative part of the phase is called
the phase-lag function. The name minimum phase system comes from the
minimum phase-lag function.
In order to prove this statement consider a system H (z) with the sam-
ple amplitude of the frequency response as a nonminimum phase system
Hmin (z). Its frequency response can be written as
z−1 − ae− jθ
H (z) = Hmin (z) H A (z) = Hmin (z)
1 − ae jθ z−1
Here we assumed the first-order allpass system without any loss of gener-
ality, since the same proof can be used for any number of allpass systems
that multiply Hmin (z). Since 0 < a1 < 1 and the system Hmin (z) is stable the
system H (z) has a zero at |z| = 1/a1 > 1.
The phases of the system are related as
e− jω − ae− jθ
arg{ H A (e jω )} = arg{ }
1 − ae jθ e− jω
1 − ae− jθ e jω
= arg{e− jω } = −ω + arg{1 − ae− jθ e jω }
1 − ae jθ e− jω
a sin(ω − θ )
− arg{1 − ae jθ e− jω } = −ω − 2 arctan .
1 − a cos(ω − θ )
d arg{ H A (e jω )} a cos(ω − θ ) − a2
τgA (ω ) = − =1+2
dω 1 − 2a cos(ω − θ ) + a2
1 − a2 1 − a2
= 2
=' ' .
1 − 2a cos(ω − θ ) + a '1 − ae j(ω −θ ) '2
250 From Continuous to Discrete Systems
τg (ω ) = τg min (ω ) + τgA (ω )
τg (ω ) ≥ τg min (ω ),
1 − ae− jθ
arg{ H A (e j0 )} = arg{ }=0 (5.16)
1 − ae jθ
"ω
arg{ H A (e jω )} = − τg (ω )dω ≤ 0 (5.17)
0
since τg (ω ) > 0 for 0 ≤ ω < π.
We can conclude that the minimum phase systems satisfy the following
conditions.
1. A minimum phase system is system of minimum group delay out
of the systems with the same amplitude of frequency response. A system
containing one or more allpass parts with uncompensated zeros outside of
the unit circle will have larger delay than the system which does not contain
zeros outside the unit circle.
2. The phase of a minimal phase system will be lower than the phase
of any other system with the same amplitude of frequency response since,
according to (5.17),
This proves the fact that the phase of any system arg { H (e jω ) is always
lower than the phase of minimal phase system arg { Hmin (e jω )}, having the
same amplitude of the frequency response.
3. Since the group delay is minimal we can conclude that
n n
∑ |hmin (m)|2 ≥ ∑ |h(m)|2
m =0 m =0
This relation may be proven in a similar way like minimal phase property,
by considering the outputs of a minimum phase system and a system
H (z) = Hmin (z) H A (z).
Ljubiša Stanković Digital Signal Processing 251
Example 5.12. A system has absolute squared amplitude of the frequency response
equal to
B C2
' '2 5
' ' 2 cos(ω ) + 2
'H (e jω )' =
(12 cos(ω ) + 13)(24 cos(ω ) + 25)
Find the corresponding minimal phase system.
⋆ For the system we can write
' '2
' '
'H (e jω )' = H (e jω ) H ∗ (e jω ) = H (e jω ) H (e− jω )
In the z−domain the system with this amplitude of the frequency response
(with real-valued coefficients) satisfies
' ' ' '2
∗ 1 '
' 1 '' ' '
H (z) H ( ∗ )' = H (z) H ( )' = 'H (e jω )' = H (e jω ) H (e− jω ).
z z=e jω z z=e jω
In this sense
B C2
'
'
'2
' e jω + e− jω + 5
2
'H (e jω )' =
(6e jω + 6e− jω + 13)(12e jω + 12e− jω + 25)
and
B C2
1 z+ 5
2 + z −1
H (z) H ( ) =
z (6z + 13 + 6z−1 )(12z + 25 + 12z−1 )
B C2
z2 + 52 z + 1
=
(6z2 + 13z + 6)(12z2 + 25z + 12)
1 (z + 2)2 (z + 12 )2 1 ( 1z + 12 )2 (z + 12 )2
= = .
36 (z + 23 )(z + 32 )(z + 34 )(z + 43 ) 36 (z + 23 )( 1z + 23 )(z + 34 )( 1z + 34 )
The minimum phase system, with the desired amplitude of the frequency
response, is a part of H (z) H ∗ ( z1∗ ) with zeros and poles inside the unit circle
1 (z + 12 )2
H (z) = .
6 (z + 23 )(z + 34 )
5.6 PROBLEMS
2s
H (s) = − .
s2 + 2s + 2
(1 + 4s)
H (s) = .
(s + 1/2)(s + 1)3
2QΩ1
H (s) =
s2 + 2Ω1 Qs + Ω21 + Q2
Ljubiša Stanković Digital Signal Processing 253
x (t) = A1 cos(Ω1 t + ϕ1 )
and to stop all other possible signal components. The parameters are Q =
0.01, Ω1 = π/2. The signal is sampled with ∆t = 1 and a discrete-time signal
x (n) is formed. Using bilinear transform design the discrete system that
corresponds to the continuous-time system with transfer function H (s).
Problem 5.7. (a) By using the bilinear transform find the transfer function of
the second-order Butterworth filter with f ac = 4kHz. The sampling interval
is ∆t = 50µ sec.
(b) Translate the discrete-time transfer function to obtain a highpass filter.
Find its corresponding critical frequency in the continuous-time domain.
Problem 5.8. Design a discrete-time lowpass Butterworth filter for the sam-
pling frequency 1/∆t = 10 kHz. The passband should be from 0 to 1 kHz,
maximal attenuation in the passband should be 3 dB and the attenuation
should be more than 10 dB for frequencies above 2 kHz.
Problem 5.9. Using the impulse invariance method design a Butterworth
filter with the passband frequency ω p = 0.1π and stopband frequency
ωn = 0.3π in the discrete domain. Maximal attenuation in the passband
region should be less than 2dB, and the minimal attenuation in the stopband
should be 20dB.
Problem 5.10. Highpass filter can be obtained from a lowpass by using
HH (s) = H (1/s). Using the bilinear transform with ∆t = 2 we can trans-
form the continuous-time domain function into discrete domain using the
relation s = (z − 1)/(z + 1). If we have a design of a lowpass filter how to
change its coefficients in order to get a highpass filter.
Problem 5.11. For filtering of a continuous-time signal a discrete-time filter
is used. Find the corresponding continuous-time filter frequencies if the
discrete-time filter is: a) a lowpass with ω p = 0.15π, b) bandpass within
0.2π ≤ ω ≤ 0.25π, c) a highpass with ω p = 0.35. Consider cases when
∆t = 0.001s and ∆t = 0.1s.
What should be the starting frequencies to design these systems in the
continuous-time domain if the impulse invariance method is used and what
are the design frequencies if the bilinear transform is used?
Problem 5.12. A transfer function of the first-order lowpass system is
1−α
H (z) = .
1 − αz−1
254 From Continuous to Discrete Systems
Problem 5.13. Using allpass system find stable systems with the same
amplitude of the frequency response as the systems:
(a)
2 − 3z−1 + 2z−2
H1 (z) =
1 − 4z−1 + 4z−2
(b)
z
H2 (z) = .
(4 − z)(1/3 − z)
Problem 5.14. The z-transform
Problem 5.15. A signal x (n) has passed trough a media whose influence
can be described by the transfer function
√
(4 − z)(1/3 − z)(z2 − 2z + 14 )
H (z) = .
z − 12
5.7 SOLUTIONS
1
LC 25
H (s) = =
s2 + s RL + LC
1 s2 + 8s + 25
25
=
(s + 4 + 3j)(s + 4 − j3)
Ljubiša Stanković Digital Signal Processing 255
j 25
6 − j 25
6
H (s) = + .
s + 4 + j3 s + 4 − j3
The poles are mapped using
s i → zi = e si .
j 25
6 − j 25
6
H (z) = +
1 − e−(4+ j3) z−1 1 − e−(4− j3) z−1
25 −4 −1
3 e z sin 3
= ,
1 − 2e cos 3z 1 + e−8 z−2
− 4 −
25 −4
y(n) = e sin(3) x (n − 1) + 2e−4 cos(3)y(n − 1) − e−8 y(n − 2).
3
The output signal values can be calculated for any input signal using
this difference equation. For x (n) = δ(n) the impulse response would follow.
The impulse response can be obtained in a closed form from
as
25 −4n − j3n
h(n) = e ( je − je j3n )u(n) =
6
25
= e−4n sin(3n)u(n).
3
Solution 5.2. The system is not of lowpass type. For s → ∞ we get H (s) → 1.
Thus, the impulse invariance method cannot be used. The bilinear trans-
form can be used. It produces
(1 − z −1 )2 −z + 3−1
4 (1+z−1 )2 − 6 11+ z −1 13z−2 − 2z−1 + 1
H (z) = = .
(1 − z −1 )2 −z + 3
4 (1+z−1 )2 + 6 11+
−1 z−2 − 2z−1 + 13
z −1
256 From Continuous to Discrete Systems
3 1
y′′ (t) + y′ (t) + y(t) = x (t)
2 2
the transfer function is
1
H (s) = 3 1
.
s2 + 2s + 2
Corresponding discrete system is obtained using
1 − z −1
s→ = 10(1 − z−1 )
∆t
as
1
H (z) =
100(1 − z−1 )2 + 32 10(1 − z−1 ) + 1
2
1
= 231
.
100z−2 − 215z−1 + 2
2 430 200
y(n) = x (n) + y ( n − 1) − y ( n − 2 ).
231 231 231
Solution 5.4. The transfer function can be written as
1+j 1−j
H (s) = − − .
s+1−j s+1+j
1 − z −2
H ( z ) = −2 .
5 − 2z−1 + z−2
Solution 5.5. (a) The transfer function
(1 + 4s)
H (s) =
(s + 1/2)(s + 1)3
Ljubiša Stanković Digital Signal Processing 257
k1 k2 k3 k4
H (s) = + + +
s + 1/2 (s + 1) (s + 1)2 ( s + 1 )3
'
with k1 = H (s)(s + 1/2)|s=−1/2 = −8 and k4 = H (s)(s + 1)3 's=−1 = 6. By
equating the coefficients with s3 to 0 we get the relation k1 + k2 = 0. Similar
relation follows for the coefficients with s2 as 3k1 + 5k2 /2 + k3 = 0 or
k1 /2 + k3 = 0. Then k2 = 8 and k3 = 4. With
ki ki
→
s − si 1 − e s i z −1
and
1 dm k i 1 dm ki
m → m { }
m! dsi s − si m! dsi 1 − esi z−1
we get the discrete system
−8 8
H (z) = +
1 − e−1/2 z−1 1 − e −1 z −1
* +' * +'
d 4 '
' d2 6 '
'
+ s − 1 ' + 2 s − 1 '
ds1 1 − e 1 z s1 =−1 ds1 1 − e 1 z s1 =−1
−8 8 4e−1 z−1 3e−2 z−2 + 3e−1 z−1
= + − −
+ − −
+
1 − e−1/2 z−1 1−e z1 1 1
(1 − e z ) 1 2 (1 − e −1 z −1 )3
−5.83819z−3 − 9.68722z−2 + 22.0531z−1
=
(z−1 − e)3 (z−1 − e1/2 )
−1
(1 + 8 11− z
+ z −1
)
H (z) = −1 −1
−z + 1/2)(2 1−z
(2 11+ + 1 )3
z −1 1 + z −1
−14z−4 − 24z−3 + 12z−2 + 40z−1 + 18
= .
3z−4 − 32z−3 + 126z−2 − 216z−1 + 135
Solution 5.6. Since we use the bilinear transform we have to pre-modify the
system according to
2 Ω ∆t
Ωd = tan( 1 ) = 2.0 = 0.6366π.
∆t 2
The frequency value is shifted from Ω1 = 0.5π to Ωd = 0.6366π. The modi-
fied system is
2QΩd
Hd (s) = 2 .
s + 2Ωd Qs + Ω2d + Q2
−z −1
Now using s = 2 11+ z −1
the corresponding discrete- system is obtained,
2QΩd
H (z) = B C2 B C .
− z −1
2 11+ + 2Ωd Q 2 11− z −1
+ Ω 2 + Q2
z −1 + z −1 d
s1 s2 4π 2 f c2
Ha ( s ) = = √ .
(s − s1 )(s − s2 ) s + 2π f c 2s + 4π 2 f c2
2
1.0548(1 + z−1 )2
H (z) = .
5.1066 − 1.8874z−1 + z−2
Ljubiša Stanković Digital Signal Processing 259
1.0548(1 − z−1 )2
H (z) = .
5.1066 + 1.8874z−1 + z−2
Solution 5.8. For the continuous-time system the design frequencies are
f p = 1 kHz
f s = 2 kHz.
They correspond to
Ω p = 2π 103 rad/s
Ωs = 4π 103 rad/s.
ω p = 0.2π
ωs = 0.4π.
The frequencies for the filter design, that will be mapped to ωs and ω p by
using the bilinear transform, are
2 0.6498
Ω pd = tan(0.2π/2) =
∆t ∆t
2 1.4531
Ωsd = tan(0.4π/2) = .
∆t ∆t
1−100.1a p
1 log 1−100.1as
N= = 1.368.
2 log Ω pd
Ωsd
We assume N = 2.
260 From Continuous to Discrete Systems
Mapping this system into the discrete-time domain by using the bilinear
transform,
2 1 − z −1
s= ,
∆t 1 + z−1
produces
0.067569(1 + z−1 )2
H (z) = .
1 − 1.14216z−1 + 0.412441z−2
Solution 5.9. The Butterworth filter order is
1−100.1a p
1 log 1−100.1as
N= = 2.335.
2 log Ω p
Ωs
Ωp
Ωc = 2N
, = 0.109345π = 0.3435.
100.1a p − 1
ki ∆tk i
→ .
s − spi 1 − es pi ∆t z−1
The discrete-time system transfer function is
−0.0253z−2 − 0.0318z−1
H (z) = .
−1.98774 + 4.61093z−1 − 3.68033z−2 + z−3
Solution 5.10. The transfer function is
1
HH (s) = H ( )
s
2 1− z −1 2 z −1
with s = ∆t 1+z−1 = ∆t z+1 and ∆t = 2. Corresponding lowpass filter would
be
z−1
HL (z) = H (s)|s= z−1 = H ( ).
z +1 z+1
The discrete highpass filter is
'
1 '
HH (z) = HH (s)|s= z−1 = H ( )''
z +1 s s = z −1
z +1
z+1
HH (z) = H (
).
z−1
Obviously HH (z) = HL (−z). It means that a discrete highpass system can be
realized by replacing z with −z in the transfer function. For ∆t ̸= 2 a scaling
is present as well.
262 From Continuous to Discrete Systems
2 − 3z−1 + 2z−2
H1 (z) =
(1 − 2z−1 )2
is not stable since it has a second-order pole at z = 2. This system may be
stabilized, keeping the same amplitude of the frequency response, using a
second-order allpass system with zero at z = 2
( )2
z −1 − 1
2
H A (z) = 1 −1
.
1− 2 z
2 − 3z−1 + 2z−2
H1 (z) = .
( z −1 − 2 )2
Causal system H2 (z) has a pole at z = 4. It can be stabilized by using allpass
system
z−1 − 14 4−z
H A (z) = 1 −1
= .
1 − 4z 4z − 1
Ljubiša Stanković Digital Signal Processing 263
(z − 14 )(z + 12 )
H (z) = .
(z + 45 )(z − 37 )
z − 12
Hi (z) = .
(4 − z)(1/3 − z)(z − 1.2071)(z − 0.2071)
These poles have to be compensated, keeping the same amplitude, by using
two first-order allpass systems. The resulting system transfer function is
z − 4 z − 1.2071
Hi (z)
1 − 4z 1 − 1.2071z
z − 12
= .
(1/3 − z)(z − 0.2071)(1 − 4z) (1 − 1.2071z)
264 From Continuous to Discrete Systems
5.8 EXERCISE
( s + 2)
H (s) = .
4s2 + s + 1
What is the corresponding discrete-time system obtained with ∆t = 1 by
using the impulse invariance method and the bilinear transform.
Exercise 5.2. A continuous system is described by a differential equation
1
y′′ (t) + 6y′ (t) − y(t) = x (t) + x ′ (t)
2
with zero initial conditions. What is the corresponding transfer function
of a discrete system obtained by using the first-order backward difference
approximation with ∆t = 1?
Exercise 5.3. (a) A continuous system
2QΩ0
H (s) =
s2 + 2Ω0 Qs + Ω20 + Q2
x (t) = A cos(Ω0 t + ϕ)
sampled with the sampling interval ∆t = 10−3 s. What would be the corre-
sponding continuous-time output signal after an ideal D/A converter.
Exercise 5.4. (a) By using the bilinear transform find the transfer function
of a third-order Butterworth filter with f ac = 3.4 kHz. The sampling step is
∆t = 40 µ sec.
(b) Translate the discrete transfer function to obtain a bandpass system
with corresponding central frequency f ac = 12.5 kHz in the continuous
domain.
Ljubiša Stanković Digital Signal Processing 265
2 − 5z−1 + 2z−2
H1 (z) = ,
1 − 4z−1 + z−2
z−1
H2 (z) = .
(2 − z)(1/4 − z)
Exercise 5.7. The z-transform
(z − 13 )(z−1 − 13 )
R(z) =
(z + 12 )(z−1 + 12 )
can can be written as
1
R(z) = H (z) H ∗ ().
z∗
Find H (z) for the minimum phase system. If h(n) is the impulse response
of H (z) and h1 (n) is the impulse response of
z−1 − a1 e− jθ1
H1 (z) = H (z)
1 − a1 e jθ1 z−1
show that |h(0)| ≤ |h1 (0)| for any θ1 and | a1 | < 1. All systems are causal.
Exercise 5.8. A signal x (n) has passed trough a media whose influence can
be described by the transfer function
(1 − z/3)(1 − 5z)(z2 − z + 34 )
H (z) =
z2 − 2/3
and the signal r (n) = x (n) ∗ h(n) is' obtained.
' Find
' a causal
' and stable system
to process r (n) in order to obtain 'Y (e jω )' = ' X (e jω )'.
266 From Continuous to Discrete Systems
Chapter 6
Realization of Discrete Systems
L
INEAR
ence equation relating the output signal with the input signal at the
considered instant and the previous values of the output and input
signal. The transfer function can be written in various forms producing dif-
ferent system realizations. Some of them will be presented next. Symbols
that are used in the realizations are presented in Fig.6.1.
a
z -1
x(n) ax(n) x(n) x(n-1) x(n) x(n)
x(n)
+ + ×
x(n) x(n)+y(n) x(n) - x(n)- y(n) x(n) x(n)y(n)
Figure 6.1 Symbols and their function in the realization of discrete-time systems.
267
268 Realization of Discrete Systems
x(n) B0 y(n)
+ +
z-1 z-1
x(n-1) + + y(n-1)
B A1
-1 1 -1
z z
x(n-2) y(n-2)
B A
2 2
z-1 z-1
+ +
B A
-1 1 1
z z-1
+ +
B A
2 2
-1 -1
z z
B A
M N
and
1
H2 (z) = .
1 − A1 z−1 − ... − A N z− N
The overall transfer function is
It means that these two blocks can interchange their positions. After the
positions are interchanged, then by using the same delay systems, we get
the resulting system in the direct realization II form, presented in Fig.6.4.
This system uses a reduced number of delay blocks in the realization.
Example 6.1. Find the transfer function of a discrete system presented in Fig.6.5.
⋆The system can be recognized as a direct realization II form. After
its blocks are separated and interchanged the system in a form presented in
Fig.6.6 is obtained.
The output of the first block is
1 1
y1 ( n ) = x ( n ) − x ( n − 1) + x ( n − 2). (6.2)
2 3
Its transfer function is
1 1
H1 (z) = 1 − z−1 + z−2 .
2 3
270 Realization of Discrete Systems
x(n) B0 y(n)
+ +
z-1
+ +
A B
1 -1 1
z
+ +
A B
2 2
z-1
A B
N M
x(n) y(n)
+ +
-1
z
+
-1/2
-1
z
+
1/2 1/3
z-1
-1/6
1
H2 (z) = .
1 − 12 z−2 + 16 z−3
Ljubiša Stanković Digital Signal Processing 271
+
-1/2
-1 -1
z z
+
1/3 1/2
-1
z
-1/6
The difference equation for the whole system is obtained after y1 (n) from
(6.2) is replaced into (6.3)
1 1 1 1
y(n) = y ( n − 2) − y ( n − 3) + x ( n ) − x ( n − 1) + x ( n − 2).
2 6 2 3
The system transfer function is
1 − 12 z−1 + 13 z−2
H (z) = H1 (z) H2 (z) = .
1 − 12 z−2 + 16 z−3
1 1
H (z) = −
=
1 + A1 z 1 1 − z p1 z−1
272 Realization of Discrete Systems
the error in coefficient A1 is the same as the error in the system pole z p1 . If
the coefficient is quantized with a step ∆ then the error in the pole location
is of order ∆. The same holds for the system zeros.
For a second-order system with real-valued coefficients and a pair of
complex-conjugated poles
1 1
H (z) = =
1 + A 1 z −1 + A 2 z −2 (1 − z p1 z−1 )(1 − z p2 z−1 )
the relation between the coefficients and the real and imaginary parts of the
poles z p1/2 = x p ± jy p is
1
H (z) =
1 − 2x p z −1 + ( x2p + y2p )z−2
A1 = −2x p
A2 = x2p + y2p .
The error in coefficient A1 defines the error in the real part of poles x p .
When the coefficient A2 assumes discrete values A2 = m∆, with A1 ∼
x p = n∆ then the imaginary part of poles may assume the values y p =
F √
± A2 − x2p = ± m∆ − n2 ∆2 with n2 ≤ mN. For small n, i.e., for small real
√
part of a pole, y p = ± ∆m. For N discretization levels, assuming that the
poles are within the unit circle x2p + y2p ≤ 1, the first discretization step is
√
changed from 1/N order to 1/ N order. The error, in this case, could be
significantly increased. The changes in y p due to the discretization of A2
may be large.
The quantization of x p and y p as a result of quantization of − A1 /2
and A2 = x2p + y2p is shown in Fig.6.7 for the case of N = 16 and N = 32
quantization levels. We see that the error in y p , when it assumes small
values, can be very large. We can conclude that the poles close to the unit
circle with larger imaginary values y p are less sensitive to the errors. The
highest error could appear if a second order real-valued pole (with y p = 0)
were implemented by using a second order system.
We have concluded that the poles close to the real axis (small y p )
are sensitive to the error in coefficients even in the second order systems.
The sensitivity increases with the system order, since the higher powers in
polynomial increase the maximal possible error.
Consider a general form of a polynomial in the transfer function,
written in two forms
y =Im{z } y =Im{z }
p p p p
1 1
0.5 0.5
0 0
-0.5 -0.5
-1 -1
-1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1
x =Re{z } x =Re{z }
p p p p
Figure 6.7 Quantization of the real and imaginary parts x p = Re{z p } and y p = Im{z p } of
poles (zeros) as a result of the quantization in 16 levels (left) and 32 levels (right) of the
coefficients A1 = −2x p and A2 = x2p + y2p .
and
P(z) = (z − z1 )(z − z2 )...(z − z M ).
If the coefficients A1 , A2 , ..., A M are changed for small ∆A1 , ∆A2 , ..., ∆A M
(due to quantization) then the pole position (without loss of generality and
for notation simplicity consider the pole z1 ) is changed for
- .
∼ ∂z1 ∂z1 ∂z1
∆z1 = ∆A1 + ∆A2 + ... + ∆A M . (6.4)
∂A1 ∂A2 ∂A M | z = z1
Since there is no a direct relation between z1 and A1 we will find ∂z1 /∂Ai
using
∂P(z) ∂P(z) ∂z1
= .
∂Ai |z=z1 ∂z1 ∂Ai |z=z1
From this relation it follows
∂P(z)
∂z1 ∂Ai |z=z1 z1M−i
= ∂P(z)
= .
∂Ai |z=z1 −(z1 − z2 )(z1 − z3 )...(z1 − z M )
∂z1 |z=z1
The coefficients ∂z1 /∂Ai|z=z1 could be large, especially in the case when
there are close poles, with a small distance (zi − zk ).
274 Realization of Discrete Systems
with
12 7 111 95
P(z) = (z − )(z − )(z − )(z − )
27 29 132 101
∼
= (z − 0.4444)(z − 0.2414)(z − 0.8409)(z − 0.9406)
In the realization of this system the coefficients are rounded to two decimal
positions, with absolute error up to 0.005. Find the poles of the system with
rounded coefficients.
⋆The system denominator is
P(z) ∼
= z4 − 2.4673z3 + 2.1200z2 − 0.7336z + 0.0849.
with poles
The poles of the function with rounded coefficients can differ significantly
from the original pole values. Maximal error in poles is 0.115. One pole is on
the unit circle making the system with rounded coefficients unstable, in this
case.
Note that if the system is written as a product of the first-order func-
tions in the denominator and each pole value is rounded to two decimals
1
H (z) = 7 12 111 95
(z − 29 )( z − 27 )( z − 132 )( z − 101 )
P(z) ∼
= (z − 0.24)(z − 0.44)(z − 0.84)(z − 0.94)
the poles will differ from the original ones for no more than 0.005.
If the poles are grouped into the second-order terms (what should be
done if the coefficients were complex-conjugate in order to avoid calculation
with complex valued coefficients), then
P(z) ∼
= (z − 0.6858z + 0.1073)(z − 1.7815z + 0.7910).
P(z)
we will get
The sensitivity analysis for this example can be done for each pole.
Assume that the poles are denoted as z1 = 12/27, z2 = 7/29, z3 = 111/132
and z4 = 95/101. Then
∆z1 ∼
= 0.0878.
The true error is ∆z1 = 0.0926. A small difference is due to the linear approx-
imation, assuming small ∆Ai . The obtained result is a good estimate of an
order of error for the pole z1 . The error in z1 is about 18.5 time greater than
the maximal error in the coefficients Ai , that is of order 0.005.
Commonly real-valued signals are processed and the poles and zeros in the
transfer function are in complex-conjugated pairs. In that case it is better to
group these pairs into second order systems to avoid complex calculations.
The transfer function is of the form
B00 + B10 z−1 + B20 z−2 B0K + B1K z−1 + B2K z−2
H (z) = × ... ×
1 − A10 z−1 − A20 z−2 1 − A1K z−1 − A2K z−2
= H0 (z) H1 (z)...HK (z),
where
B0i + B1i z−1 + B2i z−2
Hi (z) =
1 − A1i z−1 − A2i z−2
are second-order systems with real-valued coefficients. The whole system
may be realized as a cascade of lower-order (first or second-order) systems,
Fig.6.9. Of course, if there are some real-valued poles then there is no need
to group them. It is better to keep the realization order of the subsystems as
low as possible.
Ljubiša Stanković Digital Signal Processing 277
z-1 z-1
+ + + +
A B A B
10 -1 10 1K -1 1K
z z
A B A B
20 20 2K 2K
H(z)
Y ( z ) = H ( z ) R ( z ) = H ( z ) X ( z ) − H 2 ( z )Y ( z ) .
Y (z) H (z)
He (z) = = .
X (z) 1 + H 2 (z)
278 Realization of Discrete Systems
H2(z)
x(n) y(n)
+ H(z) +
- +
z -1
xp
H(z)
Figure 6.11 Complete second-order subsystem with complex-conjugate pair of poles realized
using the first-order systems.
y pL z−1
Qi ( z ) =
1 − 2x pL z−1 + x2pL z−2 + y2pL z−2
y pL z−1
=
(1 − x pL z−1 )2 + y2pL z−2
1 1
= y pL z−1 * +2
(1 − x pL z−1 )2 y pL z−1
1+ 1− x pL z−1
H (z) H2 (z)
=
1 + H 2 (z)
where
y pL z−1 1
H (z) = −
and H2 (z) = .
1 − x pL z 1 1 − x pL z−1
H(z)
x(n) yp y(n)
-1
+ z
+
xp
Figure 6.12 First-order system for the realization of the second-order system with complex-
conjugate pair of poles.
The error in one coefficient (real or imaginary part of a pole) does not
influence the other coefficients. However if an error in the signal calculation
happens in one cascade, then it will propagate as an input to the following
cascades. In that sense it would be the best to order cascades in such a way
that the lowest probability of an error appears in early cascades. From the
analysis of error we can conclude that the cascades with the poles and zeros
close to the origin are more sensitive to the error and should be used in later
stages.
1.4533(1 + z−1 )3
H (z) =
(−0.8673z−1 + 3.1327)(3.0177z−2 − 5.434z−1 + 7.54)
1 + z −1 1 + 2z−1 + z−2
= 0.0615 ×
1 − 0.2769z−1 1 − 0.7207z−1 + 0.4002z−2
⋆(a) Realization of the system H (z) when both the first and the
second-order subsystems can used is done according to system transfer func-
tion as in Fig.6.13.
(b) For the first-order systems the realization should be done based on
1 + z −1
H (z) = 0.0615 × (1 + z −1 ) × (1 + z −1 )
1 − 0.2769z−1
1
× ,
1 − 0.7207z−1 + 0.4002z−2
280 Realization of Discrete Systems
+ +
1 0.2769 0.7207 -1 2
z
-0.4002 1
with
1
1 − 0.7207z−1 + 0.4002z−2
1
=
(1 − (0.3603 + j0.5199)z−1 )(1 − (0.3603 − j0.5199)z−1 )
1
=
1 − 2 × 0.3603z−1 + 0.36032 z−2 + 0.51992 z−2
1 1 1
= = .
0.51992 z−2 + (1 − 0.3603z−1 )2 (1 − 0.3603z−1 )2 1 + ( 0.5199z−1−1 )2
1−0.3603z
In this way the system can be written and realized in terms of the first-order
subsystems
1 + z −1 1 + z −1
H (z) = 0.0615 ×
1 − 0.2769z−1 1 − 0.3603z−1
1+z − 1 1
× × −1 .
1 − 0.3603z−1 −1
1 + 0.5199z −1 × 0.5199z −1
1−0.3603z 1−0.3603z
x(n) 0.0615
+ + + +
-1 -1
z z
0.2769 0.3603
y(n)
+ + +
-1
z yp=0.5199
y yp
p
-1 + -1 +
z z
0.3603
0.3603 0.3603
In the case of a parallel realization the error in one subsystem does not
influence the other subsystems. If an error in the signal calculation appears
in one parallel subsystem, then it will influence the output signal, but will
not influence the outputs of other parallel subsystems.
x(n) y(n)
+ + +
B00
z-1
+ +
A10 B10
-1
z
A B
20 20
+ +
B01
z-1
+ +
A B
11 11
z-1
A21 B21
+ +
B0K
z-1
+ +
A1K B1K
-1
z
A2K B2K
z-1
1.1078 -1 0.2542
z
-0.5482
0.7256
+ +
-1
z
0.9246 -0.084
z-1
-0.2343
+ + +
1.1078 1 0.9246 0.0858
z-1 z
-1
x(n) B y(n)
0
+ +
z-1 z-1
+ +
A B
1 -1 -1 1
z z
+ +
A B
2 2
-1 -1
z z
AN BM
transfer functions it follows that the inverse realization has the same transfer
function as the original realization.
y(n) B x(n)
0
+ +
z-1 z-1
+ +
A B
1 1
z-1 z-1
+ +
A B
2 2
-1 -1
z z
A B
N M
in the pervious section. Systems that would not have recursions, when the
output signal is a linear combination of the input signal and its delayed
versions only,
are the FIR systems. These systems are always stable. The FIR systems can
also have a linear phase.
Im{ H (e jω }
arg{ H (e jω } = arctan{ } = −ωq (6.5)
Re{ H (e jω }
is also acceptable in these systems. They will have a constant group delay
d(arg{ H (e jω })
τg = − =q
dω
286 Realization of Discrete Systems
and will not distort the impulse response with respect to the zero-phase
system. The impulse response will only be delayed in time for q.
Example 6.6. Consider an input signal of the form
M
x (n) = ∑ A m e j ( ωm n + θ m ).
m =1
arg{ H (e jω }
τϕ = − =q
ω
and the group delay τg are the same. In general, the group delay and the
phase delays are different. The group delay, as notion dual to the instanta-
neous frequency, is introduced and discussed in the first chapter.
N −1 N −1 N −1
H (e jω ) = ∑ h(n)e− jωn = ∑ h(n) cos(ωn) − j ∑ h(n) sin(ωn). (6.6)
n =0 n =0 n =0
Combining the linear phase condition (6.5) with form (6.6), we get
or
N −1
∑ h(n)[sin(ωq) cos(ωn) − cos(ωq)sin (ωn)] = 0.
n =0
N −1
∑ h(n) sin(ω (n − q)) = 0. (6.7)
n =0
N−1
q=
2
h(n) = h( N − 1 − n), 0 ≤ n ≤ N − 1.
Since the Fourier transform is unique, this is the unique solution for the
linear phase condition. It is illustrated for an even and odd N in Fig.6.20.
From the symmetry condition it is easy to conclude that there is no a causal
linear phase system with infinite impulse response.
6.2.2 Windows
When a system obtained from the design procedure is an IIR system and
the requirement is to implement it as an FIR system, in order to get a linear
phase or to guaranty the system stability (when small changes of coefficients
are possible), then the most obvious way is to truncate the desired impulse
response hd (n) of the resulting IIR system. The impulse response of the FIR
system is
!
hd (n) for 0 ≤ n ≤ N − 1
h(n) =
0 elsewhere.
This form can be written as
h ( n ) = h d ( n ) w ( n ),
where !
1 for 0≤n≤ N−1
w(n) =
0 elsewhere
288 Realization of Discrete Systems
q=16
h(n)
N=32
0 16 32
n
q=16.5
h(n)
N=33
0 16.5 33
n
Figure 6.20 Impulse response of a system with a linear phase for an even and odd N.
is the rectangular window function. In the Fourier domain the desired im-
pulse response truncation by a window function will mean a convolution of
the desired frequency response with the frequency response of the window
function
H (e jω ) = Hd (e jω ) ∗ W (e jω ).
Without loss of generality, assume that the most significant values of hd (n)
are within − N/2 ≤ n ≤ N/2 − 1. The impulse response hc (n) can assume
nonzero values only within − N/2 ≤ n ≤ N/2 − 1. Therefore,
N/2−1 − N/2−1 ∞
e2 = ∑ |hd (n) − hc (n)|2 + ∑ |hd (n)|2 + ∑ |hd (n)|2 .
n=− N/2 n=−∞ n= N/2
Since the last two terms are hc (n) independent and all three terms are
non negative, the error e2 is minimal if