##### Document Text Contents

Page 1

INTRODUCTION TO

Signal Processing

Solutions Manual

Sophocles J. Orfanidis

Department of Electrical & Computer Engineering

Rutgers University, Piscataway, NJ 08855

[email protected]

Copyright © 2010 by Sophocles J. Orfanidis

Web page: www.ece.rutgers.edu/~orfanidi/i2sp

Chapter 1 Problems

Problem 1.1

The Nyquist interval is [�fs/2, fs/2]= [�4,4] Hz. The 6 Hz frequency of the wheel lies outside

it, therefore, it will be aliased with f � fs = 6 � 8 = �2 Hz. Thus, the wheel will appear to be

turning at 2 Hz in the opposite direction.

If fs = 12, the Nyquist interval is [�6,6]. Points on the wheel will appear to be moving up and

down, with either positive or negative sense of rotation.

For the other two sampling frequencies, the Nyquist interval is [�8,8] or [�12,12] Hz, and

therefore, the original 6 Hz frequency lies in it and no aliasing will be perceived.

Problem 1.2

The three terms of the signal correspond to frequencies f1 = 1, f2 = 4, and f3 = 6 Hz. Of

these, f2 and f3 lie outside the Nyquist interval [�2.5,2.5]. Therefore, they will be aliased with

f2 � fs = 4� 5 = �1 and f3 � fs = 6� 5 = 1 and the aliased signal will be:

xa(t)= 10 sin(2πt)+10 sin(2π(�1)t)+5 sin(2πt)= 5 sin(2πt)

To show that they have the same sample values, we set t = nT, with T = 1/fs = 1/5 sec. Then,

x(nT)= 10 sin(2πn/5)+10 sin(8πn/5)+5 sin(12πn/5)

But,

sin(8πn/5)= sin(2πn� 2πn/5)= � sin(2πn/5)

and

sin(12πn/5)= sin(2πn+ 2πn/5)= sin(2πn/5)

Thus,

x(nT) = 10 sin(2πn/5)�10 sin(2πn/5)+5 sin(2πn/5)

= 5 sin(2πn/5)= xa(nT).

If fs = 10 Hz, then the Nyquist interval is [�5,5] Hz and only f3 lies outside it. It will be aliased

with f3 � fs = 6� 10 = �4 resulting in the aliased signal:

xa(t)= 10 sin(2πt)+10 sin(8πt)+5 sin(2π(�4)t)= 10 sin(2πt)+5 sin(8πt)

Problem 1.3

Using the trig identity 2 sinα sinβ = cos(α� β)� cos(α+ β), we find:

x(t) = cos(5πt)+4 sin(2πt)sin(3πt)= cos(5πt)+2[cos(πt)� cos(5πt)]

= 2 cos(πt)� cos(5πt)

The frequencies in the signal are f1 = 0.5 and f2 = 2.5 kHz. The Nyquist interval is [�1.5,1.5]

kHz, and f2 lies outside it. Thus, it will be aliased with f2a = 2.5 � 3 = �0.5 giving rise to the

signal:

1

Page 2

xa(t)= 2 cos(2πf1t)− cos(2πf2at)= 2 cos(πt)− cos(−πt)= cos(πt)

A class of signals aliased with x(t) and xa(t) is obtained by replacing f1 and f2 by their shifted

versions: f1 +mfs, f2 + nfs resulting in:

xmn(t)= 2 cos(πt + 6πmt)− cos(πt − 6πnt)

Problem 1.4

Using a trig identity, we write x(t)= cos(8πt)+ cos(10πt)+ cos(2πt). The frequencies con-

tained in x(t) are thus, 1 Hz, 4 Hz, and 5 Hz. If the sampling rate is 5 Hz, then the Nyquist

interval is [−2.5,2.5] Hz, and therefore, the 4 Hz and 5 Hz components lie outside it and will be

aliased with the frequencies 4− 5 = −1 Hz and 5− 5 = 0 Hz, etc.

Problem 1.5

Using trig identities we find:

x(t)= sin(6πt)[1+ 2 cos(4πt)] = sin(6πt)+ sin(10πt)+ sin(2πt)

with frequency content: f1 = 3, f2 = 5, f3 = 1 kHz. The Nyquist interval is [−2,2] kHz, and the

aliased frequencies are:

