Because we are interested in actual computations
rather than analytic calculations, we must consider the details
of the discrete Fourier transform. To compute the
length-NN DFT, we assume that the
signal has a duration less than or equal to
NN. Because frequency responses
have an explicit frequency-domain specification in terms of
filter coefficients, we don't have a direct handle on which
signal has a Fourier transform equaling a given frequency
response. Finding this signal is quite easy. First of all, note
that the discrete-time Fourier transform of a unit sample equals
one for all frequencies. Since the input and output of
linear, shift-invariant systems are related to each other by
Yei2πf=Hei2πfXei2πf
Y
2
f
H
2
f
X
2
f
, a unit-sample input, which has
Xei2πf=1
X
2
f
1
, results in the output's Fourier transform
equaling the system's transfer function.
This statement is a very important result. Derive it
yourself.
The DTFT of the unit sample equals a constant (equaling
1). Thus, the Fourier transform of the output equals the
transfer function.
In the time-domain, the output for a unit-sample input is known as the
system's unit-sample response, and is denoted by
hn
h n
.
Combining the frequency-domain and time-domain interpretations of a
linear, shift-invariant system's unit-sample response, we have that
hn
h n
and the transfer function are Fourier transform pairs
in terms of the discrete-time Fourier
transform.
hn↔Hei2πf
↔
h n
H
2
f
(1)
Returning to the issue of how to use the DFT to perform
filtering, we can analytically specify the frequency response,
and derive the corresponding
length-NN DFT by sampling the
frequency response.
∀
k
,k=0…N−1:Hk=Hei2πkN
k
k
0
…
N
1
H
k
H
2
k
N
(2)
Computing the inverse DFT yields a
length-
NN signal
no
matter what the actual duration of the unit-sample response
might be. If the unit-sample response has a duration
less than or equal to
NN (it's a
FIR filter), computing the inverse DFT of the sampled frequency
response indeed yields the unit-sample response. If, however,
the duration exceeds
NN, errors are
encountered. The nature of these errors is easily explained by
appealing to the Sampling Theorem. By sampling in the frequency
domain, we have the potential for aliasing in the time domain
(sampling in one domain, be it time or frequency, can result in
aliasing in the other) unless we sample fast enough. Here, the
duration of the unit-sample response determines the minimal
sampling rate that prevents aliasing. For FIR systems —
they by definition have finite-duration unit sample responses
— the number of required DFT samples equals the
unit-sample response's duration:
N≥q
N
q
.
Derive the minimal DFT length for a
length-qq unit-sample response
using the Sampling Theorem. Because sampling in the
frequency domain causes repetitions of the unit-sample
response in the time domain, sketch the time-domain result
for various choices of the DFT length
NN.
In sampling a discrete-time signal's Fourier transform
LL times equally over
0
2π
0
2
to form the DFT, the corresponding signal equals the periodic
repetition of the original signal.
Sk↔∑
i
=−∞∞sn−iL
↔
S
k
i
s
n
i
L
(3) To avoid aliasing (in the time domain), the
transform length must equal or exceed the signal's duration.
Express the unit-sample response of a FIR filter in terms of
difference equation coefficients. Note that the
corresponding question for IIR filters is far more difficult
to answer: Consider the example.
The difference equation for an FIR filter has the form
yn=∑
m
=0q
b
m
xn−m
y
n
m
0
q
b
m
x
n
m
(4)
The unit-sample response equals
hn=∑
m
=0q
b
m
δn−m
h
n
m
0
q
b
m
δ
n
m
(5)
which corresponds to the representation described in
a
problem of a length-
qq
boxcar filter.
For IIR systems, we cannot use the DFT to find
the system's unit-sample response: aliasing of the unit-sample
response will always occur. Consequently,
we can only implement an IIR filter accurately in the time
domain with the system's difference
equation. Frequency-domain implementations are
restricted to FIR filters.
Another issue arises in frequency-domain filtering that is
related to time-domain aliasing, this time when we consider the
output. Assume we have an input signal having duration
N x
N x
that we pass through a FIR filter having a
length-
q+1
q
1
unit-sample response. What is the duration of the output signal? The
difference equation for this filter is
yn=
b0
xn+…+
b q
xn−q
y
n
b0
x
n
…
b q
x
n
q
(6)
This equation says that the output depends on current and past
input values, with the input value
qq samples previous defining the
extent of the filter's
memory of past input
values. For example, the output at index
N x
N x
depends on
x N x
x
N x
(which equals zero),
x N x −1
x
N x
1
, through
x N x −q
x
N x
q
. Thus, the output returns to zero only after the last input value passes
through the filter's memory. As the input signal's last value occurs at
index
N x −1
N x
1
, the last nonzero output value occurs when
n−q= N x −1
n
q
N x
1
or
n=q+ N x −1
n
q
N x
1
. Thus, the output signal's duration equals
q+ N x
q
N x
.
In words, we express this result as "The output's
duration equals the input's duration plus the filter's duration minus
one.". Demonstrate the accuracy of this statement.
The unit-sample response's duration is
q+1
q
1
and the signal's
N x
N x . Thus the statement is correct.
The main theme of this result is that a
filter's output extends longer than either its input or its
unit-sample response. Thus, to avoid aliasing when we use DFTs,
the dominant factor is not the duration of input or of the
unit-sample response, but of the output. Thus, the number of
values at which we must evaluate the frequency response's DFT
must be at least
q+ N x
q
N x
and we must compute the same length DFT of the input. To
accommodate a shorter signal than DFT length, we simply
zero-pad the input: Ensure that for indices
extending beyond the signal's duration that the signal is
zero. Frequency-domain filtering, diagrammed in Figure 1, is accomplished by storing the
filter's frequency response as the DFT
Hk
H
k
, computing the input's DFT
Xk
X
k
, multiplying them to create the output's DFT
Yk=HkXk
Y
k
H
k
X
k
, and computing the inverse DFT of the result to yield
yn
y
n
.
Before detailing this procedure, let's clarify
why so many new issues arose in trying to develop a
frequency-domain implementation of linear filtering. The
frequency-domain relationship between a filter's input
and output is always true:
Yei2πf=Hei2πfXei2πf
Y
2
f
H
2
f
X
2
f
. The Fourier transforms in this result are
discrete-time Fourier transforms; for example,
Xei2πf=∑nxne−(i2πfn)
X
2
f
n
x
n
2
f
n
. Unfortunately, using this relationship to perform
filtering is restricted to the situation when we have analytic
formulas for the frequency response and the input signal. The
reason why we had to "invent" the discrete Fourier transform
(DFT) has the same origin: The spectrum resulting from the
discrete-time Fourier transform depends on the
continuous frequency variable
ff. That's fine for analytic
calculation, but computationally we would have to make an
uncountably infinite number of computations.
Did you know that two kinds of infinities can be
meaningfully defined? A countably infinite quantity
means that it can be associated with a limiting process
associated with integers. An uncountably infinite
quantity cannot be so associated. The number of rational
numbers is countably infinite (the numerator and denominator
correspond to locating the rational by row and column; the total
number so-located can be counted, voila!); the number of
irrational numbers is uncountably infinite. Guess which is
"bigger?"
The DFT computes the Fourier transform at a
finite set of frequencies — samples the true spectrum
— which can lead to aliasing in the time-domain unless we
sample sufficiently fast. The sampling interval here is
1K
1
K
for a length-
KK DFT: faster
sampling to avoid aliasing thus requires a longer transform
calculation. Since the longest signal among the input,
unit-sample response and output is the output, it is that
signal's duration that determines the transform length. We
simply extend the other two signals with zeros (zero-pad) to
compute their DFTs.
Suppose we want to average daily stock prices
taken over last year to yield a running weekly average
(average over five trading sessions). The filter we want is a
length-5 averager (as shown in the unit-sample response),
and the input's duration is 253 (365 calendar days minus
weekend days and holidays). The output duration will be
253+5−1=257
253
5
1
257
, and this determines the transform length we need to
use. Because we want to use the FFT, we are restricted to
power-of-two transform lengths. We need to choose any FFT
length that exceeds the required DFT length. As it turns out,
256 is a power of two (
28=256
2
8
256
), and this length just undershoots our required
length. To use frequency domain techniques, we must use
length-512 fast Fourier transforms.
Figure 2 shows the input and the
filtered output. The MATLAB programs that compute the filtered
output in the time and frequency domains are
Time Domain
h = [1 1 1 1 1]/5;
y = filter(h,1,[djia zeros(1,4)]);
Frequency Domain
h = [1 1 1 1 1]/5;
DJIA = fft(djia, 512);
H = fft(h, 512);
Y = H.*X;
y = ifft(Y);
The filter program has the
feature that the length of its output equals the length of
its input. To force it to produce a signal having the
proper length, the program zero-pads the input
appropriately.
MATLAB's
fft function
automatically zero-pads its input if the specified transform
length (its second argument) exceeds the signal's
length. The frequency domain result will have a small
imaginary component — largest value is
2.2×10-112.2-11 —
because of the inherent finite precision nature of computer
arithmetic. Because of the unfortunate misfit between signal
lengths and favored FFT lengths, the number of arithmetic
operations in the time-domain implementation is far less
than those required by the frequency domain version: 514
versus 62,271. If the input signal had been one sample
shorter, the frequency-domain computations would have been
more than a factor of two less (28,696), but far more than
in the time-domain implementation.
An interesting signal processing aspect of
this example is demonstrated at the beginning and end of the
output. The ramping up and down that occurs can be traced to
assuming the input is zero before it begins and after it
ends. The filter "sees" these initial and final values as the
difference equation passes over the input. These artifacts can
be handled in two ways: we can just ignore the edge effects or
the data from previous and succeeding years' last and first
week, respectively, can be placed at the ends.
"Electrical Engineering Digital Processing Systems in Braille."