Digital Signal Processing L03 Fourier

Download as pdf or txt
Download as pdf or txt
You are on page 1of 59

ELEN E4810: Digital Signal Processing

Topic 3: Fourier domain


1. The Fourier domain
2. Discrete-Time Fourier Transform (DTFT)
3. Discrete Fourier Transform (DFT)
4. Convolution with the DFT

Dan Ellis 2003-09-16 1


1. The Fourier Transform
 Basic observation (continuous time):
A periodic signal can be decomposed
into sinusoids at integer multiples of the
fundamental frequency
 i.e. if x (t) = x (t +T )
we can approach x with Harmonics
M
2k of the
x (t) ak cos t + k fundamental
k=0
T
Dan Ellis 2003-09-16 2
M
2k
Fourier Series ak cos T t + k
 For a square wave, k=0

k1
1
k = 0; ak = (1) 2 k = 1,3,5,...
k
0 o.w.
2 1 2 1 2
i.e. x(t) = cos t cos 3t + cos 5t ...
1
T 3 T 5 T
0.5

0.5

1
1.5 1 0.5 0 0.5 1 1.5

Dan Ellis 2003-09-16 3


Fourier domain
 x is equivalently described
by its Fourier Series ak 1.0

parameters:
1 2 3 4 5 6 7 k
k1
1
ak = (1) 2 k = 1,3,5... k

k Negative ak is
equivalent to phase of 1 2 3 4 5 6 7 k
M 2k
j t
 Complex form: x (t) ck e T

k=M

Dan Ellis 2003-09-16 4


Fourier analysis
 How to find {|ck|}, {arg[ck]}?
Inner product with complex sinusoids:
T /2 2 k
1 j

t
ck = x(t)e dt T
2 T / 2
1
= ( x (t ) cos j 2Tkt dt j x (t ) sin j 2Tkt dt )
2

Dan Ellis 2003-09-16 5


Fourier analysis
Works if k1, k2 are positive integers,



1
1 k1 = k 2


cos(k1t ) cos(k 2t ) dt =
0 k1 k 2
1
=
2
cos(k 1 + k 2 )t + cos(k 1 k 2 )t dt

1 sin(k1 + k 2 )t sin(k1 k 2 )t
= +
2 k1 + k 2 k1 k 2
= sinc (k1 + k 2 ) + sinc (k1 k 2 )
Dan Ellis 2003-09-16 6
sinc
sin x
 sincx =
x

sin(x)/x
1.0

4 3 2 2 3 4 x
sin(x)
y=x

 = 1 when x = 0
= 0 when x = r, r 0, r = 1, 2, 3,...
Dan Ellis 2003-09-16 7
Fourier Analysis
2 k
1 T /2 j

t
 Thus, ck = x(t)e T
dt
2 T / 2 2 k
j t
because real & imag sinusoids in e T

pick out corresponding sinusoidal


components linearly combined in
M
2k
x (t) = ak cos t + k
k=0
T

Dan Ellis 2003-09-16 8


Fourier Transform
 Fourier series for periodic signals
extends naturally to Fourier Transform
for any (CT) signal (not just periodic):
1 Fourier
X() =
2
x(t)e jt
dt Transform (FT)


x(t) = X()e jt
d Inverse Fourier
Transform (IFT)

 Discrete index k continuous freq.


Dan Ellis 2003-09-16 9
Fourier Transform
 Mapping between two continuous
functions: |X( )| 0
level
/ dB -20

-40

0.02
x(t)
-60

-80
0.01 0 2000 4000 6000 8000
freq / Hz

0
arg{X( )}

-0.01
0 0.002 0.004 0.006 0.008 /2
time / sec
0

/2


2 ambiguity 0 2000 4000 6000 8000
freq / Hz

Dan Ellis 2003-09-16 10


Fourier Transform of a sine
j 0 t
 Assume x(t) = e
Now, since x(t) = X()e jt
d

...we know X() = ( 0 )
...where (x) is the Dirac delta function
(x-x0)
(continuous time) i.e.
( x x0 ) f ( x)dt = f ( x0 ) x0
f(x)
x
j 0 t
 x(t) = Ae X() = A ( 0 )
Dan Ellis 2003-09-16 11
Fourier Transforms
Time Frequency
Fourier Series Continuous Discrete
(FS) periodic ~x(t) infinite ck
Fourier Continuous Continuous
Transform (FT) infinite x(t) infinite X()
Discrete-Time Discrete Continuous
FT (DTFT) infinite x[n] periodic X(ej)
Discrete FT Discrete Discrete
(DFT) finite/pdc ~x[n] finite/pdc X[k]