f1a = f1 − fs = 3− 4 = −1, f2a = f2 − fs = 5− 4 = 1, f3a = f3 = 1

Thus,

xa(t)= sin(−2πt)+ sin(2πt)+ sin(2πt)= sin(2πt)

Problem 1.6

Use the trigonometric identity 2 cosa cosb = cos(a+ b)+ cos(a− b) three times to get

x(t) = 2[cos(10πt)+ cos(6πt)]cos(12πt)

= cos(22πt)+ cos(2πt)+ cos(18πt)+ cos(6πt)

The frequencies present in this signal are f1 = 11, f2 = 1, f3 = 9, and f4 = 3 Hz. With a sampling

rate of 10 Hz, only f1 and f3 lie outside the Nyquist interval [−5,5] Hz, and they will be aliased

with f1 − fs = 11− 10 = 1 Hz and f3 − fs = 9− 10 = −1 Hz. The aliased signal will be:

xa(t) = cos(2π(1)t)+ cos(2πt)+ cos(2π(−1)t)+ cos(6πt)

= 3 cos(2πt)+ cos(6πt)

To prove the equality of the samples, replace t = nT = n/10, because T = 1/fs = 1/10. Then,

x(nT)= cos(22πn/10)+ cos(2πn/10)+ cos(18πn/10)+ cos(6πn/10)

But, cos(22πn/10)= cos(2πn+2πn/10)= cos(2πn/10) and similarly, cos(18πn/10)= cos(2πn−

2πn/10)= cos(2πn/10). Therefore, the sample values become

2

x(nT) = cos(2πn/10)+ cos(2πn/10)+ cos(2πn/10)+ cos(6πn/10)

= 3 cos(2πn/10)+ cos(6πn/10)= xa(nT)

If fs = 12 Hz, then f1 and f3 will lie outside of [−6,6] and will be aliased with 11− 12 = −1 and

9− 12 = −3. The aliased signal will be:

xa(t) = cos(2π(−1)t)+ cos(2πt)+ cos(2π(−3)t)+ cos(6πt)

= 2 cos(2πt)+2 cos(6πt)

Problem 1.7

We use the same technique as in the square-wave Example 1.4.6. At a sampling rate of 8 Hz, the

signal frequencies of

{1,3,5,7,9,11,13,15, . . . }

will be aliased with:

{1,3,−3,−1,1,3,−3,−1, . . . }

Therefore only sin(2πt) and sin(6πt) terms will appear in the aliased signal. Thus, we write it

in the form:

xa(t)= B sin(2πt)+C sin(6πt)

To determine B and C, we demand that xa(t) and x(t) agree at the sampling instants t = nT =

n/8, because T = 1/fs = 1/8 sec. Therefore, we demand

B sin(2πn/8)+C sin(6πn/8)= x(n/8)

Setting n = 1,2, we get two equations

B sin(2π/8)+C sin(6π/8)= x(1/8)= 0.5

B sin(4π/8)+C sin(12π/8)= x(2/8)= 1

⇒

B

1√

2

+C 1√

2

= 0.5

B−C = 1

The values for x(1/8) and x(2/8) were obtained by inspecting the triangular waveform. Solving

for B and C, we find:

B =

√

2+ 2

4

, C =

√

2− 2

4

Problem 1.8

For fs = 5 kHz, we have:

xa(t)= sin(2πf1t)

For fs = 10 kHz, we have:

xa(t)= 2 sin(2πf1t)+ sin(2πf2t)

3

Page 74

Eliminating F(z) and E(z) in favor of X(z) and Y(z) gives:

Y(z)= z−1H1(z)E(z)= z−1H1(z)

(

X(z)+F(z)) = z−1H1(z)(X(z)+H2(z)Y(z))

and (

1− z−1H1(z)H2(z)

)

Y(z)= z−1H1(z)X(z)

Solving for the closed-loop transfer function H(z)= Y(z)/X(z), we find:

H(z)= z

−1H1(z)

1− z−1H1(z)H2(z)

Case B has the same transfer function H(z). Cases C and D have:

H(z)= H1(z)

1− z−1H1(z)H2(z)

Problem 7.29

The overall transfer function of this system can be computed using the results of Problem 7.28

and the individual transfer functions:

H1(z)=

1+ z−1

1− 0.5z−1 , H2(z)= z

