Moving into the Frequency Domain

CM0268
Frequency domains can be obtained through the MATLAB
DSP
GRAPHICS
transformation from one (Time or Spatial) domain to the other
261
(Frequency) via Fourier Transform (FT)

Fourier Transform (FT) MPEG Audio

Related Discrete Cosine Transform (DCT) Heart of JPEG and
MPEG Video, (alt.) MPEG Audio.
Not Studied here CM0340 Multimedia (YEAR 3) 1

JJ
II
J
I
Back
Close
1D Example
Lets consider a 1D (e.g. Audio) example to see what the CM0268
MATLAB
different domains mean: DSP
GRAPHICS

262
Consider a complicated sound such as the noise of a car
horn. We can describe this sound in two related ways:
Sample the amplitude of the sound many times a second,
which gives an approximation to the sound as a function of
time.
1

Analyse the sound in terms of the pitches of the notes, or

frequencies, which make the sound up, recording the JJ
amplitude of each frequency. II
J
I
Back
Close
An 8 Hz Sine Wave
In the example (next slide): CM0268
MATLAB
DSP
A signal that consists of a sinusoidal wave at 8 Hz. GRAPHICS

263
8 Hz means that wave is completing 8 cycles in 1 second
The frequency of that wave (8 Hz).
From the frequency domain we can see that the composition
of our signal is
one wave (one peak) occurring with a frequency of 8 Hz 1

with a magnitude/fraction of 1.0 i.e. it is the whole signal.

JJ
II
J
I
Back
Close
An 8 Hz Sine Wave (Cont.)

CM0268
MATLAB
DSP
GRAPHICS

264

JJ
II
J
I
Back
Close
2D Image Example
Now images are no more complex really: CM0268
MATLAB
DSP
Brightness along a line can be recorded as a set of values GRAPHICS

265
measured at equally spaced distances apart,
Or equivalently, at a set of spatial frequency values.
Each of these frequency values is a frequency component.
An image is a 2D array of pixel measurements.
We form a 2D grid of spatial frequencies. 1

A given frequency component now specifies what contribution

is made by data which is changing with specified x and y JJ
direction spatial frequencies. II
J
I
Back
Close
What do frequencies mean in an image?
Large values at high frequency components then the data CM0268
MATLAB
is changing rapidly on a short distance scale. DSP
GRAPHICS

e.g. a page of text 266

Large low frequency components then the large scale features

of the picture are more important.
e.g. a single fairly simple object which occupies most of the
image.
1

JJ
II
J
I
Back
Close
How to Filter?
Low pass filter CM0268
MATLAB
DSP
Ignore high frequency noise components make zero or GRAPHICS

a very low value. 267

Only store lower frequency components

JJ
II
J
I
Back
Close
Visualising Frequency Domain Transforms
Any function (signal) can be decomposed into purely CM0268
MATLAB
sinusoidal components (sine waves of different size/shape) DSP
GRAPHICS

When added together make up our original signal. 268

Fourier transform is the tool that performs such an operation

JJ
II
J
I
Back
Close
Summing Sine Waves
Digital signals are composite signals made up of many sinusoidal CM0268

frequencies MATLAB
DSP
GRAPHICS

269

JJ
II
J
I
Back
Close
Summing Sine Waves to give a Square(ish)
Wave CM0268
MATLAB
DSP
We can take the previous example a step further: GRAPHICS

270
A 200Hz digital signal (square(ish) wave) may be a composed
of 200, 600, 1000, 1400, 1800, 2200, 2600, 3000, 3400 and
3800 sinusoidal signals which sum to give:

JJ
II
J
I
Back
Close
So What Does All This Mean?
Transforming a signal into the frequency domain allows us CM0268
MATLAB
DSP
To see what sine waves make up our underlying signal GRAPHICS

271
E.g.
One part sinusoidal wave at 50 Hz and
Second part sinusoidal wave at 200 Hz.
More complex signals will give more complex graphs but the
idea is exactly the same. 1