Dan Ellis 2003-09-16 12


2. Discrete Time FT (DTFT)
 FT defined for discrete sequences:

j
X(e ) = x[n]e jn
DTFT
n=
 Summation (not integral)
 Discrete (normalized)
frequency variable
 Argument is ej, not plain

Dan Ellis 2003-09-16 13


DTFT example
 e.g. x[n] = n[n], || < 1 -1 1 2 3 4 5 6 7 n

X(e ) = j
[n]e n jn
n=

= n=0 (e )
j n
S = c cS = c
n n

n=0 n=1
1
= j S cS = c = 1 0

1 e
1
3
|X(ej )|

arg{X(ej )}
S=
2
2 3 4 5
1 c
1


0 2 3 4 5

Dan Ellis 2003-09-16 14


Periodicity of X(ej)
 X(ej) has periodicity 2 in :
X(e j ( +2 )
) = x[n]e j ( +2 )n

= x[n]e jn
e j 2 n
= X(e ) j

|X(ej )| arg{X(ej )}
3

2
2 3 4 5
1


0 2 3 4 5

 Phase ambiguity of ej makes it implicit

Dan Ellis 2003-09-16 15


Inverse DTFT (IDTFT)
 Same basic form as other IFTs:
1
x[n] =
2
X(e j
)e jn
d IDTFT

 Note: continuous, periodic X(ej)


discrete, infinite x[n] ...
 IDTFT is actually forward Fourier Series
(except for sign of )

Dan Ellis 2003-09-16 16


IDTFT
 Verify by substituting in DTFT:
1
x[n] =
2
X(e j
)e jn
d
1
(l x[l]e )e
jl jn
= d
2
1 j (nl ) = 0 unless
= l x[l] e d n=l
2
= l x[l]sinc (n l)= x[n]
Dan Ellis 2003-09-16 17
sinc again
1 1 e
j (nl )

2
e j (nl )
d =
2 j (n l )

1 e j (nl )
e j (nl )
=
2 j (n l )
1 2 j sin (n l)
= = sinc (n l)
2 j (n l )
Same as cos imag jsin part cancels



Dan Ellis 2003-09-16 18


DTFTs of simple sequences

 x[n] = [n] X(e j ) = x[n]e jn

n=
j 0
=e =1 (for all )

 i.e. x[n] X(ej)


[n] 1
x[n] X(ej)

-3 -2 -1 1 2 3 n -

Dan Ellis 2003-09-16 19


DTFTs of simple sequences
1 IDTFT
 x[n] = e : x[n] =
2

j 0 n

X(e j
)e jn
d
j
X(e ) = 2 ( 0 ) over - < <
but X(ej) must be periodic in
e j 0 n
k 2 ( 0 2k)
 If 0 = 0 then x[n] = 1 n
so 1 k 2 ( 2k)
Dan Ellis 2003-09-16 20
DTFTs of simple sequences
 From before:
n 1
[n] j ( || < 1)
1 e

 [n] tricky - not finite


1
[n] j
+ k ( + 2k)
1 e
DTFT of 1/2
Dan Ellis 2003-09-16 21
DTFT properties
 Linear:
j j
g[n]+ h[n] G(e ) + H (e )
 Time shift:
jn 0 j
g[n n0 ] e G(e )
 Frequency shift: delay
e j 0 n
g[n] G (e j ( 0 )
) in
frequency

Dan Ellis 2003-09-16 22


DTFT example
 x[n] = [n] + n[n-1] ?
= [n] + (n-1[n-1])
j j 1 1
X(e ) = 1+ e
1 e
j

e j 1 e j + e j
= 1+ j
=
1 e 1 e j
1
= x[n] = n[n]
1 e j
Dan Ellis 2003-09-16 23
DTFT symmetry
If x[n] X(ej) then...
x[-n] X(e-j) from summation
(e -j)* = ej
x*[n] X*(e-j)
1
Re{x[n]} XCS(e )
j =
2
[ X (e j
) + X *
(e j
)]
conjugate symmetry cancels Im parts on IDTFT
1
jIm{x[n]} XCA = [ X (e ) X (e )]
(ej)j * j

2
xcs[n] Re{X(ej)}
xca[n] jIm{X(ej)}
Dan Ellis 2003-09-16 24
DTFT of real x[n]
 When x[n] is pure real, X(ej) = X*(e-j)
xcs[n] xev[n] = xev[-n] XR(ej) = XR(e-j)
xca[n] xod[n] = -xod[-n] XI(ej) = -XI(e-j)

Imag
xim[n]

x[n] real, even Real

X(ej) even, real xre[n]

Dan Ellis 2003-09-16 25