−1 0.4(1− z−1)

1− 0.4z−1

H(z)= H1(z)

1−H1(z)H2(z)

=

1+ z−1

1− 0.5z−1

1− 1+ z

−1

1− 0.5z−1 ·

0.4(1− z−1)z−1

1− 0.4z−1

which gives

H(z)= 1+ 0.6z

−1 − 0.4z−2

1− 1.3z−1 + 0.2z−2 + 0.4z−3

In case A, we assign internal state variables as shown in Fig. P7.33.

The difference equations and sample processing algorithm, written in the right computational

order, are:

e(n) = x(n)+f(n− 1)

w(n) = 0.5w(n− 1)+e(n)

y(n) = w(n)+w(n− 1)

v(n) = 0.4v(n− 1)+y(n)

f(n) = 0.4v(n)−0.4v(n− 1)

for each input sample x do:

e = x+ f1

w0 = 0.5w1 + e

y = w0 +w1

v0 = 0.4v1 + y

f0 = 0.4(v0 − v1)

f1 = f0 (update delays)

w1 = w0

v1 = v0

Notice how the extra delay factor z−1 inH2(z) delays the output f(n) and allows the start of the

computational cycle with e(n)= x(n)+f(n−1). The output f(n) is computed last. If this delay

were not there, we would have e(n)= x(n)+f(n), which is not possible because f(n) is not yet

146

x(n)

e(n) w(n)

w0

f0

v0

w1

v1

f1 v(n)

w(n-1)

v(n-1)

f(n)f(n-1)

y(n)

0.5

z-1

0.4

z-1

z-1

-0.4

0.4

Fig. P7.33 Case A of Problem 7.29.

computed.

In case B, shown in Fig. P7.34, the overall transfer function remains the same. The change causes

a re-ordering of the computational cycle. The difference equations and corresponding sample

processing algorithm are now:

v(n) = 0.4v(n− 1)+y(n− 1)

f(n) = 0.4v(n)−0.4v(n− 1)

e(n) = x(n)+f(n)

w(n) = 0.5w(n− 1)+e(n)

y(n) = w(n)+w(n− 1)

for each input sample x do:

v0 = 0.4v1 + y1

f = 0.4(v0 − v1)

e = x+ f

w0 = 0.5w1 + e

y = w0 +w1

y1 = y (update delays)

w1 = w0

v1 = v0

x(n)

e(n)

f(n)

w(n)

w0

v0

w1

v1

y1

v(n)

w(n-1)

v(n-1)

y(n-1)

y(n)

0.5

z-1

0.4

z-1

-0.4

0.4

z-1

Fig. P7.34 Case B of Problem 7.29.

Now, the start of the computational cycle is v(n)= 0.4v(n−1)+y(n−1), where both quantities

147

Page 75

v(n − 1) and y(n − 1) are available from the previous time step. Note that y(n) is computed

last and then put into the feedback delay to be used in the next time step.

Problem 7.30

The filter G(z) can be implemented by the input/output sample processing algorithm of its

canonical form:

for each x do:

y = 0.2x+ 0.8w1

w1 = y

where w1 is its internal state, that is, w1(n)= y(n − 1) in the difference equation: y(n)=

0.2x(n)+0.8y(n − 1). Fig. P7.35 shows the three cases with explicitly assigned signal labels.

Identifying the proper input and output to the filterG(z) and using the above sample processing

algorithm, we obtain the implementation of the three cases:

for each x do:

u = 0.2v4 + 0.8w1

w1 = u

y = v0 = x+ u

delay(4,v)

for each x do:

y = x+ v4

v0 = 0.2y + 0.8w1

w1 = v0

delay(4,v)

for each x do:

y = u0 = x+ v2

v0 = 0.2u2 + 0.8w1

w1 = v0

delay(2,v), delay(2,u)

z-4

G(z)

z-4

G(z)

z-2 z

-2

G(z)

(a) (b) (c)

x

u

v0

v4

y x x

v0 v0

u0

u2

v4 v2

y y

Fig. P7.35 Feedback systems of Problem 7.30.

148

Chapter 8 Problems

Problem 8.1

The transfer function is obtained by the period-3 repetition of the sequence b(n)= [0,1,2], that

is, in the time- and z-domains:

h(n) = b(n)+b(n− 3)+b(n− 6)+· · ·