Filtering now involves attenuating or removing certain

frequencies easily performed. JJ
II
The graph of the frequency domain is called the frequency J
spectrum more soon I
Back
Close
Visualising Frequency Domain:
Think Graphic Equaliser CM0268
MATLAB
DSP
GRAPHICS

An easy way to visualise what is happening is to think of a 272

graphic equaliser on a stereo system (or some software audio

players, e.g. iTunes).

JJ
II
J
I
Back
Close
Fourier Theory
The tool which converts a spatial (real space) description of
audio/image data into one in terms of its frequency components
CM0268
is called the Fourier transform MATLAB
DSP
GRAPHICS

273
The new version is usually referred to as the Fourier space
description of the data.
We then essentially process the data:
E.g. for filtering basically this means attenuating or setting
certain frequencies to zero
1

We then need to convert data back to real audio/imagery to use

in our applications.
JJ
II
The corresponding inverse transformation which turns a Fourier J
space description back into a real space one is called the I
inverse Fourier transform. Back
Close
Fourier Transform
1D Case (e.g. Audio Signal) CM0268
MATLAB
DSP
GRAPHICS
Considering a continuous function f (x) of a single variable x
274
representing distance.
The Fourier transform of that function is denoted F (u), where
u represents spatial frequency is defined by
Z
F (u) = f (x)e2ixu dx.

1

Note: In general F (u) will be a complex quantity even though

the original data is purely real. JJ
The meaning of this is that not only is the magnitude of each II
J
frequency present important, but that its phase relationship is I
too. Back
Close
Inverse 1D Fourier Transform

The inverse Fourier transform for regenerating f (x) from F (u)

CM0268
is given by MATLAB
DSP
GRAPHICS

275
Z
f (x) = F (u)e2ixu du,

which is rather similar, except that theexponential term has

the opposite sign. not negative
1

JJ
II
J
I
Back
Close
Example Fourier Transform
Lets see how we compute a Fourier Transform: consider a CM0268
MATLAB
particular function f (x) defined as DSP
GRAPHICS

276

1 if |x| 1
f (x) =
0 otherwise,

1
JJ
II
1 J
I
Back
Close
 So its Fourier transform is:
Z
F (u) = f (x)e2ixu dx
Z1 CM0268

1 e2ixu dx
MATLAB
= DSP
GRAPHICS

1
277
1 2iu
= (e e2iu)
2iu
sin 2u
= .
u
In this case F (u) is purely real, which is a consequence of the
original data being symmetric in x and x.
1

A graph of F (u) is shown overleaf. JJ

II
This function is often referred to as the Sinc function. J
I
Back
Close
The Sync Function
The Fourier transform of a top hat function:
CM0268
MATLAB
DSP
GRAPHICS

278

JJ
II
J
I
Back
Close
2D Case (e.g. Image data)
CM0268
MATLAB
If f (x, y) is a function, for example the brightness in an image, DSP
GRAPHICS

its Fourier transform is given by 279

Z Z
F (u, v) = f (x, y)e2i(xu+yv) dx dy,

and the inverse transform, as might be expected, is

1

Z Z
f (x, y) = F (u, v)e2i(xu+yv) du dv. JJ

II
J
I
Back
Close
But All Our Audio and Image data are
Digitised!! CM0268
MATLAB
DSP
Thus, we need a discrete formulation of the Fourier transform: GRAPHICS

280
Which takes such regularly spaced data values, and
Returns the value of the Fourier transform for a set of values
in frequency space which are equally spaced.

This is done quite naturally by replacing the integral by a 1

summation, to give the discrete Fourier transform or DFT for

short. JJ
II
In 1D it is convenient now to assume that x goes up in steps J
of 1, and that there are N samples, at values of x from 0 to N 1. I
Back
Close
1D Discrete Fourier transform
So the DFT takes the form CM0268
MATLAB
N 1 DSP
GRAPHICS
1 X
F (u) = f (x)e2ixu/N , 281
N x=0

