Computer Science 121

1
Computer Science 121

• Scientific Computing
• Winter 2007
• Advanced Topics
• Digital Signal Processing

2
DSP I The Fourier Transform

• Weve seen that a signal can be synthesized from
  a set of component frequencies and associated
  amplitudes (and phases)

• What if we already have a signal and want to
  analyze it to get the amplitude Ai (and phase fi)
  for each frequency fi ?
• This (without phase?) is essentially what the
  cochlea (inner ear) does
• So its useful for speech recognition, noise
  filtering, and related signal-processing tasks

3
DSP I The Fourier Transform Confusing
Mathematical Background

• Fourier (ca. 1800) showed that any function
  (signal) can be represented as a sum of sines and
  cosines

• I.e., we have an infinite number of frequencies
  n, and want to get the amplitudes an, bn for
  each.
• Use integration

4

• Keeping the sine and cosine components separate
  allows us to preserve phase information. We can
  do this by using complex numbers

• Euler (ca. 1750) showed that cos(x) i sin(x)
  eix. So we can rewrite the Fourier series as

5

• This complex form is also nice because it gives
  us just one set of coefficients cn, which we can
  obtain by integration

6

• We now have a Fourier Transform that we can use
  to determine the amount cn that each frequency
  contributes to our signal g(x)

• Of course, we cant look at an infinite number of
  frequencies, and we (probably) dont know how to
  integrate the signal analytically (as we could
  with, e.g., g(x) x2). So we replace our
  continuous functions and transforms by discrete
  versions

7
The Discrete Fourier Transform

• where N is the length of our signal vector

8
The Discrete Fourier Transform Some Cool Results

• For N sample values, we get N frequencies the
  signal and its DFT are exactly the same thing,
  expressed in different domains (time vs.
  frequency).
• Information about each sample is smeared
  (convolved) across the entire DFT.
• The component ei2pnk/N (eip)2nk/N is an almost
  magical number, containing all three
  fundamental mathematical constants.

9
The Fast Fourier Transform

• With N values to look at for each of the N
  frequencies, the DFT takes N2 time to compute.
• Cooley Tukey (1965 c.f. Gauss 1805) invented a
  Fast Fourier Transform using divide-and-conquer
  to compute it in N log2(N) time (what does this
  remind us of?)
• The FFT is just a faster version of the DFT,
  producing the same output in much less time
• FFT is what everyone uses in practice.

10
FFT Dealing with Imaginary Numbers

• The DFT/FFT produces amplitude coefficients of
  the form a ib essentially, a magnitude, b
  phase
• To visualize the transformed signal, we want real
  (not complex) numbers
• Two options
• Throw out the imaginary part a real(fft(x))
• Convert amplitude to power

11
FFT Converting Amplitude to Power

• Every complex number a ib has a conjugate a -
  ib
• Multiplying by conjugate always yields a positive
  real (a ib) (a - ib) (a2
  - i2b2) (a2 b2)
• This operation gives us the power at the
  corresponding frequency
• In Matlab p fft(x) . conj(fft(x))
• p is called the power spectrum of the signal

12
FFT DEMO
gtgt FS 10000 gtgt t linspace(0,1,FS) gtgt Pt1
sin(2pi500t) gtgt Pt2 sin(2pi1000t) gtgt Pt
Pt1 Pt2 gtgt plot(fft(Pt).conj(fft(Pt)))
13
FFT Applications

• Possibly the most widely used numerical algorithm
  of all time
• DSP (speech, seismography, )
• Image processing (PhotoShop)
• Solving partial differential equations
• Multiplying large numbers
• Holography
• Allows us to visualize the frequencies that make
  up a signal
• Visualization is often easier if we plot the log
  of the power instead of its actual value.
• E.g., what can happen if we record in a room with
  noisy fluorescent lights (60 Hz buzz)?

14
Visualizing the Effect of Noise

• gtgt x,fs wavread('g-fugue-mono.wav')
• gtgt sound(x,fs)
• gtgt dur length(x)/fs
• gtgt nsamp fsdur
• gtgt t linspace(0,dur,nsamp)
• gtgt hum sin(2pi60t)' 60 Hz buzz
• gtgt sound(hum,fs)
• gtgt xh x hum
• gtgt sound(xh, fs)
• gtgt continues

15
Visualizing the Effect of Noise

• gtgt fx fft(x) clean recording
• gtgt sound(ifft(fx), fs) ifft(fft(x)) x
• gtgt fxh fft(xh) recorded in present of 60Hz
  noise
• gtgt plot(log(fx.conj(fx))) log Power plot
• gtgt figure
• gtgt plot(log(fxh.conj(fxh)))
• gtgt figure(1), xlim(1 1000), ylim(0 25)
• gtgt figure(2), xlim(1 1000), ylim(0 25)

16
DSP II Digital Filters

• Real-world signals will often contain components
  that we want to remove, or filter out, before
  using the signal - e.g., background noise in an
  audio recording
