Swiss Army 2004 May IEEEXplore 01296550, Semestr 3
[ Pobierz całość w formacie PDF ]
dsp
tips & tricks
Richard Lyons and Amy Bell
The Swiss Army Knife of Digital Networks
T
his article describes a
general discrete-sig-
nal network that
appears, in various
forms, inside many
digital signal processing (DSP)
applications. So the “DSP Tip” for
this column is for every DSP engi-
neer to become acquainted with this
network. Figure 1 shows how the
network’s structure has the distinct
look of a digital filter, a comb filter
followed by a second-order recur-
sive network. However, we do not
call this unique general network a
filter because its capabilities extend
far beyond simple filtering.
Through a series of examples, we
illustrate the fundamental strength
of the network: its ability to be
reconfigured to perform a surpris-
ingly large number of useful func-
tions based on the values of its
seven control parameters.
The general network has a trans-
fer function of
magnitude and phase responses.
The frequency axis in those tables is
normalized such that a value of 0.5
represents a frequency of
f
s
/
2,
where
f
s
is the sample rate in hertz.
c
1
z
−
N
H
(
z
)
=
(
1
−
)
b
1
z
−
1
b
2
z
−
2
b
0
+
+
×
a
2
z
−
2
.
(
Moving Averager
Referring to the first entry in Table
1, this network configuration is a
computationally efficient method
for computing the
N
-point moving
average of
x
1
/
a
0
−
a
1
z
−
1
−
1
)
From this point on, we’ll use DSP
filter lingo and call the second-order
recursive network a “biquad”
because its transfer function is the
ratio of two quadratic polynomials.
The tables in this article list various
signal-processing functions per-
formed by the network based on the
a
n
,
b
n
, and
c
1
coefficients. The vari-
able
N
is the order of the comb fil-
ter. Included in the tables are
depictions of the network’s impulse
response,
z
-plane pole/zero loca-
tions, as well as frequency-domain
. Also called a
recur-
sive running sum
or
boxcar
averager
, this structure is equivalent
to an
N
-tap direct convolution finite
impulse response (FIR) filter with
all the coefficients having a value of
1
(
n
)
However, this moving aver-
ager is efficient because it performs
only one add and one subtract per
output sample regardless of the
value of
N
(whereas an
N
-tap direct
convolution FIR filter must per-
form
N
/
N
.
1 additions per output
sample). The moving averager’s
transfer function is
H
ma
−
(
)
=
(
Comb
Second-Order Recursive Network (Biquad)
z
z
−
N
z
−
1
1
/
N
)(
1
−
)/(
1
−
).
+
+
x(n)
+
y(n)
+
–
“DSP Tips and Tricks” introduces prac-
tical tips and tricks of design and
implementation of signal processing
algorithms so that you may be able to
incorporate them into your designs.
We welcome readers who enjoy read-
ing this column to submit their contri-
butions. Contact Associate Editors Rick
Lyons (r.lyons@ieee.org) or Amy Bell
(abell@vt.edu).
a
0
b
0
z
–1
z
–N
a
1
b
1
z
–1
c
1
b
2
a
2
▲
1. General discrete-signal processing network.
90
IEEE SIGNAL PROCESSING MAGAZINE
MAY 2004
(
)(
+
1
)
’s numerator results in
N
zeros equally spaced around the
z
-
plane’s unit circle located at
z
H
ma
(
z
)
D
total
=
+
delay
samples where
delay
, typically in the range of
1
delay
delay
R
2
=
.
.(
3
)
−
0
.
5
(
delay
+
3
)(
delay
+
4
)
to 0.5, is a fraction of the 1
f
s
sam-
ple period. For example, when
delay
is 0.2, the network delay (at
low frequencies) is 1.2 samples. The
real-valued
R
coefficient is
/
=
e
j
2
π
k
/
N
, where integer
k
is
(
k
)
The band, centered at dc, over
which the group delay varies no
more than
0
’s denominator
places a single pole at
z
≤
k
<
N
.
H
ma
(
z
)
1 on the
unit circle, canceling the zero at
that location.
=
10 from the
specified
D
total
value, the bar in the
group delay plot, ranges roughly
from 0
|
delay
|
/
−
delay
delay
+
R
=
2
,
(
2
)
38
f
s
for this sec-
ond-order network. Performance
for
.
26
f
s
to 0
.
Differencer
This is a discrete version of a first-
order differentiator. An ideal differ-
entiator has a frequency magnitude
response that is a linear function of
frequency, and this network only
approaches that ideal at low fre-
quencies relative to
f
s
.
producing a
z
-plane transfer func-
tion of
H
1
,
del
(
delay
=
0
.
3
(
R
1
=−
0
.
182 and
z
−
1
z
)
=
(
R
+
)/
is shown in Table 1,
where we see the magnitude
response being constant.
The
flat
group delay band is
wider for negative
R
2
=
0
.
28
)
Rz
−
1
(
1
+
)
with a pole at
z
=−
R
and a zero at
z
=−
1
/
R
.
Performance for
delay
=
0.2
(
R
91) is shown in Table 1,
where we see the magnitude
response being constant. The band,
centered at dc, over which the
group delay varies no more than
|
delay
|
/
=
0
.
delay
than when
delay
is positive. For example, this
means if you desire a group delay of
D
total
=
Integrator
This structure performs the
run-
ning summation
of the
x
5 samples, it is better to
use an external unit delay and set
delay
to
2
.
inputs
samples, making it the discrete-
time equivalent of a continuous-
time integrator.
(
n
)
10 from the specified
D
total
value, the bar in the group delay
plot, ranges roughly from 0.1
f
s
to
0.2
f
s
for first-order networks. If
your signal is oversampled, making
it low in frequency relative to
f
s
,
this first-order all-pass delay net-
work may be of some use. If you
propose its use in a new design, you
can impress your colleagues by say-
ing this network is based on the
Thiran approximation [1].
5 rather than letting
delay
be 0.5. To ensure stability,
delay
must be greater than
−
0
.
1.
Reference [1] provides methods for
designing higher-order allpass delay
networks.
−
Leaky Integrator
This network configuration, also
called an
exponential averager
, is a
venerable structure used in low-pass
filter implementations for random
noise reduction. It is a first-order
infinite impulse response (IIR) filter
where, for stable low-pass operation,
the constant
Goertzel Network
Referring to the first entry in Table
2, this traditional Goertzel network
is used for single-tone detection
because it computes a single-bin
N
-
point discrete Fourier transform
(DFT) centered at an angle of
θ
=
Second-Order
Delay Network
A subclass of a second-order IIR fil-
ter, the coefficients in Table 1 yield
an all-pass network having a relative-
ly constant group at low frequencies
(over a wider frequency range, by
the way, than the first-order delay
network.) This network’s delay is
D
total
=
α
lies in the range
0
1.
This nonlinear-phase filter has a
single pole at
z
<α<
N
rad on the unit circle,
corresponding to a cyclic frequency
of
kf
s
/
2
π
k
/
on the
z
-
plane and a transfer function of
H
li
(
=
1
−
α
N
Hz. Frequency variable
k
,
in the range 0
need not be
an integer. The behavior of the net-
work is shown by the solid curves in
Table 2. However, the frequency
magnitude response of the Goertzel
algorithm, for
N
≤
k
<
N
,
z
−
1
]. Small
z
)
=
α/
[1
−
(
1
−
α)
values for
yield narrow passbands
at the expense of increased filter
response time. Table 1 shows the fil-
ter’s behavior for
α
1 as solid
curves. For comparison, the frequen-
cy domain performance for
α
=
0
.
+
delay
samples, where
delay
is typically in the range of
−
2
=
8 and
k
=
1, is
shown as the dashed curve.
After
N
α
=
0
.
5
For example, when
delay
is 0.3, the network delay (at
low frequencies) is 2.3 samples. The
real-valued coefficients are
0
.
5to0
.
5
.
+
1 input samples are
is indicated by the dashed curves.
applied,
y
is a single-bin DFT
result. The DFT computational
workload is
N
(
n
)
First-Order Delay Network
A subclass of a first-order IIR filter,
the coefficients in Table 1 yield an
all-pass
network having a relatively
constant group delay at low fre-
quencies. The network’s delay is
+
2 real multiplies and
2
N
1 real adds. The network is
typically stable because
N
is kept fair-
ly low (in the hundreds) in practice
before the network is reinitialized
[2], [3].
+
−
2
delay
delay
+
R
1
=
3
and
MAY 2004
IEEE SIGNAL PROCESSING MAGAZINE
91
dsp
tips & tricks
continued
Table 1. General functions.
Functions and
Network Behavior
Coefficients
Impulse Response
z
-plane
Magnitude (dB)
Phase (rad.)
Moving Averager
a
0
=
0
5
0.2
1
1
,
a
1
=
1
,
a
2
=
0
,
1/N
0.1
b
0
=
1
/
N
,
b
1
=
0
,
b
2
=
0
,
0
–10
0
0
c
1
=
1
,
N
=
8
–1
–0.1
–20
–5
0
5
10
–1
0
1
–0.5
0
0.5
–0.5
0
0.5
Time
Real Part
Frequency
Frequency
0
2
1
Differencer
a
0
=
1
1
,
a
1
=
0
,
a
2
=
0
,
0
0
–10
0
b
0
=
1
,
b
1
=−
1
,
b
2
=
0
,
c
1
=
0
–1
–1
–20
–2
0
5
10
–1
0
1
–0.5
0
0.5
–0.5
0
0.5
Time
Real Part
Frequency
Frequency
Integrator
a
0
=
0
2
1
1
1
,
a
1
=
1
,
a
2
=
0
,
0.5
0
–10
0
b
0
=
1
,
b
1
=
0
,
b
2
=
0
,
c
i
=
0
0
–1
–2
–20
0
5
10
–1
0
1
–0.5
0
0.5
0.5
0
0.5
Time
Real Part
Frequency
Frequency
0.2
0
2
Leaky Integrator
a
0
=
1
α
= 0.5
α
= 0.1
= 0.1
α
1
,
a
1
=
1
−
α,
a
2
=
0
,
0.1
α
= 0.5
0
–10
0
b
0
=
α,
b
1
=
0
,
b
2
=
0
,
0
α
= 0.1
c
1
=
0
,α
=
0
.
1
–1
–20
–2
0
10
20
–1
0
1
–0.5
0
0.5
–0.5
0
0.5
Real Part
Frequency
Frequency
Time
Group Delay
1
1
First-Order Delay Network
a
0
=
1.2
0
R = 0.91
Delay
∆
= 0.2
∆
= 0.2
Delay
1
,
a
1
=−
R
,
a
2
=
0
,
0.2/f
s
0.5
0
1.1
b
0
=
R
,
b
1
=
1
,
b
2
=
0
,
Zero at
z
= 11
0
c
1
=
0
–1
1
0
5
10
–1
0
1
–0.5
0
0.5
–0.5
0
0.5
Time
Real Part
Frequency
Frequency
Second-Order Delay Network
a
0
=
Group Delay
1
1
2.5
0
R
1
= –0.182
R
2
= 0.028
Delay
1
,
a
1
=−
R
1
,
a
2
=−
R
2
,
2.3
0.3/f
s
0.5
∆
= 0.3
0
b
0
=
R
2
,
b
1
=
R
1
,
b
2
=
1
,
2
Zeros
Not
Shown
c
1
=
0
∆
= 0.3
0
Delay
–1
1.5
–0.5
0
0.5
0
5
10
–1
0
1
–0.5
0
0.5
Time
Real Part
Frequency
Frequency
Sliding DFT Network
This structure computes a single-
bin
N
-point DFT centered at an
angle of
size and integer
k
is 0
N
.
The real damping factor
r
is kept as
close to, but less than, unity as pos-
sible to maintain network stability.
After
N
input samples have been
applied, this network will compute
≤
k
<
a new follow-on DFT result based
on each new
x
sample (thus the
term
sliding
) at a computational
workload of only four real multi-
plies and four real adds per input
sample [2], [3]. Setting coefficient
(
n
)
N
rad on the unit
circle, corresponding to a cyclic fre-
quency of
kf
s
/
θ
=
2
π
k
/
N
Hz.
N
is the DFT
92
IEEE SIGNAL PROCESSING MAGAZINE
MAY 2004
Table 2. Analysis and synthesis functions.
Functions and
Network Behavior
Coefficients
Impulse Response
z
-plane
Magnitude (dB)
Phase (rad.)
0
10
Goertzel Network
a
0
=
1
1
N = 8
k = 1
1
,
a
1
=
2cos
(θ),
5
0
0
–10
e
−
j
θ
,
a
2
=−
1
,
b
0
=
1
,
b
1
=−
0
Real
Part
b
2
=
0
,
c
1
=
0
,θ
=
2
π
k
/
N
–1
–1
–20
–5
0
5
10
–1
0
1
–0.5
0
0.5
–0.5
0
0.5
Time
Real Part
Frequency
Frequency
0
5
1
Sliding DFT Network
a
0
=
1
N = 22
re
j
θ
,
k = 2
a
1
=
1
,
a
2
=
0
,
θ
= 4
π
/22
0
0
–10
0
b
0
=
1
,
b
1
=
0
,
b
2
=
0
,
r = 0.999
Real
Part
r
N
c
1
=
,θ
=
2
π
k
/
N
–1
–1
–20
–5
0
10
20
–1
0
1
–0.5
0
0.5
–0.5
0
0.5
Time
Real Part
Frequency
Frequency
0
15
2
Real Oscillator
a
0
=
1
1
,
a
1
=
2cos
(θ),
θ
=
π
/4
10
0
0
–10
a
2
=−
1
,
b
0
=
1
,
5
b
1
=
0
,
b
2
=−
1
–1
–2
–20
0
0
20
–1
0
1
–0.5
0
0.5
–0.5
0
0.5
Time
Real Part
Frequency
Frequency
0
10
1
Quadrature Oscillator
a
0
=
1
G = 1
e
j
θ
,
G
(
n
),
a
1
=
θ
=
π
/4
5
Real Part
0
0
–10
a
2
=
0
,
b
0
=
1
,
0
b
1
=
0
,
b
2
=
0
–1
–1
–20
–5
0
10
20
–1
0
1
–0.5
0
0.5
–0.5
0
0.5
Time
Real Part
Frequency
Frequency
0
0.5
Audio Comb
a
0
=
1
1
1
,
a
1
=
0
,
a
2
=
α,
0
0.5
–2
0
b
0
=
1
,
b
1
=
0
,
b
2
=
0
,
c
1
=
0
,α
=
0
.
2
0
–1
–0.5
–4
0
5
10
–1
0
1
–0.5
0
0.5
–0.5
0
0.5
Time
Real Part
Frequency
Frequency
r
N
allows the analysis band
to be centered at an angle of
θ
=
c
1
=−
amplitude is not a function of the
output frequency. The argument for
coefficient
a
1
may need to be scaled so that all
intermediate results are in the prop-
er numerical range [4].
N
rad, corre-
sponding to a cyclic frequency of
(
2
π(
k
+
1
/
2
)/
in Table 2
is
θ
=
f
s
rad, where
f
t
is the
oscillator’s tuned frequency in
hertz. To start the oscillator we set
the
y
2
π
f
t
/
Quadrature Oscillator
Called the
coupled quadrature oscil-
lator,
this structure provides
y
k
+
1
/
2
)
f
s
/
N
Hz.
Real Oscillator
There are many possible digital
oscillator structures, but this net-
work generates a real-valued sinu-
soidal
sample driving the
a
1
multiplier equal to 1 and compute
new output samples as the time
index
n
advances. For fixed-point
implementations, filter coefficients
(
n
−
1
)
(
n
)
=
cos
outputs for
a complex exponential sequence
whose tuned frequency is
f
t
Hz.
The exponent for
a
1
in Table 2 is
(
n
θ)
+
j
sin
(
n
θ)
y
(
n
)
sequence whose
MAY 2004
IEEE SIGNAL PROCESSING MAGAZINE
93
dsp
tips & tricks
continued
θ
=
f
s
rad. To start the oscil-
lator, we set the complex
y
2
π
f
t
/
values are the
N
roots of unity when
we set
H
comb
(
write home about; its transition
regions are so wide that they don’t
actually have distinct passbands and
stopbands. Of course, to ensure sta-
bility,
R
p
must be between zero and
one to keep the pole inside the
z
-
plane’s unit circle, and the closer
R
p
is to unity, the more narrowband is
the filter.
)
sample, driving the
a
1
multiplier,
equal to 1
(
n
−
1
z
)
equal to zero yield-
N
e
j
2
π
k
/
N
N
1 . The
N
zeros on the unit circle result in
frequency response nulls (infinite
attenuation) located at cyclic fre-
quencies of
mf
s
/
ing
z
(
k
)
=
(
)
=
j
0 and begin comput-
ing output samples as the time
index
n
advances. To ensure oscilla-
tor output stability in fixed-point
arithmetic implementations, instan-
taneous gain correction
G
+
N
, where integer
m
is 0
2. The peak gain of
this linear-phase filter is two.
If we set coefficient
c
1
to
≤
m
≤
N
/
must
be computed for each output sam-
ple. The
G
(
n
)
1 in
the comb filter, making its transfer
function
H
alt
,
comb
(
−
First-Order Equalizer
This structure has a frequency mag-
nitude response that is constant
across the entire frequency band (an
allpass filter). It has a pole at
z
sample values will be
very close to unity [5], [6].
(
n
)
z
−
N
, we
obtain an alternate linear-phase comb
filter having zeros rotated counter-
clockwise around the unit circle by an
angle of
z
)
=
1
+
Audio Comb
This structure is a second-order
(the simplest) version of an IIR
comb filter used by audio folks to
synthesize the sound of a plucked-
string instrument. The input to the
filter is random noise samples. The
filter has frequency response peaks
at dc and
R
on the
z
-plane and a zero located at
1
=
N
rad positioning the
zeros at angles of 2
π/
N
rad
on the
z
-plane’s unit circle. The
rotated zeros result in frequency
response nulls located at cyclic fre-
quencies of (
m
π(
k
+
1
/
2
)/
R
*, where * means conjugate.
The value of
R
, which can be real
or complex but whose magnitude
must be less than unity to ensure
stability, controls the nonlinear-
phase response. The equalizer has a
transfer function of
H
1
,
eq
(
/
N
. With
this filter a frequency magnitude peak
is located at 0 Hz (dc).
+
1
/
2
)
f
s
/
2, with dips in the
response located at
±
f
s
/
4. The fil-
ter’s transfer function is
H
ac
(
±
f
s
/
z
)
=
R
∗
+
z
−
1
Rz
−
1
z
)
=
(
−
.
These networks can be used as
phase equalizers
by cascading them
after a filter or network whose non-
linear phase response requires crude
linearization. The goal is to make
the cascaded filters’ combined phase
as linear as possible. Table 3 shows
the filter’s behavior for
R
)/(
1
−
)
Bandpass Filter at
fff
s
4
This network is a bandpass filter cen-
tered at
f
s
/
/
z
−
2
1
/(
1
−
α
)
, resu
lt
ing in two poles
=±
√
α
on the
z
-plane.
To maintain stability the real-valued
α
located at
z
x
-like
frequency response and linear-phase
over the passband. It has poles at
z
4 having a sin
(
x
)/
must be less than unity, and the
closer
is to unity, the more nar-
row the frequency response peaks.
For a more realistic-sounding
synthesis, we can set
a
1
=
α
α
zero cancella-
tion the comb filter’s delay
=±
j
, so for pole
/
must
be an integer multiple of four. This
guaranteed-stable, multiplierless,
bandpass filter’s transfer function is
H
bp
(
(
N
)
and the
top delay element of the biquad in
Figure 1 may have its delay
increased to, say, eight instead of
one, yielding more frequency
response peaks between 0 and
f
s
/
7 as
solid curves. For comparison, the
phase response for
R
=
0
.
3 is
indicated by the dashed curve.
These first-order allpass filters can
also be used for interpolation and
audio reverberation for low-fre-
quency signals.
=−
0
.
z
−
N
z
−
2
z
)
=
(
1
−
)/(
1
+
)
.
First-Order IIR Filter
This is the direct form II version of
a simple first-order IIR filter having
a single pole located at a radius of
R
p
from the
z
-plane’s origin at an
angle of
2
Hz. In this music application, the
filter’s input is Gaussian white noise
samples. Other plucked-string
instrument synthesis networks have
been used with success [7], [8].
Second-Order IIR Filter
This is the direct form II version of a
second-order IIR filter, the
workhorse of IIR filter implementa-
tions. Conjugate pole and zero pairs
may be positioned anywhere on the
z
-plane to control the filter’s fre-
quency response. [There’s a terrific
piece of MATLAB code (PEZ, cre-
ated by the talented Dr. Craig
Ulmer) that allows us to see the fre-
quency-domain effect of moving
multiple poles and zeros, manually
θ
p
rad and a single zero at
a radius of
R
z
at an angle of
π
+
θ
z
.
Comb Filter
Referring to the first entry in Table
3, this standard comb filter is a key
component on many filtering appli-
cations, as we shall see. Its transfer
function,
For
real-valued
coefficients
(
0), the filter can only
exhibit either a low-pass or a
high-pass frequency response; no
bandpass or bandstop filters are pos-
sible. The filter’s transfer function is
H
1
,
iir
(
θ
p
=
θ
z
=
z
−
N
,
results in
N
zeros equally spaced
around the
z
-plane’s unit circle
located at
z
H
comb
(
z
)
=
1
−
R
p
e
j
θ
z
z
−
1
z
)
=
(
1
+
)/
R
p
e
j
θ
p
z
−
1
(
.
The shape of the filter’s frequency
magnitude responses are nothing to
1
−
)
e
j
2
π
k
/
N
, where
(
k
)
=
integer
k
is 0
≤
k
<
N
. Those
z
(
k
)
94
IEEE SIGNAL PROCESSING MAGAZINE
MAY 2004
[ Pobierz całość w formacie PDF ]