while the inverse DFT is

N
X 1
f (x) = F (u)e2ixu/N .
x=0
1

NOTE: Minor changes from the continuous case are a factor JJ

of 1/N in the exponential terms, and also the factor 1/N in front II
J
of the forward transform which does not appear in the inverse
I
transform. Back
Close
2D Discrete Fourier transform
CM0268
MATLAB
The 2D DFT works is similar. So for an N M grid in x and y DSP
GRAPHICS

we have 282

N 1 M 1
1 XX
F (u, v) = f (x, y)e2i(xu/N +yv/M ),
N M x=0 y=0

and
N
X 1 M
X 1
f (x, y) = F (u, v)e2i(xu/N +yv/M ).
1

u=0 v=0
JJ
II
J
I
Back
Close
Balancing the 2D DFT
Often N = M , and it is then it is more convenient to redefine CM0268
MATLAB
F (u, v) by multiplying it by a factor of N , so that the forward and DSP
GRAPHICS

inverse transforms are more symmetrical: 283

N 1 N 1
1 XX
F (u, v) = f (x, y)e2i(xu+yv)/N ,
N x=0 y=0
and
1

N 1 N 1 JJ
1 XX II
f (x, y) = F (u, v)e2i(xu+yv)/N .
N u=0 v=0 J
I
Back
Close
Visualising the Fourier Transform
CM0268
MATLAB
DSP
Its useful to visualise the Fourier Transform GRAPHICS

284
Standard tools
Easily plotted in MATLAB

JJ
II
J
I
Back
Close
The Magnitude Spectrum of Fourier
Transform CM0268
MATLAB
DSP
Recall that the Fourier Transform of even our real audio/image GRAPHICS

data is always complex. 285

How can we visualise a complex data array?

Compute the absolute value of the complex data:
q
|F (k)| = FR2 (k) + FI2(k) for k = 0, 1, . . . , N 1

where FR (k) is the real part and FI (k) of the N sampled Fourier 1

Transform, F (k).
JJ
II
This is called the magnitude spectrum of the Fourier Transform J
Easy in MATLAB: Sp = abs(fft(X,N))/N; I
(Normalised form) Back
Close
The Phase Spectrum of the Fourier
Transform CM0268
MATLAB
DSP
GRAPHICS

The Fourier Transform also represent phase, the phase 286

spectrum is given by:

FI (k)
= arctan for k = 0, 1, . . . , N 1
FR (k)

JJ
II
J
I
Back
Close
Relating a Sample Point to a Frequency
Point CM0268
MATLAB
DSP
When plotting graphs of FT Spectra and doing other FT GRAPHICS

processing we may wish to plot the x-axis in Hz (Frequency) 287

rather than sample point number k = 0, 1, . . . , N 1

There is a simple relation between the two:

The sample points go in steps k = 0, 1, . . . , N 1
For a given sample point k the frequency relating to this is given
by: 1

fs
fk = k JJ
N
where fs is the sampling frequency and N the number of samples. II
J
Thus we have equidistant frequency steps of fNs ranging from
I
0 Hz to NN1 fs Hz Back
Close
MATLAB Fourier Frequency Spectra
Example CM0268
MATLAB
DSP
The following code (fourierspectraeg.m): GRAPHICS

N=16; 288
FS=40000;
x=cos(2*pi*2*(0:1:N-1)/N);
f=((0:N-1)/N)*FS;
X=abs(fft(x,N))/N;
figure(1)
subplot(3,1,3);plot(f,X);
subplot(3,1,1);stem(0:N-1,x,.);
axis([-0.2*44100/16 max(f) -0.1 1.1]);
axis([-0.2 N -1.2 1.2]);
legend(Magnitude spectrum |X(f)|);
legend(Cosine signal x(n));
ylabel(c));
ylabel(a));
xlabel(f in Hz \rightarrow)
xlabel(n \rightarrow);
figure(2)
X=abs(fft(x,N))/N;
subplot(3,1,1);
subplot(3,1,2);stem(0:N-1,X,.);
plot(f,20*log10(X./(0.5)));
axis([-0.2 N -0.1 1.1]); 1