H(z) = B(z)+z−3B(z)+z−6B(z)+· · · = B(z)

1− z−3 =

z−1 + 2z−2

1− z−3

The direct and canonical realizations are shown in Fig. P8.1. Their sample processing algorithms

are:

for each input x do:

v0 = x

y = w0 = v1 + 2v2 +w3

delay(2,v)

delay(3,w)

for each input x do:

w0 = x+w3

y = w1 + 2w2

delay(3,w)

x y x y

z

-1

z

-1

z

-1

z

-1

z

-1

z

-1

w0 w0

w1 w1

w2 w2

w3 w3

z

-1

z

-1

v0

v1

v2

2 2

Fig. P8.1 Direct and canonical realizations of Problem 8.1.

For the direct form realization, if the input is a unit impulse so that x = 1 at the first iteration of

the algorithm and zero for all other iterations, then the first three iterations will fill the feedback

delay buffer w with the desired waveform samples. From then on, the feed-forward part of the

algorithm will be zero, giving rise to the following generation algorithm, expressed in its of the

linear and circular buffer versions:

repeat forever:

w0 = w3

delay(3,w)

repeat forever:

∗p = tap(3,w, p,3)

cdelay(3,w,&p)

The following table shows the initialization and steady part of the linear buffer case. The columns

v0, v1, v2 are the input impulse and its delays: δ(n), δ(n−1), δ(n−2). At each iteration, w0 is

replaced by the value of w3 and only then is the row of w’s shifted to the right to give the values

of w1, w2, w3 of the next row. The w-columns are the successively delayed output signals y(n),

y(n− 1), y(n− 2), y(n− 3):

149

Page 147

The reason for using N = 2000 was because we wanted to evaluate the sample autocorrelation

at M = 150 lags. Thus, the number of lags is 150/2000 ≡ 7.5% of the block length, which falls

within the recommended range for getting statistically reliable lags using the routine corr, as

was mentioned in Section A.1.

Problem B.8

We look at the cases ωmin = 0.001π and ωmax = 0.1π, 0.3π, 0.6π. The designed model param-

eters are:

For ωmax = 0.1π, we have B = 3 and

c = 10, a = [0.6858,0.9686,0.9969]

For ωmax = 0.3π, we have B = 4 and

c = 6.6943, a = [0.0575,0.8592,0.9790,0.9969]

For ωmax = 0.6π, we have B = 6 and

c = 3.5944, a = [−0.8850,0.4756,0.8541,0.9594,0.9887,0.9969]

The spectra are shown in Figs. P14.6—P14.8.

Fig. P14.6 Model spectrum and periodogram estimate. ωmax = 0.1π.

Problem B.9

Solving the equation (0.99574)c= 0.94791, we find c = 12.5309. Then, the filter of Eq. (B.30)

becomes:

H(z)= G(1− 0.98444z

−1)(1− 0.82159z−1)(1− 0.08522z−1)

(1− 0.99574z−1)(1− 0.94791z−1)(1− 0.51153z−1)

The NRR = hTh of the filter (B.29) without the factor G can be computed analytically or much

faster using MATLAB to compute a fairly long portion of the impulse response h (e.g., of length

2000) and then setting G = 1/√NRR = 1/

√

hTh. The magnitude response squared of Eqs. (B.29)

and (B.30) are shown in Fig. P14.9. The cascade realization is shown in Fig. P14.10. The corre-

sponding sample processing algorithm will be:

292

Fig. P14.7 Model spectrum and periodogram estimate. ωmax = 0.3π.

Fig. P14.8 Model spectrum and periodogram estimate. ωmax = 0.6π.

for each white noise input sample x do:

u0 = 0.57534x+ 0.99574u1

x1 = u0 − 0.98444u1

u1 = u0

v0 = x1 + 0.94791v1

x2 = v0 − 0.83392v1

v1 = v0

w0 = x2 + 0.53568w1

y = w0 − 0.07568w1

w1 = w0

293

Page 148

Fig. P14.9 1/f noise filters for Problem B.9.

z-1

x

u0 v0 w0

u1

x1 x2

v1 w1

y

z-1 z-1

0.57534

0.99574 -0.98444 0.94791 -0.83392 0.53568 -0.07568

Fig. P14.10 Cascade realization of 1/f noise filter.