DTFT and convolution
 Convolution: x[n] = g[n] h[n]

X(e ) = n= ( g[n] h[n])e
j jn

= n (
k g[k]h[n k] e jn
)
= k g[k]e jk
n h[n k]e j (nk )

j j
= G(e ) H (e )
Convolution
g[n] h[n] G(e j )H (e j ) becomes
multiplication
Dan Ellis 2003-09-16 26
DTFT modulation
Modulation: x[n] = g[n] h[n]

Could solve if g[n] was just sinusoids...
1
X(e ) = n
j
G(e j
)e jn
d h[n]e jn

2
1
=
2

G(e ) n h[n]e
j
[ j ( )n
] d
1
g[n] h[n]
2
G(e j
)H (e j ( )
)d
Dual of convolution in time
Dan Ellis 2003-09-16 27
Parsevals relation
 Energy in time and frequency domains
are equal:
1
g[n]h *
[n] =
2

G(e j
)H *
(e j
)d
n

 If g = h, then gg* = |g|2 = energy...

Dan Ellis 2003-09-16 28


Energy density spectrum
2
 Energy of sequence g = g[n]
n

1 2
 By Parseval g =
2
G(e j
) d

 Define Energy Density Spectrum (EDS)


j j 2
Sgg (e ) = G(e )

Dan Ellis 2003-09-16 29


EDS and autocorrelation
 Autocorrelation of g[n]:

rgg [l] = g[n]g[n l] = g[n] g[n]
n=

DTFT {rgg [l]} = G(e j )G(e j )


 If g[n] is real, G(e-j) = G*(ej), so
2 no phase
DTFT {rgg [l]} = G(e ) = Sgg (e )
j j
info.

 Mag-sq of spectrum is DTFT of autoco


Dan Ellis 2003-09-16 30
Convolution with DTFT
j j
 Since g[n] h[n] G(e )H (e )
we can calculate a convolution by:
 finding DTFTs of g, h G, H
 multiply them: GH
 IDTFT of product is result, g[n] h[n]
G(e j )
g[n] DTFT
Y (e j )
IDTFT y[n]
h[n] DTFT
H (e j )

Dan Ellis 2003-09-16 31


DTFT convolution example
1
 x[n] = n[n] X(e ) = j

1 e j
 h[n] = [n] - [n-1]
H (e j ) = 1 (e j 1 ) 1

 y[n] = x[n] h[n]


Y (e j ) = H (e j )X(e j )
1
j ( )
j
= 1 e =1
1 e
y[n] = [n] i.e. ...
Dan Ellis 2003-09-16 32
3. Discrete FT (DFT)
Discrete FT Discrete Discrete
(DFT) finite/pdc x[n] finite/pdc X[k]

 A finite or periodic sequence has only


N unique values, x[n] for 0 n < N
 Spectrum is completely defined by N
distinct frequency samples
 Divide 0..2 into N equal steps,
{k} = 2k/N

Dan Ellis 2003-09-16 33


DFT and IDFT
 Uniform sampling of DTFT spectrum:
N 1 2 k
j n
X[k] = X(e j ) = 2 k = x[n]e N

N n=0

N 1
 DFT: X[k] = x[n]WN
kn

n=0

2
j
where WN = e N i.e. 1/Nth of a revolution

Dan Ellis 2003-09-16 34


IDFT
1 N 1
Inverse DFT IDFT x[n] = X[k]WN
nk

N k=0
 Check:
1
N
(
x[n] = k l x[l]WN WN kl nk
)
N 1 N 1 Sum of complete set
1
= x[l] WN k (ln ) of rotated vectors

N l=0 k=0
= 0 if l n; = N if l = n
im
= x[n]
WN2
WN re

Dan Ellis 2003-09-16 35


DFT examples
1 n=0
 Finite impulse x[n] =
0 n = 1..N 1
N 1
X[k] = x[n]WNkn = WN0 = 1 k
n=0

Periodic sinusoid:
2rn 1 rn
x[n] = cos (r I) = (WN + WNrn )
N 2
X[k] = 2 (WNrn + WNrn )WNkn
N 1
1
n=0

N /2 k = r,k = N r
=
0 o.w.
Dan Ellis 2003-09-16 36
DFT: Matrix form
N 1
 X[k] = n=0 x[n] W N
kn
as a matrix multiply:
X[0] 1 1 1 L 1 x[0]
X[1] 1 WN1 WN 2
L WN (N 1)
x[1]
X[2] = 1 W 2 WN 4
L WN2(N 1)
x[2]
N

M M M M O M M
(N 1) 2
X[N 1] 1 WN
(N 1)
WN
2(N 1)
L WN x[N 1]

 i.e. X = DN x