axis([-0.2*44100/16 max(f) ...

legend(Magnitude spectrum |X(k)|);
-45 20]);
ylabel(b));
legend(Magnitude spectrum |X(f)| ...
xlabel(k \rightarrow)
in dB); JJ
ylabel(|X(f)| in dB \rightarrow);
N=1024;
x=cos(2*pi*(2*1024/16)*(0:1:N-1)/N);
xlabel(f in Hz \rightarrow) II
J
I
Back
Close
MATLAB Fourier Frequency Spectra Example (Cont.)
The above code produces the following:
1
Cosine signal x(n) CM0268
MATLAB
0 DSP
a)

GRAPHICS

1 289
0 2 4 6 8 10 12 14 16
n

1
Magnitude spectrum |X(k)|

0.5
b)

0
0 2 4 6 8 10 12 14 16
k 1

1
Magnitude spectrum |X(f)|
JJ
0.5
c)

II
0 J
0 0.5 1 1.5 2 2.5 3 3.5
f in Hz x 10
4 I
Back
Close
Magnitude Spectrum in dB
CM0268
MATLAB
Note: It is common to plot both spectra magnitude (also DSP
GRAPHICS

frequency ranges not show here) on a dB/log scale: 290

(Last Plot in above code)

20
|X(f)| in dB

Magnitude spectrum |X(f)| in dB

0

20

40
1

0 0.5 1 1.5 2 2.5 3 3.5

f in Hz 4
x 10

JJ
II
J
I
Back
Close
Time-Frequency Representation:
Spectrogram CM0268
MATLAB
DSP
It is often useful to look at the frequency distribution over a GRAPHICS

short-time: 291

Split signal into N segments

Do a windowed Fourier Transform
Window needed to reduce leakage effect of doing a short
sample FFT.
1

Apply a Blackman, Hamming or Hanning Window

MATLAB function does the job: Spectrogram see help JJ
spectrogram II
J
See also MATLABs specgramdemo I
Back
Close
MATLAB Example
The code: CM0268
MATLAB
load(handel) DSP
GRAPHICS
[N M] = size(y);
figure(1) 292
spectrogram(fft(y,N),512,20,1024,Fs);

Produces the following:

JJ
II
J
I
Back
Close
Filtering in the Frequency Domain
Low Pass Filter CM0268
MATLAB
DSP
Example: Frequencies above the Nyquist Limit, GRAPHICS

293
Noise:
The idea with noise Filtering is to reduce various spurious effects
of a local nature in the image, caused perhaps by
noise in the acquisition system,
arising as a result of transmission of the data, for example
from a space probe utilising a low-power transmitter.
1

JJ
II
J
I
Back
Close
Frequency Space Filtering Methods
CM0268
Noise = High Frequencies: MATLAB
DSP
GRAPHICS

In audio data many spurious peaks in over a short timescale. 294

In an image means there are many rapid transitions (over a short

distance) in intensity from high to low and back again or vice
versa, as faulty pixels are encountered.
Not all high frequency data noise though!

Therefore noise will contribute heavily to the high frequency

components of the image when it is considered in Fourier space.
JJ
II
Thus if we reduce the high frequency components Low-Pass
J
Filter, we should reduce the amount of noise in the data.
I
Back
Close
(Low-pass) Filtering in the Fourier Space
We thus create a new version of the image in Fourier space by
CM0268
computing MATLAB
DSP
GRAPHICS
G(u, v) = H(u, v)F (u, v)
295
where:
F (u, v) is the Fourier transform of the original image,
H(u, v) is a filter function, designed to reduce high frequencies,
and
G(u, v) is the Fourier transform of the improved image.
1

