Digital Signal Processing CS3291

1
UNIVERSITY of MANCHESTER Department of Computer
Science CS3291 Digital Signal Processing
05 Section 8 Introduction to the DFT
2
DTFT of xn is
If xn obtained from xa(t), correctly
bandlimited then
for -? lt ? lt??? ? ?T DTFT is
convenient way of calculating Xa(j?) by computer.
3

• Two difficulties
• (i)Infinite range of summation
• (ii) X(e j? ) is continuous function of ?
• Solutions
• (i) Time-domain windowing
• Restrict xn to x0, x1,
  xN-1 ? xn 0, N-1
• (ii) Frequency-domain sampling
• Store values of X(ej?) for -? lt ? lt
  ?
• For real signals we only need 0 ? ? lt
  ? but generalise
  to
  complex signals
• Instead of -? lt ? lt ? take 0 ? ? lt
  2?

4
Why take 0 ? ? lt 2? ?

• X(e j ? ) X( e j (? 2? ) ) for any ?
• So for X(e - j ? / 3 ) look up X(e j 5 ?
  / 3 ) .
• Same information, it is convenient for ? to
  start off at 0.
• In many cases, not interested in ? gt ?

5
Taking M equally spaced samples over 0 ? ? lt 2?
we get Xk 0, M-1 ? X0,
X1,, XM-1 where Xk
X(exp(j?k)) with ?k 2?k/M
6

• For spectral analysis, the larger M, the better
  for drawing accurate graphs etc.
• But, if we need minimum M for storage of
  unambiguous spectrum, take MN.
• DFT xn 0, N-1 ? Xk 0, N-1
• ?
  ?
• (complex) (complex)

?k 2?k/N
7

• DFT transforms a finite sequence to another
  finite sequence.
• DTFT transforms infinite sequence to continuous
  functn of ?

Inverse DFT Xk0, N-1 ? xn0, N-1
Note Similarity with DFT
8
Programming the DFT its inverse
?k 2?k/N

• Similarity exploited by programs able to
  perform DFT or its inverse using same code.
• Programs to implement these equations in a
  direct manner given in MATLAB (using
  complex arith) C (using real arith only).
• These direct programs are very slow FFT is
  much faster.

9
Direct DFT using complex arithmetic in MATLAB
Given N complex time-domain samples in array
x1N E 2pi/N for k0 N-1 X(1k) 0
j0 Wk kE for L 0 N-1 C
cos(LWk) j sin(LWk) X(1k) X(1k)
x(1L) C end end Now have N complex
freq-dom samples in array X1N
10
Direct inverse DFT using complex arithmetic in
MATLAB
Given N complex freq-domain samples in array
x1N E -2pi/N for k0 N-1 X(1k) 0
j0 Wk kE for L 0 N-1 C
cos(LWk) j sin(LWk) X(1k)
X(1k) x(1L) C end X(1k) X(1k)/N
end Now have N complex time-dom samples in
array X1N
11
Direct forward/inverse DFT using complex arith in
MATLAB
Given N complex samples in array x1N if
(Invers 1) E -2pi/N else E 2pi/N
for k0 N-1 X(1k) 0 j0 Wk kE
for L 0 N-1 C cos(LWk) j
sin(LWk) X(1k) X(1k) x(1L) C
end if (Inverse 1) X(1k) X(1k)/N
end Now have N complex samples in array
X1N
12
// Direct fwd/inverse DFT using real arith only
in C void directdft(void) // DFT or Inverse
DFT by direct method. // OrderN, Real
imag parts of input in arrays xr xi //
Output- Real part in array X, Imag part in Y
// Invers is 0 for DFT, 1 for IDFT int k, L
float c,e,s,wk if(Invers1) e
-2.0PI/(float)N else e 2.0PI/(float)N
for(k0kltNk) xk0.0 yk0.0
wk(float)ke for(L0 LltN L)
ccos(Lwk) ssin(Lwk)
xk xk xrLc xiLs
yk yk xiLc - xrLs
if(Invers1) xkxk/(float)N
ykyk/(float)N
13

• Fast Fourier Transform (FFT)
• An FFT algorithm, programmed in "C ", is
  presented in Table 2.
• Gives exactly same results as DFT only much
  faster.
• Its detail how speed is achieved is outside
  our syllabus.
• We are interested in how to use DFT interpret
  its results.
• MATLAB has efficient fft procedure in its SP
  tool-box'.
• Dont need to know how its programmed, only how
  to use it!
• Direct DFT programs of academic interest only.
• Table 4 is MATLAB program which reads music from
  a 'wav' file, splits it up into 512 sample
  sections
• performs a DFT (by FFT) analysis on each
  section.

14