Dan Ellis 2003-09-16 37
Matrix IDFT
 If X = DN x
1
then x = DN X
 i.e. inverse DFT is also just a matrix,

1 1 1 L 1
1 WN1 WN2 L WN(N 1)
1
1
D N = 1 WN2 WN4 L WN2(N 1)
N M M M O M

(N 1) 2(N 1) (N 1) 2
1 W N WN L WN
= /NDN
1 *

Dan Ellis 2003-09-16 38


DFT and DTFT

DTFT
j
X(e ) = x[n]e jn continuous freq
infinite x[n], -<n<
n=
N 1
discrete freq k=N/2
DFT X[k] = x[n]W N
kn
finite x[n], 0n<N
n=0

 DFT samples DTFT at discrete freqs:


j X[k]
X[k] = X(e ) = 2 k X(ej)
N
k=1...
Dan Ellis 2003-09-16 39
DFT and MATLAB
 MATLAB is concerned with sequences
not continuous functions like X(ej)
 Instead, we use the DFT to sample X(ej)
on an (arbitrarily-fine) grid:
 X = freqz(x,1,w); samples the DTFT
of sequence x at angular frequencies in w
 X = fft(x); calculates the N-point DFT
of an N-point sequence x
M

Dan Ellis 2003-09-16 40


DTFT from DFT
 N-point DFT completely specifies the
continuous DTFT of the finite sequence
N 1
1 N 1
jn
X(e ) = X[k]WN e
j kn

n=0 N k=0
N 1 N 1
1 j ( 2 k ) n
= X[k] e N periodic
N k=0 n=0
sinc
k = 2Nk
N 1 k
1 sin N 2 j ( N21) k
= X[k] k
e
interpolation N k=0 sin 2
Dan Ellis 2003-09-16 41
Periodic sinc
N 1 jN k
1 e
e j k n
=
1 e j k
n=0
jN k / 2 jN k / 2 jN k / 2
e e e
= j k n / 2 j k / 2 j k / 2
e e e
k
j 2 k sin N 2
( N 1)
=e k
sin 2
 = N when k = 0; = (-)N when k/2 =
= 0 when k/2 = r/N, r = 1, 2, ...
other values in-between...
Dan Ellis 2003-09-16 42
Periodic sinc N
sinNx/
sinx

sin Nx sinNx
0
2 x
sin x sinx

-N

sin N k/2
X[k] = X(ej2k/N) X(ej 0) = X[k]
N sin k/2
1.5

Periodic 1 X[3]
sin N 3/2
sinc N sin 3/2
interpolation 0.5

X[k]X(ej) 0
freq
-0.5
-1 0 k=1 k=3 k=4
= 2/N 0 = 6/N = 8/N

Dan Ellis 2003-09-16 43


DFT from overlength DTFT
 If x[n] has more than N points, can still
j
form X[k] = X(e ) = 2 k
N

 IDFT of X[k] will give N point x[n]


 How does x[n] relate to x[n] ?

Dan Ellis 2003-09-16 44


DFT from overlength DTFT
DTFT sample IDFT
x[n]
X(e )
j X[k] x[n]
-A n < B 0n<N
N 1

nk
x[n] = N x[l]WN WN
1 kl
=1 for n-l = rN, rI
k=0 l=
= 0 otherwise

1 N 1

= x[l] N WN k (nl)

l= k=0
all values shifted by
x[n] = x[n rN] exact multiples of N pts
r= to lie in 0 n < N
Dan Ellis 2003-09-16 45
DFT from DTFT example
 If x[n] = { 8, 5, 4, 3, 2, 2, 1, 1} (8 point)
 We form X[k] for k = 0, 1, 2, 3
by sampling X(ej) at = 0, /2, , 3/2

 IDFT of X[k] gives 4 pt x[n] = x[n rN]
r=
 Overlap only for r = -1: (N = 4)
8 5 4 3
x[n] = + + + + = {10 7 5 4}
2 2 1 1
Dan Ellis 2003-09-16 46
DFT from DTFT example
 x[n]
-1 1 2 3 4 5 6 7 8 n

 x[n+N]
(r = -1) -5 -4 -3 -2 -1 1 2 3 4 5 n

 x[n]
1 2 3 n

 x[n] is the time aliased or folded down


version of x[n].
Dan Ellis 2003-09-16 47
Properties: Circular time shift
 DFT properties mirror DTFT, with twists:
 Time shift must stay within N-pt window
g [ n n0 N ] W G[k] N
kn 0

 Modulo-N indexing keeps index between