Inverse Fourier transform G(u, v) to get g(x, y) our improved

image
JJ
II
Note: Discrete Cosine Transform approach identical, sub. FT with J
DCT I
Back
Close
Ideal Low-Pass Filter
The simplest sort of filter to use is an ideal low-pass filter, which in
CM0268
one dimension appears as : MATLAB
DSP
GRAPHICS

2.0 296
H(u)

JJ
II
u0 u
J
I
Back
Close
Ideal Low-Pass Filter (Cont.)
2.0
H(u)

CM0268
MATLAB
DSP
GRAPHICS

297

u0 u

This is a top hat function which is 1 for u between 0 and u0 , the cut-off
1

frequency, and zero elsewhere.

JJ
So All frequency space space information above u0 is thrown II
away, and all information below u0 is kept. J
A very simple computational process. I
Back
Close
Ideal 2D Low-Pass Filter
The two dimensional analogue of this is the function
CM0268

 MATLAB
DSP
1 if u2 + v 2 w0 GRAPHICS

H(u, v) = 298
0 otherwise,
where w0 is now the cut-off frequency.

Thus, all frequencies inside a radius w0 are kept, and all others
discarded.

w0

JJ
II
J
I
Back
Close
Not So Ideal Low-Pass Filter?
CM0268
The problem with this filter is that as well as the noise: MATLAB
DSP
GRAPHICS

In audio: plenty of other high frequency content 299

In Images: edges (places of rapid transition from light to dark)

also significantly contribute to the high frequency
components.

Thus an ideal low-pass filter will tend to blur the data:

1

High audio frequencies become muffled

Edges in images become blurred. JJ
II
The lower the cut-off frequency is made, the more pronounced this
J
effect becomes in useful data content
I
Back
Close
Ideal Low Pass Filter Example 1

CM0268
MATLAB
DSP
GRAPHICS

300

(a) Input Image (b) Image Spectra

JJ
II
J
(c) Ideal Low Pass Filter (d) Filtered Image I
Back
Close
Ideal Low-Pass Filter Example 1 MATLAB Code
low pass.m:
% Create a white box on a black background image
CM0268
M = 256; N = 256; MATLAB
image = zeros(M,N) DSP
GRAPHICS
box = ones(64,64);
%box at centre 301
image(97:160,97:160) = box;

% Show Image

figure(1);
imshow(image);

% compute fft and display its spectra

1

F=fft2(double(image));
figure(2);
imshow(abs(fftshift(F))); JJ
II
J
I
Back
Close
Ideal Low-Pass Filter Example 1 MATLAB Code (Cont.)
%compute Ideal Low Pass Filter
u0 = 20; % set cut off frequency

u=0:(M-1); CM0268
MATLAB
v=0:(N-1); DSP
GRAPHICS
idx=find(u>M/2);
u(idx)=u(idx)-M; 302
idy=find(v>N/2);
v(idy)=v(idy)-N;
[V,U]=meshgrid(v,u);
D=sqrt(U.2+V.2);
H=double(D<=u0);

% display
figure(3);
imshow(fftshift(H)); 1

% Apply filter and do inverse FFT

G=H.*F; JJ
g=real(ifft2(double(G))); II
% Show Result
J
figure(4); I
imshow(g); Back
Close
Ideal Low-Pass Filter Example 2

CM0268
MATLAB
DSP
GRAPHICS

303

(a) Input Image (b) Image Spectra

JJ
II
J
(c) Ideal Low-Pass Filter (d) Filtered Image I
Back
Close
Ideal Low-Pass Filter Example 2 MATLAB Code
lowpass2.m:
% read in MATLAB demo text image
CM0268
image = imread(text.png); MATLAB
[M N] = size(image) DSP
GRAPHICS

304
% Show Image

figure(1);
imshow(image);

% compute fft and display its spectra

F=fft2(double(image));
figure(2); 1