294

INTRODUCTION TO

Signal Processing

Solutions Manual

Sophocles J. Orfanidis

Department of Electrical & Computer Engineering

Rutgers University, Piscataway, NJ 08855

[email protected]

Copyright © 2010 by Sophocles J. Orfanidis

Web page: www.ece.rutgers.edu/~orfanidi/i2sp

Chapter 1 Problems

Problem 1.1

The Nyquist interval is [�fs/2, fs/2]= [�4,4] Hz. The 6 Hz frequency of the wheel lies outside

it, therefore, it will be aliased with f � fs = 6 � 8 = �2 Hz. Thus, the wheel will appear to be

turning at 2 Hz in the opposite direction.

If fs = 12, the Nyquist interval is [�6,6]. Points on the wheel will appear to be moving up and

down, with either positive or negative sense of rotation.

For the other two sampling frequencies, the Nyquist interval is [�8,8] or [�12,12] Hz, and

therefore, the original 6 Hz frequency lies in it and no aliasing will be perceived.

Problem 1.2

The three terms of the signal correspond to frequencies f1 = 1, f2 = 4, and f3 = 6 Hz. Of

these, f2 and f3 lie outside the Nyquist interval [�2.5,2.5]. Therefore, they will be aliased with

f2 � fs = 4� 5 = �1 and f3 � fs = 6� 5 = 1 and the aliased signal will be:

xa(t)= 10 sin(2πt)+10 sin(2π(�1)t)+5 sin(2πt)= 5 sin(2πt)

To show that they have the same sample values, we set t = nT, with T = 1/fs = 1/5 sec. Then,

x(nT)= 10 sin(2πn/5)+10 sin(8πn/5)+5 sin(12πn/5)

But,

sin(8πn/5)= sin(2πn� 2πn/5)= � sin(2πn/5)

and

sin(12πn/5)= sin(2πn+ 2πn/5)= sin(2πn/5)

Thus,

x(nT) = 10 sin(2πn/5)�10 sin(2πn/5)+5 sin(2πn/5)

= 5 sin(2πn/5)= xa(nT).

If fs = 10 Hz, then the Nyquist interval is [�5,5] Hz and only f3 lies outside it. It will be aliased

with f3 � fs = 6� 10 = �4 resulting in the aliased signal:

xa(t)= 10 sin(2πt)+10 sin(8πt)+5 sin(2π(�4)t)= 10 sin(2πt)+5 sin(8πt)

Problem 1.3

Using the trig identity 2 sinα sinβ = cos(α� β)� cos(α+ β), we find:

x(t) = cos(5πt)+4 sin(2πt)sin(3πt)= cos(5πt)+2[cos(πt)� cos(5πt)]

= 2 cos(πt)� cos(5πt)

The frequencies in the signal are f1 = 0.5 and f2 = 2.5 kHz. The Nyquist interval is [�1.5,1.5]

kHz, and f2 lies outside it. Thus, it will be aliased with f2a = 2.5 � 3 = �0.5 giving rise to the

signal:

1

Page 2

xa(t)= 2 cos(2πf1t)− cos(2πf2at)= 2 cos(πt)− cos(−πt)= cos(πt)

A class of signals aliased with x(t) and xa(t) is obtained by replacing f1 and f2 by their shifted

versions: f1 +mfs, f2 + nfs resulting in:

xmn(t)= 2 cos(πt + 6πmt)− cos(πt − 6πnt)

Problem 1.4

Using a trig identity, we write x(t)= cos(8πt)+ cos(10πt)+ cos(2πt). The frequencies con-

tained in x(t) are thus, 1 Hz, 4 Hz, and 5 Hz. If the sampling rate is 5 Hz, then the Nyquist

interval is [−2.5,2.5] Hz, and therefore, the 4 Hz and 5 Hz components lie outside it and will be

aliased with the frequencies 4− 5 = −1 Hz and 5− 5 = 0 Hz, etc.

Problem 1.5

Using trig identities we find:

x(t)= sin(6πt)[1+ 2 cos(4πt)] = sin(6πt)+ sin(10πt)+ sin(2πt)

with frequency content: f1 = 3, f2 = 5, f3 = 1 kHz. The Nyquist interval is [−2,2] kHz, and the

aliased frequencies are:

f1a = f1 − fs = 3− 4 = −1, f2a = f2 − fs = 5− 4 = 1, f3a = f3 = 1

Thus,

xa(t)= sin(−2πt)+ sin(2πt)+ sin(2πt)= sin(2πt)

Problem 1.6

Use the trigonometric identity 2 cosa cosb = cos(a+ b)+ cos(a− b) three times to get

x(t) = 2[cos(10πt)+ cos(6πt)]cos(12πt)

= cos(22πt)+ cos(2πt)+ cos(18πt)+ cos(6πt)

The frequencies present in this signal are f1 = 11, f2 = 1, f3 = 9, and f4 = 3 Hz. With a sampling

rate of 10 Hz, only f1 and f3 lie outside the Nyquist interval [−5,5] Hz, and they will be aliased

with f1 − fs = 11− 10 = 1 Hz and f3 − fs = 9− 10 = −1 Hz. The aliased signal will be:

xa(t) = cos(2π(1)t)+ cos(2πt)+ cos(2π(−1)t)+ cos(6πt)

= 3 cos(2πt)+ cos(6πt)

To prove the equality of the samples, replace t = nT = n/10, because T = 1/fs = 1/10. Then,

x(nT)= cos(22πn/10)+ cos(2πn/10)+ cos(18πn/10)+ cos(6πn/10)

But, cos(22πn/10)= cos(2πn+2πn/10)= cos(2πn/10) and similarly, cos(18πn/10)= cos(2πn−

2πn/10)= cos(2πn/10). Therefore, the sample values become

2

x(nT) = cos(2πn/10)+ cos(2πn/10)+ cos(2πn/10)+ cos(6πn/10)

= 3 cos(2πn/10)+ cos(6πn/10)= xa(nT)

If fs = 12 Hz, then f1 and f3 will lie outside of [−6,6] and will be aliased with 11− 12 = −1 and

9− 12 = −3. The aliased signal will be:

xa(t) = cos(2π(−1)t)+ cos(2πt)+ cos(2π(−3)t)+ cos(6πt)

= 2 cos(2πt)+2 cos(6πt)

Problem 1.7

We use the same technique as in the square-wave Example 1.4.6. At a sampling rate of 8 Hz, the

signal frequencies of

{1,3,5,7,9,11,13,15, . . . }

will be aliased with:

{1,3,−3,−1,1,3,−3,−1, . . . }

Therefore only sin(2πt) and sin(6πt) terms will appear in the aliased signal. Thus, we write it

in the form:

xa(t)= B sin(2πt)+C sin(6πt)

To determine B and C, we demand that xa(t) and x(t) agree at the sampling instants t = nT =

n/8, because T = 1/fs = 1/8 sec. Therefore, we demand

B sin(2πn/8)+C sin(6πn/8)= x(n/8)

Setting n = 1,2, we get two equations

B sin(2π/8)+C sin(6π/8)= x(1/8)= 0.5

B sin(4π/8)+C sin(12π/8)= x(2/8)= 1

⇒

B

1√

2

+C 1√

2

= 0.5

B−C = 1

The values for x(1/8) and x(2/8) were obtained by inspecting the triangular waveform. Solving

for B and C, we find:

B =

√

2+ 2

4

, C =

√

2− 2

4

Problem 1.8

For fs = 5 kHz, we have:

xa(t)= sin(2πf1t)

For fs = 10 kHz, we have:

xa(t)= 2 sin(2πf1t)+ sin(2πf2t)

3

Page 74

Eliminating F(z) and E(z) in favor of X(z) and Y(z) gives:

Y(z)= z−1H1(z)E(z)= z−1H1(z)

(

X(z)+F(z)) = z−1H1(z)(X(z)+H2(z)Y(z))

and (

1− z−1H1(z)H2(z)

)

Y(z)= z−1H1(z)X(z)

Solving for the closed-loop transfer function H(z)= Y(z)/X(z), we find:

H(z)= z

−1H1(z)

1− z−1H1(z)H2(z)

Case B has the same transfer function H(z). Cases C and D have:

H(z)= H1(z)

1− z−1H1(z)H2(z)

Problem 7.29

The overall transfer function of this system can be computed using the results of Problem 7.28

and the individual transfer functions:

H1(z)=

1+ z−1

1− 0.5z−1 , H2(z)= z