• Or we may want to emphasize (filter in) certain
  aspects of a signal or image at the expense of
  others e.g., show only edges of objects in an
  image
• To see how we can do this, it is easiest to start
  with a familiar idea

17
Averaging as Low-Pass Filtering

• Q .What information do we lose/discard when we
  take an average (mean) of a set of numbers?
• A .We discard information about individual
  differences.
• Replacing a signal value with the mean of itself
  and its neighbors likewise smoothes out the
  differences

18
Averaging as Low-Pass Filtering

• Extreme points become less extreme

19
Averaging as Low-Pass Filtering

• Repeating this for all the samples gives us a
  smoothed signal

gtgt FS 10000 gtgt t linspace(0,1,FS) gtgt Pt
sin(2pi440t) gtgt Pt Pt .05
randn(size(Pt)) gtgt plot(t, Pt), xlim(0 .01)
gtgt Pt2 Pt necessary? gtgt for i
3length(Pt)-2, ... Pt2(i)
mean(Pt(i-2i2)) end gtgt plot(t, Pt2),
xlim(0 .01)
20
Averaging as Low-Pass Filtering

• Smooth also means slowly-changing, or lower
  frequency
• Rough also means quickly changing, or higher
  frequency
• So smoothing is the same as removing
  higher-frequency blips in the signal, allowing
  the lower frequencies to pass through unmodified.
• So averaging/smoothing is also called low-pass
  filtering.

21
Averaging as Low-Pass Filtering
gtgt FS 10000 gtgt t linspace(0,1,FS) gtgt Pt1
sin(2pi440t) gtgt Pt2 sin(2pi1000t) gtgt Pt
Pt1 Pt2 gtgt sound(Pt, FS) gtgt
plot(fft(Pt).conj(fft(Pt))) gtgt xlim(0 5000)
gtgt Pt3 Pt gtgt for i 5length(Pt3)-4,
Pt3(i) mean(Pt(i-4i4)) end gtgt
sound(Pt3, FS) gtgt plot(fft(Pt3).conj(fft(Pt3))) gt
gt xlim(0 5000), ylim(0 2.5e7)
22
First-Differencing as High-Pass Filtering

• Instead of smoothing out the successive
  differences in signal values, we can keep the
  differences and discard the actual values.
• This called first-differencing and is the
  discrete-time version of the first derivative
  from the calculus.
• I.e., we can keep the rough part and throw away
  the smooth
• This produces a kind of high-pass filter.

23
First-Differencing as High-Pass Filtering
gtgt FS 10000 gtgt t linspace(0,1,FS) gtgt Pt
sin(2pi440t) gtgt Pt Pt .05
randn(size(Pt)) gtgt plot(t, Pt), xlim(0 .01)
gtgt Pt2 Pt gtgt for i 2length(Pt), ...
Pt2(i) Pt(i) Pt(i-1) end gtgt plot(t,
Pt2) gtgt xlim(0 .01), ylim(-1 1)
24
First-Differencing as High-Pass Filtering
gtgt plot(log(fft(Pt).conj(fft(Pt)))) gtgt xlim(0
5000), ylim(-20 20)
gtgt plot(log(fft(Pt2).conj(fft(Pt2)))) gtgt xlim(0
5000), ylim(-20 20)
25
Filter Coefficients

• Wed like a general mathematical form for all
  such filters.
• Given a signal x, we can express the filtered
  version y as

• Where the ai are the coefficients of the filter

26

• In this form, the neighborhood includes just
  the previous M -1 samples and current sample (not
  the next couple of samples), but it doesnt
  really make a difference.
• E.g., our low-pass filter has ai 1/M for all i
• Q. What are the M and ai for our
  first-differencing high-pass filter?

27
Can we avoid looping?
Pt2 Pt for i 3length(Pt)-2, ... Pt2(i)
mean(Pt(i-2i2)) end
Pt2 Pt Pt2(2end-1) 1/3Pt(1end-2) ...
1/3Pt(2end-1) ...
1/3Pt(3end)
Pt2 zeros(size(Pt)) for i 13 Pt2(2end-1)
Pt2(2end-1) 1/3Pt(iend-(3-i)) end
28
Filtering in Matlab

• You shouldnt have to design your own digital
  filters from scratch use Matlabs
  signal-processing toolkit
• http//www.mathworks.com/access/helpdesk/help/too
  lbox/signal
• gtgt help signal

• Idea is to specify desired frequency
  characteristics (high-pass, low-pass) and have
  Matlab figure out the filter coefficients.

29
Example Edge Detection

• Recall that we can describe an image as a
  two-dimensional signal

P(x, y) Sj,k Aj,k sin(2p( fj x fk y fj,k))

• Edges (lines between objects) in the image
  correspond to high-frequency components thats
  where the signal (image) changes most rapidly
• So a high-pass image filter acts as an edge
  detector

30
Example Edge Detection

• Matlabs image-processing toolkit has a built-in
  edge function to do this for us

http//www.mathworks.com/access/helpdesk/help/too
lbox/images gtgt help images gtgt help edge