imshow(abs(fftshift(F))/256);
JJ
II
J
I
Back
Close
Ideal Low-Pass Filter Example 2 MATLAB Code (Cont.)
%compute Ideal Low Pass Filter
u0 = 50; % set cut off frequency

u=0:(M-1); CM0268
MATLAB
v=0:(N-1); DSP
GRAPHICS
idx=find(u>M/2);
u(idx)=u(idx)-M; 305
idy=find(v>N/2);
v(idy)=v(idy)-N;
[V,U]=meshgrid(v,u);
D=sqrt(U.2+V.2);
H=double(D<=u0);

% display
figure(3);
imshow(fftshift(H)); 1

% Apply filter and do inverse FFT

G=H.*F; JJ
g=real(ifft2(double(G))); II
% Show Result J
figure(4); I
imshow(g); Back
Close
Low-Pass Butterworth Filter
CM0268
Another filter sometimes used is the Butterworth low pass filter. MATLAB
DSP
GRAPHICS

306
In the 2D case, H(u, v) takes the form

1
H(u, v) = n,
1 + [(u2 + v 2)/w02]
where n is called the order of the filter.
1

JJ
II
J
I
Back
Close
Low-Pass Butterworth Filter (Cont.)

This keeps some of the high frequency information, as illustrated

by the second order one dimensional Butterworth filter: CM0268
MATLAB
DSP
GRAPHICS

.0
307
H(u)

.0
u0 u
JJ
Consequently reduces the blurring. II
J
I
Back
Close
Low-Pass Butterworth Filter (Cont.)
The 2D second order Butterworth filter looks like this:

CM0268
MATLAB
DSP
GRAPHICS

308

w0

JJ
II
J
I
Back
Close
Butterworth Low Pass Filter Example 1

CM0268
MATLAB
DSP
GRAPHICS

309

(a) Input Image (b) Image Spectra

JJ
II
J
(c) Butterworth Low-Pass Filter (d) Filtered Image I
Back
Close
Butterworth Low-Pass Filter Example 1 (Cont.)

Comparison of Ideal and Butterworth Low Pass Filter:

CM0268
MATLAB
DSP
GRAPHICS

310

JJ
II
J
Ideal Low-Pass Butterworth Low Pass I
Back
Close
Butterworth Low-Pass Filter Example 1 MATLAB Code
butterworth.m:
% Load Image and Compute FFT as in Ideal Low Pass Filter
% Example 1 CM0268
....... MATLAB
DSP
GRAPHICS
% Compute Butterworth Low Pass Filter
u0 = 20; % set cut off frequency 311

u=0:(M-1);
v=0:(N-1);
idx=find(u>M/2);
u(idx)=u(idx)-M;
idy=find(v>N/2);
v(idy)=v(idy)-N;
[V,U]=meshgrid(v,u);
1

for i = 1: M
for j = 1:N
%Apply a 2nd order Butterworth JJ
UVw = double((U(i,j)*U(i,j) + V(i,j)*V(i,j))/(u0*u0));
H(i,j) = 1/(1 + UVw*UVw);
II
end J
end I
% Display Filter and Filtered Image as before
Back
Close
Butterworth Low-Pass Butterworth Filter Example 2

CM0268
MATLAB
DSP
GRAPHICS

312

(a) Input Image (b) Image Spectra

JJ
II
J
(c) Butterworth Low-Pass Filter (d) Filtered Image I
Back
Close
Butterworth Low-Pass Filter Example 2 (Cont.)

Comparison of Ideal and Butterworth Low-Pass Filter:

CM0268
MATLAB
DSP
GRAPHICS

313

JJ
II
J
Ideal Low Pass Butterworth Low Pass I
Back
Close
Butterworth Low Pass Filter Example 2 MATLAB Code
butterworth2.m:
% Load Image and Compute FFT as in Ideal Low Pass Filter
% Example 2 CM0268
....... MATLAB
DSP
GRAPHICS
% Compute Butterworth Low Pass Filter
u0 = 50; % set cut off frequency 314