• Effect of windowing (old)
• DFT of x0, x1, ..., xN-1 is
  frequency-sampled version of DTFT of infinite
  sequence xn with all samples outside range
  n 0 to N-1 set to zero.
• xn effectively multiplied by "rectangular
  window" rn

When N is even, the DTFT, R(ej?) of rn is
Dirichlet Kernel
15

• Effect of windowing (newer)
• Given xn0,N-1 ? x0, x1, ..., xN-1
• Assumed obtained by multiplying infinite
  sequence xn by "rectangular window"
  rn

When N is even, the DTFT, R(ej?) of rn is
Dirichlet Kernel
16

• Effect of windowing
• Given xn0,N-1 ? x0, x1, ..., xN-1
• Assumed obtained by multiplying infinite
  sequence xn by "rectangular window"
  rn

rn
17
When N is even, the DTFT, R(ej?) of rn is
Dirichlet Kernel
18
sincs10(?/(2?)) plotted against ?
?
-?
19
-?
?
20

• Notation sincsM(?/(2?)) is not in textbooks.
• This function is also encountered when designing
  FIR filters.
• Causes the stop-band ripples which appear when
  FIR low-pass filters are designed with
  rectangular windows.
• Magnitude of R(ej ? ) shown above when M10.
• Note relatively narrow main lobe side-lobes.
• Zero-crossings occur at ? ?2? / M , ?4? / M,
  etc.
• Like a "sinc" function in many ways.

21

• Frequency-domain convolution
• DTFT of product of xn rn is complex (
  freq.-domain) convolution between X(ej?)
  R(ej?) denoted X(ej?) ? R(ej?)

• Observe the form of this expression its
  similarity with time-domain convolution
  formulae.
• Also note the (1/2?) factor limits -? to ?
  which may remind us of the inverse DTFT.

22
Duality of time- frequency-domain convolution
(i) Fourier transform of h(t) ? x(t) is
H(j?).X(j?). (ii) DTFT of hn? xn is
H(e j ? ).X(e j ? ) Time-domain
Frequency-domain Convolution (filtering)
Multiplication Multiplication (windowing)
Complex convolution
23
Duality of time- frequency-domain convolution

Time-domain Frequency-domain Convolution
(filtering) Multiplication Multiplication (
windowing) Complex convolution
24

• Spectral analysis of sampled sine-waves
• DTFT of cos(?0n) cannot be found directly.
• But inverse DTFT of
• X(e j ? ) ? ?(? - ?0) ? ?(? ?0) is

? ? (1/2?)? ( ??(?-?0) ??(??0) ) ej?n
d? ??? (1/2) e j ?o n e -j
?o n cos (? 0 n) for -? lt n lt ? It
may be inferred that DTFT of cos(?0n)
is ? ?(? - ?0) ? ?(? ?0)
25
(No Transcript)
26

• DTFT of rectangularly windowed sampled sine-wave
• DTFT of cos?0n0,N-1 DTFT of xn.rn
• X(ej?) ? R(ej?) P(ej?) say.
• By frequency-domain convolution formula,

?
27
?
?
X(e j? )
?
?0
-?0
?
-?
R(e j(? ? ?) )
?
28
X(e j ?) ? ?(?? - ?0) ? ?( ? ?0)
Now draw this its modulus
29
P(e j ? )
-? 0
? 0
30
P(e j ? )
? 0
? 0
31

• When rectangular window is of width N samples
• Ampl of main peak ampl of sine-waveN/2 check!
• N-1 zero-crossings between 0 and ?.
• 1 main peak N-2 ripples in magnitude spectrum.
• Increasing N gives sharper peak more ripples.

32

• Effect of windowing on DFT
• Effect on DTFT is frequency spreading side
  lobes.
• How does windowing with frequency sampling
  affect the DFT?
• Effect is rather confusing
• Consider 64 pt DFT of cos(0.7363n) in Fig 1.
• ?0 0.7363 lies between ?7 (2?/64)7
  ??8  (2?/64)8
• Samples of rectangular window seen.
• X7 and X8 strongly affected by sinusoid.

33
Amplitude ? 20
Fig 1 Magnitude of 64 pt DFT spectrum of
cos(0.7363n)
34

• Now consider 64 pt DFT of cos (n?/4) in Fig
  2
• ??0 ?/4 coincides exactly with ??8.
• Only X8 affected.

35
Ampl. 32 Approx. 36 difference in ampl. For
same ampl. of sinusoid.
Fig2 Magnitude of 64 point DFT of cos(pn/4).
36

• What has happened to the sincs function in
  Figure 2 ?
• We dont see it because all the samples apart
  from one
• occur at the zero-crossings of the sincs
  function.
• Effect of the rectangular window is no longer
  seen because the sampling points happen to
  coincide exactly with the zero-crossings of the
  sincs function.