0 and N-1:
g[n n0 ] n n0
g [ n n0 N ] =
g [ N + n n0 ] n < n0
0 n0 < N

Dan Ellis 2003-09-16 48


Circular time shift
 Points shifted out to the right dont
disappear they come in from the left
g[n] g[<n-2>5]
delay by 2

1 2 3 4
n 1 2 3 4
n
5-pt sequence

 Like a barrel shifter:


origin pointer

Dan Ellis 2003-09-16 49


Circular time reversal
 Time reversal is tricky in modulo-N
indexing - not reversing the sequence:
x [n]
5-pt sequence
made periodic n
-7 -6 -5 -4 -3 -2 -1 1 2 3 4 5 6 7 8 9 10 11

x [ n N ]
Time-reversed
periodic sequence n
-7 -6 -5 -4 -3 -2 -1 1 2 3 4 5 6 7 8 9 10 11

 Zero point stays fixed; remainder flips


Dan Ellis 2003-09-16 50
4. Convolution with the DFT
 IDTFT of product of DTFTs of two N-pt
sequences is their 2N-1 pt convolution
 IDFT of the product of two N-pt DFTs
can only give N points!
 Equivalent of 2N-1 pt result time aliased:

 i.e. yc [n ] = y l [n + rN]
r=
 must be, because G[k]H[k] are exact
samples of G(ej)H(ej)
 This is known as circular convolution
Dan Ellis 2003-09-16 51
Circular convolution
 Can also do entire convolution with
modulo-N indexing
 Hence, Circular Convolution:
N 1

g[m ]h[ n m N ] G[k]H[k]


m=0

 Written as g[n]
N h[n]

Dan Ellis 2003-09-16 52


Circular convolution example
 4 pt sequences: g[n]={1 2 0 1} h[n]={2 2 1 0}
N 1

g[m ]h[ n m N ] 1 2 3 n 1 2 3 n
m=0
g[n] 4 h[n]={4 7 5 4}
h[<n - 0>4]
1 2 3
n 1
h[<n - 1>4]
1 2 3
n 2
1 2 3 n
h[<n - 2>4]
1 2 3
n 0
h[<n - 3>4] check: g[n] * h[n]
1 2 3
n 1 ={2 6 5 4 2 1 0}
Dan Ellis 2003-09-16 53
Duality
 DFT and IDFT are very similar
 both map an N-pt vector to an N-pt vector

 Duality:
Circular
if g[n] G [k ] time reversal
then G [n] N g[ k N ]
 i.e. if you treat DFT sequence as a time
sequence, result is almost symmetric

Dan Ellis 2003-09-16 54


DFT properties summary
 Circular convolution
N 1
m=0 g[m ]h[ n m N ] G[k]H[k]
 Modulation N 1
g[n] h[n] 1
N m=0 G[m]H[ k m N ]
Duality
G [n] N g[ k ]

N

N 1 N 1
Parseval
n=0 x[n] k=0 X [k ]
 2 2
= 1
N

Dan Ellis 2003-09-16 55


Linear convolution w/ the DFT
 DFT fast circular convolution
 .. but we need linear convolution
 Circular conv. is time-aliased linear
conv.; can aliasing be avoided?
 e.g. convolving L-pt g[n] with M-pt h[n]:
y[n] = g[n] * h[n] has L+M-1 nonzero pts
 Set DFT size N L+M-1 no aliasing

Dan Ellis 2003-09-16 56


Linear convolution w/ the DFT
 Procedure (N = L + M - 1): g[n]
 pad L-pt g[n] with (at least)
M-1 zeros L n
N-pt DFT G[k], k = 0..N-1 h[n]
 pad M-pt h[n] with (at least)
M n
L-1 zeros
N-pt DFT H[k], k = 0..N-1 yc[n]

 Y[k] = G[k]H[k], k = 0..N-1


Nn
 IDFT{Y[k]} = yL [n + rN ] = yL [n ]
r=

Dan Ellis 2003-09-16 57


Overlap-Add convolution
 Very long g[n] break up into
segments, convolve piecewise, overlap
bound size of DFT, processing delay
g[n] i N n < (i +1) N
 Make gi [n] =
0 o.w.
g[n] = igi [n]
h[n] * g[n] = ih[n] * gi [n]
 Called Overlap-Add (OLA) convolution...
Dan Ellis 2003-09-16 58
Overlap-Add convolution
g[n] h[n]

n n
g0[n] h[n] g0[n]

n n
g1[n] h[n] g1[n]
n n
g2[n] h[n] g2[n]
n n
N 2N 3N valid OLA sum
h[n] g[n]
n
N 2N 3N
Dan Ellis 2003-09-16 59

You might also like