u=0:(M-1);
v=0:(N-1);
idx=find(u>M/2);
u(idx)=u(idx)-M;
idy=find(v>N/2);
v(idy)=v(idy)-N;
[V,U]=meshgrid(v,u);
1

for i = 1: M
for j = 1:N
%Apply a 2nd order Butterworth JJ
UVw = double((U(i,j)*U(i,j) + V(i,j)*V(i,j))/(u0*u0));
H(i,j) = 1/(1 + UVw*UVw);
II
end J
end I
% Display Filter and Filtered Image as before
Back
Close
Low Pass Filtering Noisy Images
Use Matlab function, imnoise() to add noise to image (lowpass.m,
lowpass2.m): Image with Noise Added Low Pass Filtered Noisy Image

CM0268
MATLAB
DSP
GRAPHICS

315

(a) Input Noisy Image

Image with Noise Added (b) Deconvolved Noisy
High Cutoff FrequencyImage (Low
Low Pass Filtered Image Cut Off)

JJ
II
J
I
(c) Input Noisy Image (d) Deconvolved Noisy Image (Higher Cut Off) Back
Close
Other Filters
High-Pass Filters opposite of low-pass, select high
CM0268
frequencies, attenuate those below u0 MATLAB
DSP
GRAPHICS

Band-pass allow frequencies in a range u0 . . . u1 attenuate those 316

outside this range
Band-reject opposite of band-pass, attenuate frequencies within
u0 . . . u1 select those outside this range
Notch attenuate frequencies in a narrow bandwidth around cut-off
frequency, u0
Resonator amplify frequencies in a narrow bandwidth around cut-off 1

frequency, u0
JJ
II
Other filters exist that are a combination of the above
J
I
Back
Close
Convolution
Several important audio and optical effects can be described in
CM0268
terms of convolutions. MATLAB
DSP
GRAPHICS

In fact the above Fourier filtering is applying convolutions of low 317

pass filter where the equations are Fourier Transforms of real
space equivalents.
Deblurring high pass filtering
Reverb more soon.

JJ
II
J
I
Back
Close
1D Convolution
CM0268
Let us examine the concepts using 1D continuous functions. MATLAB
DSP
GRAPHICS

318
The convolution of two functions f (x) and g(x), written f (x) g(x),
is defined by the integral
Z
f (x) g(x) = f ()g(x ) d.

JJ
II
J
I
Back
Close
1D Convolution Example
For example, let us take two top hat functions of the type described
CM0268
earlier. MATLAB
DSP
GRAPHICS

319
Let f () be the top hat function shown:

1 if || 1
f () =
0 otherwise,
and let g() be as shown in next slide, defined by

1/2 if 0 1
g() = 1

0 otherwise.

JJ
II
J
I
Back
Close
1D Convolution Example (Cont.)

1.0 1.0
CM0268
MATLAB
DSP
GRAPHICS

320

-5.0 0.0 5.0 -5.0 0.0 5.0

  JJ
1 if || 1 1/2 if 0 1 II
f () = g() =
0 otherwise, 0 otherwise. J
I
Back
Close
1D Convolution Example (Cont.)

g() is the reflection of CM0268

1.0 MATLAB
this function in the vertical DSP
GRAPHICS

axis, 321
g(x ) is the latter shifted
to the right by a distance x.
Thus for a given value of
x, f ()g(x ) integrated
over all is the area of
overlap of these two top 1

hats, as f () has unit

height.
JJ
An example is shown for x II
in the range 1 x 0 J
opposite I
-5.0 x 0.0 5.0 Back
Close
1D Convolution Example (cont.)

If we now consider x moving from to +, we can see that

CM0268
For x 1 or x 2, there is no overlap; MATLAB
DSP
GRAPHICS

As x goes from 1 to 0 the area of overlap steadily increases 322