• This always occurs when the frequency of the
  cosine wave coincides exactly with a frequency
  sampling point.

37
Difference between Figs 1 2 undesirable. Use
non-rectangular windows e.g. Hann wn0,N-1
with
(Slightly different definition from the one we
had for FIR filters).
38
(No Transcript)
39

• In frequency-domain
• (i) broader main-lobe for Hann window whose
  width is approx doubled as
  compared with rectr window
• (ii) reduced side-lobe levels.
• For sine-wave, at least 3 spectral 'bins'
  strongly affected even when
  its frequency coincides with a sampling point.
• Often take highest of the 3 peaks as amplitude
  of sine-wave.
• With rectr just one bin affected when
  frequency of sine-wave
  coincides with a bin. If frequency then shifted
  towards next bin, ? 35
  variation amplitude seen.
• With Hann Variation of only about 15 seen.
• Still undesirable, but
  not as bad.

40

• Reduction in amplitude estimation error for
  sine-waves is at expense of some loss of
  spectral resolution
• With Hann window, 3 bins strongly affected by
  one sine-wave
• We will only know that the frequency of the
  sine-wave lies within the range of these 3
  frequencies.

41

• DFT FFT have many applications in DSP esp in
  comms.
• e.g. filtering by (1) performing an FFT,
• (2) zeroing unwanted
  spectral components
• (3) performing an
  inverse FFT.
• FFT is Swiss army knife' of signal processing.
  
• Most spectrum analysers use an FFT algorithm.
  
• Some applications of spectral estimation are
• determining frequency loudness of a musical
  note,
• deciding whether a signal is periodic or
  non-periodic,
• analysing speech to recognise spoken words,
• looking for sine-waves buried in noise,
• measuring frequency distribution of power or
  energy.

42

• Spectral analysis of 'power signals' xn
• These exist for all time.
• Extract segment xn0,N-1 to represent xn.
  
• Energy of xn0,N-1 is

• Mean square value (MSV) of this segment is
  (1/N) E
• This is power of a periodic discrete time
  signal, of period N samples, for which a
  single cycle is xn0,N-1 .
• It may be used as an estimate of the power of
  xn.

43
Parseval's Theorem for the DFT It may be shown
that for a real signal segment xn0,N-1
Proof
44

• Parsevals Theorem allows energy power
  estimates to be calculated in
  frequency-domain instead of in time-domain.
• Allows us to see how power is distributed in
  frequency.
• Is there more power at high frequencies than at
  low frequencies or vice-versa?
• By Parseval's thm, estimate of power of xn,
  obtained by analysing xn0,N-1 , is

• Usefulness of this estimate illustrated by
  following example.

45
Example Real periodic signal xn is rect
windowed to give xn0,39 . 40-point DFT gives
magnitude spectrum below. Estimate power of
xn comment on reliability of estimate. If
xn is passed thro ideal digital low-pass
filter with cut-off ?/2 radians/sample how is
power likely to be affected?
46

• Ans
• MSV of xn0,39 ? power of xn
• (1/1600)2402 2302
  22022102 3.75 watts.
• Reduces to 3.125 watts, i.e. by 0.8 dB
• Care needed to interpret such results as power
  estimates.
• For periodic or deterministic (non-random)
  signals
• estimates from segments extracted from
  different parts of xn may be similar,
  estimates could be fairly reliable.
• For random signals may be considerable
  variation from estimate to estimate.
  Averaging may be necessary.

47

• Power spectral density (PSD) estimate
• For N-point DFT, ?Xk?2 /
  N2 is estimate of PSD in Watts per bin.
• A bin is a band-width 2?/N radians/sample
  centred on ?k
• ( fS / N
  Hz centred on k fS / N )
• Instead of Xk often plot 10 log10
  (Xk2/N2) dB. against k.
• PSD estimate graph.
• Careful with random signals each spectral
  estimate different.
• Can average several PSD estimates.

48

• Spectral analysis of signals
• File OPERATOR.pcm contains sampled speech.
• SNR-12dB.pcm contains sine-wave
  corrupted with noise.
• Sampled at 8 kHz using 12-bit A/D converter.
• May be read into "MATLAB" program in Table 5
  spectrally
  analysed using the FFT.
• Meaningless to analyse a large speech file at
  once.
• Divide into blocks of 128 or 256 samples
  analyse separately.
• Blocks of N ( 512) samples may be read in and
  displayed.
• Programs in notes available on
  barry/mydocs/CS3291
• In each case, a rectangular window is used.
• Exercise, modify programs to have a Hann window.
  
• Notice effectiveness of DFT in locating
  sine-wave in noise
• even when it cannot be seen
  in time-domain graph.