−1 0.4(1− z−1)

1− 0.4z−1

H(z)= H1(z)

1−H1(z)H2(z)

=

1+ z−1

1− 0.5z−1

1− 1+ z

−1

1− 0.5z−1 ·

0.4(1− z−1)z−1

1− 0.4z−1

which gives

H(z)= 1+ 0.6z

−1 − 0.4z−2

1− 1.3z−1 + 0.2z−2 + 0.4z−3

In case A, we assign internal state variables as shown in Fig. P7.33.

The difference equations and sample processing algorithm, written in the right computational

order, are:

e(n) = x(n)+f(n− 1)

w(n) = 0.5w(n− 1)+e(n)

y(n) = w(n)+w(n− 1)

v(n) = 0.4v(n− 1)+y(n)

f(n) = 0.4v(n)−0.4v(n− 1)

for each input sample x do:

e = x+ f1

w0 = 0.5w1 + e

y = w0 +w1

v0 = 0.4v1 + y

f0 = 0.4(v0 − v1)

f1 = f0 (update delays)

w1 = w0

v1 = v0

Notice how the extra delay factor z−1 inH2(z) delays the output f(n) and allows the start of the

computational cycle with e(n)= x(n)+f(n−1). The output f(n) is computed last. If this delay

were not there, we would have e(n)= x(n)+f(n), which is not possible because f(n) is not yet

146

x(n)

e(n) w(n)

w0

f0

v0

w1

v1

f1 v(n)

w(n-1)

v(n-1)

f(n)f(n-1)

y(n)

0.5

z-1

0.4

z-1

z-1

-0.4

0.4

Fig. P7.33 Case A of Problem 7.29.

computed.

In case B, shown in Fig. P7.34, the overall transfer function remains the same. The change causes

a re-ordering of the computational cycle. The difference equations and corresponding sample

processing algorithm are now:

v(n) = 0.4v(n− 1)+y(n− 1)

f(n) = 0.4v(n)−0.4v(n− 1)

e(n) = x(n)+f(n)

w(n) = 0.5w(n− 1)+e(n)

y(n) = w(n)+w(n− 1)

for each input sample x do:

v0 = 0.4v1 + y1

f = 0.4(v0 − v1)

e = x+ f

w0 = 0.5w1 + e

y = w0 +w1

y1 = y (update delays)

w1 = w0

v1 = v0

x(n)

e(n)

f(n)

w(n)

w0

v0

w1

v1

y1

v(n)

w(n-1)

v(n-1)

y(n-1)

y(n)

0.5

z-1

0.4

z-1

-0.4

0.4

z-1

Fig. P7.34 Case B of Problem 7.29.

Now, the start of the computational cycle is v(n)= 0.4v(n−1)+y(n−1), where both quantities

147

Page 75

v(n − 1) and y(n − 1) are available from the previous time step. Note that y(n) is computed

last and then put into the feedback delay to be used in the next time step.

Problem 7.30

The filter G(z) can be implemented by the input/output sample processing algorithm of its

canonical form:

for each x do:

y = 0.2x+ 0.8w1

w1 = y

where w1 is its internal state, that is, w1(n)= y(n − 1) in the difference equation: y(n)=

0.2x(n)+0.8y(n − 1). Fig. P7.35 shows the three cases with explicitly assigned signal labels.

Identifying the proper input and output to the filterG(z) and using the above sample processing

algorithm, we obtain the implementation of the three cases:

for each x do:

u = 0.2v4 + 0.8w1

w1 = u

y = v0 = x+ u

delay(4,v)

for each x do:

y = x+ v4

v0 = 0.2y + 0.8w1

w1 = v0

delay(4,v)

for each x do:

y = u0 = x+ v2

v0 = 0.2u2 + 0.8w1

w1 = v0

delay(2,v), delay(2,u)

z-4

G(z)

z-4

G(z)

z-2 z

-2

G(z)

(a) (b) (c)

x

u

v0

v4

y x x

v0 v0

u0

u2

v4 v2

y y

Fig. P7.35 Feedback systems of Problem 7.30.

148

Chapter 8 Problems

Problem 8.1

The transfer function is obtained by the period-3 repetition of the sequence b(n)= [0,1,2], that

is, in the time- and z-domains:

h(n) = b(n)+b(n− 3)+b(n− 6)+· · ·