from 0 to 1/2;
As x increases from 0 to 1, the overlap area remains at 1/2;
Finally as x increases from 1 to 2, the overlap area steadily
decreases again from 1/2 to 0.
Thus the convolution of f (x) and g(x), f (x) g(x), in this case
has the form shown on next slide 1

JJ
II
J
I
Back
Close
1D Convolution Example (cont.)

1.0
CM0268
MATLAB
DSP
GRAPHICS

323

JJ
-5.0 0.0 5.0 II
J
Result of f (x) g(x) I
Back
Close
1D Convolution Example (cont.)

Mathematically the convolution is expressed by:

CM0268
MATLAB
DSP
GRAPHICS

(x + 1)/2 if 1 x 0 324
1/2 if 0 x 1

f (x) g(x) =

1 x/2 if 1 x 2
0 otherwise.

1.0

JJ
II
J
I
-5.0 0.0 5.0 Back
Close
Fourier Transforms and Convolutions
CM0268
One major reason that Fourier transforms are so important in signal/image
MATLAB
DSP
GRAPHICS
processing is the convolution theorem which states that:
325

If f (x) and g(x) are two functions with Fourier transforms F (u)
and G(u), then the Fourier transform of the convolution f (x) g(x)
is simply the product of the Fourier transforms of the two functions,
F (u)G(u).

Recall our Low Pass Filter Example (MATLAB CODE) 1

% Apply filter JJ
G=H.*F; II
Where F was the Fourier transform of the image, H the filter J
I
Back
Close
Computing Convolutions with the Fourier
Transform CM0268
MATLAB
E.g.: DSP
GRAPHICS

To apply some reverb to an audio signal, example later 326

To compensate for a less than ideal image capture system:

To do this fast convolution we simply:
Take the Fourier transform of the audio/imperfect image,
Take the Fourier transform of the function describing the effect of
the system, 1

Multiply by the effect to apply effect to audio data

To remove/compensate for effect: Divide by the effect to obtain JJ
the Fourier transform of the ideal image. II
J
Inverse Fourier transform to recover the new audio/ideal image. I
This process is sometimes referred to as deconvolution. Back
Close
Image Deblurring Deconvolution Example
Recall our Low Pass (Butterworth) Filter example of a few slides ago:
butterworth.m:
deconv.m and deconv2.m reuses this code and adds a deconvolution CM0268
MATLAB
stage: DSP
GRAPHICS

327
Our computed butterworth low pass filter, H is our blurring function
So to simply invert this we can divide (as opposed to multiply) by
H with the blurred image G effectively a high pass filter
Ghigh = G./H;
ghigh=real(ifft2(double(Ghigh)));
figure(5) 1

imshow(ghigh)

in this ideal example we clearly get F back and to get the image JJ
simply to inverse Fourier Transfer. II
J
In the real world we dont really know the exact blurring function I
H so things are not so easy. Back
Close
deconv.m results

CM0268
MATLAB
DSP
GRAPHICS

328

(a) Input Image (b) Blurred Low-Pass Filtered Image (c) Deconvolved Image
JJ
II
J
I
Back
Close
deconv2.m results

CM0268
MATLAB
DSP
GRAPHICS

329

(a) Input Image (b) Blurred Low-Pass Filtered Image (c) Deconvolved Image
JJ
II
J
I
Back
Close
Deconvolution is not always that simple!
Origial Image Deconvolved

CM0268
MATLAB
DSP
GRAPHICS

330

(a) Input Image (b) Deconvolved Image

Image with Noise Added Deconvolved Noisy Image

JJ
II
J
I
(c) Input Noisy Image (d) Deconvolved Noisy Image
Back
Close
High Pass Filtering
A High Pass Filter is usually defined as 1 - low pass = 1 H :
Original image High Pass Filtered

CM0268
MATLAB
DSP
GRAPHICS

331

(a) Input Image (b) High Pass Filtered Image

Image with Noise Added High Pass Filter Noisy Image

JJ
II
J
(c) Input Noisy Image (d) High Pass Filtered Noisy Image
I
Back
Close