MUS421/EE367B Lecture 2 Review of the Discrete Fourier Transform (DFT) Julius O. Smith III (
[email protected]) Center for Computer Research in Music and Acoustics (CCRMA) Department of Music, Stanford University Stanford, California 94305 April 10, 2007
Outline • Domains of Definition • Discrete Fourier Transform • Properties of the Fourier Transform For more details, see • Mathematics of the DFT (Music 320 text): http://ccrma.stanford.edu/~jos/mdft/ • Appendices A-D of Spectral Audio Signal Processing (our text): http://ccrma.stanford.edu/~jos/sasp/ 1
Domains of Definition The Fourier Transform can be defined for signals that are • Discrete or Continuous Time • Finite or Infinite Duration This results in four cases: Time Duration Finite Infinite Fourier Series (FS) Fourier Transform (FT) cont. Z P Z +∞ X(k) = x(t)e−jωk tdt X(ω) = x(t)e−jωtdt time 0
−∞
k = −∞, . . . , +∞ ω ∈ (−∞, +∞) t Discrete FT (DFT) Discrete Time FT (DTFT) discr. N −1 +∞ X X x(n)e−jωk n X(ω) = x(n)e−jωn time X(k) = n=−∞
n=0
k = 0, 1, . . . , N − 1 discrete freq. k
ω ∈ (−π, +π) continuous freq. ω
2
n
Signal and Transform Notation • n, k ∈ Z (integers) or ZN (integers modulo N ) • x(n) ∈ R (reals) or C (complex numbers) • x ∈ CN means x is a length N complex sequence • x = x(·) • X = DFT(x) ∈ CN , or x↔X where “↔” is read as “corresponds to”. • X(k) = DFTk (x) = DFTN,k (x) ∈ C • x(n) = IDFTn(X) = IDFTN,n(X) • For x ∈ C∞, X = DTFT(x) = DFT∞(x) ∈ C∞ 2π • x = conjugate of x • ∠x = phase of x The notation XY or X · Y denotes the vector containing (XY )k = X(k)Y (k), k = 0, . . . , N − 1. This is denoted by ‘X .* Y’ in Matlab, where X and Y may a pair of column vectors, or a pair of row vectors.
3
The Discrete Fourier Transform The “kth bin” of the Discrete Fourier Transform (DFT) is defined as ∆
∆
∆
X(k) = DFTk (x) = hx, sk i =
N −1 X
x(tn)e−jωk tn
n=0
∆
sk (n) = ejωk tn ;
k = 0, 1, . . . , N − 1
k 2πk ωk = 2π fs = ; N NT ∆
∆
tn = nT
We may interpret the DFT as the coefficients of projection of the signal vector x onto the N sinusoidal basis signals sk , k = 0, 1, . . . , N − 1: X(k) = hx, sk i
4
Inverse DFT The inverse DFT is given by x(tn) =
N −1 X k=0
N −1 1 X hx, sk i jωk tn s (t ) = X(ω )e k n k N k s k k2 k=0
It can be interpreted as the superposition of the projections, i.e., the sum of the sinusoidal basis signals weighted by their respective coefficients of projection: X hx, sk i x= 2 sk k sk k k
5
The DFT, Cont’d There are several ways to think about the DFT: 1. Sampled uniform filter bank output 2. Projection onto the set of “basis” sinusoids (frequencies at N roots of unity) 3. Coordinate transformation (“natural” RN basis to “sinusoidal” basis) 4. Matrix multiplication X = W∗x, where W∗[k, n] = e−jωk tn This course will emphasize the first two.
6
Properties of the DFT We are going to be performing manipulations on signals and their Fourier Transform throughout this class. It is important to understand how changes we make in one domain affect the other domain. The Fourier theorems are helpful for this purpose. Derivations of the Fourier theorems for the DTFT case may be found in Appendix D1 of the text, and in Mathematics of the DFT2 (Music 320 text) for the DFT case.
1 2
http://ccrma.stanford.edu/~jos/sasp/Fourier_Theorems_DTFT.html http://ccrma.stanford.edu/~jos/mdft/Fourier_Theorems.html
7
Linearity αx1 + βx2 ↔ αX1 + βX2 or DFT(αx1 + βx2) = α · DFT(x1) + β · DFT(x2) α, β ∈ C x1, x2, X1, X2 ∈ CN The Fourier Transform “commutes with mixing.”
8
Time Reversal Definition: ∆
∆
Flipn(x) = x(−n) = x(N − n) ∆
Note: x(n) = x(n mod N ) for signals in CN (DFT case). When computing a sampled DTFT using the DFT, we interpret time indices n = 1, 2, . . . , N/2 − 1 as positive time indices, and n = N − 1, N − 2, . . . , N/2 as the negative time indices n = −1, −2, . . . , −N/2. Under this interpretation, the Flip operator simply reverses a signal in time. Fourier theorems: Flip(x) ↔ X for x ∈ CN . In the typical special case of real signals (x ∈ RN ), Flip(x) ↔ X Thus, time-reversing a signal conjugates its spectrum.
9
Symmetry The Fourier transform of real signals exhibits special symmetry which is important to us in many cases. Basically, Real ↔ Hermitian The Fourier transform of a real signal x is therefore • Conjugate Symmetric (Hermitian symmetric) X(−k) = X(k) • Real part Symmetric (even) re {X(−k)} = re {X(k)} • Imaginary part Antisymmetric (skew-symmetric, odd) im {X(−k)} = −im {X(k)} • Magnitude Symmetric (even) |X(−k)| = |X(k)| • Phase Antisymmetric (odd) ∠X(−k) = −∠X(k) 10
Shift Theorem ∆
The Shift operator is defined as Shiftl,n(y) = y(n − l). Since indexing is defined modulo N , Shiftl (y) is a circular right-shift by l samples. Shiftl (y) ↔ e−j(·)l Y or, more loosely, y(n − l) ↔ e−jωl Y (ω) i.e., −jωk l
DFTk [Shiftl (y)] = e
Y (ωk )
e−jωk l = Linear Phase Term, slope = −l • ∠Y (ωk ) += − ωk l • Multiplying a spectrum Y by a linear phase term e−jωk l with phase slope −l corresponds to a circular right-shift in the time domain by l samples: • negative slope ⇒ time delay • positive slope ⇒ time advance
11
Convolution The cyclic convolution of x and y is defined as ∆
(x ∗ y)(n) =
N −1 X
x(m)y(n − m),
x, y ∈ CN
m=0
Cyclic convolution is also called circular convolution, ∆ since y(n − m) = y(n − m (mod N )). Convolution is cyclic in the time domain for the DFT and FS cases, and acyclic for the DTFT and FT cases. The Convolution Theorem is then (x ∗ y) ↔ X · Y
12
Correlation The cross-correlation of x and y in CN is defined as: ∆
(x ⋆ y)(n) =
N −1 X
x(m)y(n + m),
x, y ∈ CN
m=0
Using this definition we have the correlation theorem: (x ⋆ y) ↔ X(ωk )Y (ωk ) The correlation theorem is often used in the context of spectral analysis of filtered noise signals. Autocorrelation The autocorrelation of a signal x ∈ CN is simply the cross-correlation of x with itself: ∆
(x ⋆ x)(n) =
N −1 X
x(m)x(m + n),
m=0
From the correlation theorem, we have (x ⋆ x) ↔ |X(ωk )|2
13
x ∈ CN
Power Theorem The inner product of two signals is defined as: X ∆ hx, yi = xn y n n
Using this notation, we have the following: 1 hx, yi = hX, Y i N When we consider the inner product of a signal with itself, we have a special case known as Parseval’s Theorem: kXk2 1 kxk = hx, xi = hX, Xi = N N (Also called the Rayleigh’s Energy Theorem.) 2
14
Stretch We define the Stretch operator such that: StretchL : CN → CN L Which means that it transforms a length N complex signal, into a length N L signal. Specifically, we do this by inserting L − 1 zeros in between each pair of samples of the signal. x
y y = Stretch2(x) →
...
...
15
Repeat or Scale Similarly, the RepeatL operator, defined on the unit circle, frequency-scales its input spectrum by the factor L: ω ← Lω The original spectrum is repeated L times as ω traverses the unit circle. This is illustrated in the following diagram for L = 2: X
Y
Y
= REPEAT3(X) →
ω
ω
Using these definitions, we have the Stretch Theorem: StretchL(x) ↔ RepeatL(X) Application: Upsampling by any integer factor L: Passing the stretched signal through an ideal lowpass filter cutting off at ω ≥ π/L yields ideal bandlimited interpolation of the original signal by the factor L. 16
Zero-Padding ↔ Interpolation Zero padding in the time domain corresponds to ideal interpolation in the frequency domain. Proof: http://ccrma.stanford.edu/~jos/mdft/Zero Padding Theorem Spectral.html
Downsampling ↔ Aliasing The downsampling operation DownsampleM selects every M th sample of a signal: ∆
DownsampleM,n(x) = x(M n) N
In the DFT case, DownsampleM maps CN to C M , while for the DTFT, DownsampleM maps C∞ to C∞. The Aliasing Theorem states that downsampling in time corresponds to aliasing in the frequency domain: DownsampleM (x) ↔
1 AliasM (X) M
where the Alias operator is defined for X ∈ CN
17
(DFT case) as ∆
AliasM,l (X) =
M −1 X
X l+k
k=0
N M
,
l = 0, 1, . . . ,
N −1 M
For X ∈ C∞ (DTFT case), the Alias operator is given by ∆
AliasM,ω (X) = ∆
=
M −1 X
k=0 M −1 X k=0
ω +k 2π ) j( M M
X e
k z X WM
1 M
,
−π ≤ ω < π
∆
where WM = ej2π/M is a common notation for the primitive M th root of unity, and z = ejω as usual. This normalization corresponds to T = 1 after downsampling. Thus, T = 1/M prior to downsampling. The summation terms above for k 6= 0 are called aliasing components. The aliasing theorem points out that in order to downsample by factor M without aliasing, we must first lowpass-filter the spectrum to [−πfs/M, πfs/M ]. This filtering essentially zeroes out the spectral regions which alias upon sampling. 18