H(z) = B(z)+z−3B(z)+z−6B(z)+· · · = B(z)

1− z−3 =

z−1 + 2z−2

1− z−3

The direct and canonical realizations are shown in Fig. P8.1. Their sample processing algorithms

are:

for each input x do:

v0 = x

y = w0 = v1 + 2v2 +w3

delay(2,v)

delay(3,w)

for each input x do:

w0 = x+w3

y = w1 + 2w2

delay(3,w)

x y x y

z

-1

z

-1

z

-1

z

-1

z

-1

z

-1

w0 w0

w1 w1

w2 w2

w3 w3

z

-1

z

-1

v0

v1

v2

2 2

Fig. P8.1 Direct and canonical realizations of Problem 8.1.

For the direct form realization, if the input is a unit impulse so that x = 1 at the first iteration of

the algorithm and zero for all other iterations, then the first three iterations will fill the feedback

delay buffer w with the desired waveform samples. From then on, the feed-forward part of the

algorithm will be zero, giving rise to the following generation algorithm, expressed in its of the

linear and circular buffer versions:

repeat forever:

w0 = w3

delay(3,w)

repeat forever:

∗p = tap(3,w, p,3)

cdelay(3,w,&p)

The following table shows the initialization and steady part of the linear buffer case. The columns

v0, v1, v2 are the input impulse and its delays: δ(n), δ(n−1), δ(n−2). At each iteration, w0 is

replaced by the value of w3 and only then is the row of w’s shifted to the right to give the values

of w1, w2, w3 of the next row. The w-columns are the successively delayed output signals y(n),

y(n− 1), y(n− 2), y(n− 3):

149

Page 147

The reason for using N = 2000 was because we wanted to evaluate the sample autocorrelation

at M = 150 lags. Thus, the number of lags is 150/2000 ≡ 7.5% of the block length, which falls

within the recommended range for getting statistically reliable lags using the routine corr, as

was mentioned in Section A.1.

Problem B.8

We look at the cases ωmin = 0.001π and ωmax = 0.1π, 0.3π, 0.6π. The designed model param-

eters are:

For ωmax = 0.1π, we have B = 3 and

c = 10, a = [0.6858,0.9686,0.9969]

For ωmax = 0.3π, we have B = 4 and

c = 6.6943, a = [0.0575,0.8592,0.9790,0.9969]

For ωmax = 0.6π, we have B = 6 and

c = 3.5944, a = [−0.8850,0.4756,0.8541,0.9594,0.9887,0.9969]

The spectra are shown in Figs. P14.6—P14.8.

Fig. P14.6 Model spectrum and periodogram estimate. ωmax = 0.1π.

Problem B.9

Solving the equation (0.99574)c= 0.94791, we find c = 12.5309. Then, the filter of Eq. (B.30)

becomes:

H(z)= G(1− 0.98444z

−1)(1− 0.82159z−1)(1− 0.08522z−1)

(1− 0.99574z−1)(1− 0.94791z−1)(1− 0.51153z−1)

The NRR = hTh of the filter (B.29) without the factor G can be computed analytically or much

faster using MATLAB to compute a fairly long portion of the impulse response h (e.g., of length

2000) and then setting G = 1/√NRR = 1/

√

hTh. The magnitude response squared of Eqs. (B.29)

and (B.30) are shown in Fig. P14.9. The cascade realization is shown in Fig. P14.10. The corre-

sponding sample processing algorithm will be:

292

Fig. P14.7 Model spectrum and periodogram estimate. ωmax = 0.3π.

Fig. P14.8 Model spectrum and periodogram estimate. ωmax = 0.6π.

for each white noise input sample x do:

u0 = 0.57534x+ 0.99574u1

x1 = u0 − 0.98444u1

u1 = u0

v0 = x1 + 0.94791v1

x2 = v0 − 0.83392v1

v1 = v0

w0 = x2 + 0.53568w1

y = w0 − 0.07568w1

w1 = w0

293

Page 148

Fig. P14.9 1/f noise filters for Problem B.9.

z-1

x

u0 v0 w0

u1

x1 x2

v1 w1

y

z-1 z-1

0.57534

0.99574 -0.98444 0.94791 -0.83392 0.53568 -0.07568

Fig. P14.10 Cascade realization of 1/f noise filter.

294