Frequency Domain Filtering

Matthew Thurley Industrial Image Analysis – E0005E based on material by Johan Carlson

1

——————————————————–

1-1

This Lecture The concept of frequency – Fourier transforms Repetition of 1D convolution, transforms and filtering Generalization to the 2D world 2D convolution Frequency domain filtering Examples of frequency selective filters

2

The concept of frequency Let’s start off in the one-dimensional world... The period of cos(t) is 2π seconds, i.e. The signal repeats itself after 2π seconds. The period of cos(2π · t) is one second. The period is usually denoted T . Frequency is measured in Hz (Hertz) and is the number of periods (or cycles) per second, denoted f . That is, f=

1 . T

The function cos(2π · f t) has a frequency of f and a period of 1/f Angular frequency is measured in rad/s, denoted ω, so 2π . ω= T 3

The concept of frequency (cont’d...) Cosine, f = 200 Hz

Cosine, f = 100 Hz 5 x (t)

0 −5 −10

2

1

x (t)

5

−5 −10

0 10 20 Time t (ms) Cosine, f = 10 Hz

0 10 20 Time t (ms) Cosine, f = 0 Hz

5 x (t)

0

3

3

x (t)

5

−5 −10

0

0 10 Time t (ms)

20

0 −5 −10

0 10 Time t (ms)

20

4

The concept of frequency (cont’d...) Some comments Frequency is measured in hertz, cycles or periods per second Higher frequency means more periods per second. Higher frequency means shorter period. Zero frequency implies infinite period, which means a constant signal. Fourier transform of a function of time has units, cycles per second Fourier transform of a function of pixels (like an image) will have units, cycles per pixel Rapid changes in a signal correspond to high frequencies, slow changes are represented by low frequencies Rapid changes in intensity in an image are high frequencies.

5

Background - Complex Numbers Recall that a complex number contains a real and an imaginary component. C = R + jI And a complex number can be represented in polar coordinates (r, θ), (magnitude, angle), and using these polar coordinates C = |C|(cosθ + jsinθ) Using Eulers formula ejθ = cosθ + jsinθ And if we use ω = 2π/T angular frequency as our angle j 2π T

e

2π 2π = cos + jsin T T

Eulers formula is a sum of cosins and sines that Fourier used to express any periodic function.

6

Sum of cosines and sines Using a sum of cosines and sines to express a periodic function

7

The 1D Fourier transforms The Fourier theory states that any periodic function can be expressed as a sum of cosines and sines of different frequencies with each multiplied by a different coefficient. Which means it can be expressed in terms of complex exponents (eulers formula) with varying frequency. A time signal or time series can be representedthis way also. The continuous-time Fourier transform Z ∞ X(ω) = x(t)e−jωt dt, −∞

where ω = 2πf , where f (or µ in the textbook) denotes frequency in Hz. The corresponding inverse transform is 1 x(t) = 2π

Z

∞

X(ω)ejωt dω.

−∞

8

The 1D Fourier transforms (cont’d...) Example Consider the following two continuous-time signals: x1 (t) = cos(2πt) + 0.3 cos(16πt) 2

x2 (t) = e−8(t−3) cos(8πt) What are the frequency components in x1 (t)?

9

The 1D Fourier transforms (cont’d...) ...and plotting the signals and the transforms we get 0.6 0.4 0.2 0 −20

|X1(f)|

1

x (t)

1 0 −1 0

2

4 time, t (s)

6

−10 0 10 frequency, f (Hz)

20

−3

1 0.5 0 −0.5 0

2 2

|X (f)|

x2(t)

x 10

2

4 time, t (s)

6

1

0 −10

−5 0 5 frequency, f (Hz)

10

10

The 1D Fourier transforms (cont’d...) There is also a discrete Fourier transform, working on sequences of numbers, x[n] where n = 0, 2, . . . , N − 1. This is defined as X[k] =

N −1 X

x[n]e−j(2πk/N )n ,

n=0

where we get the discrete frequencies k = 0, 1, . . . , N − 1. The corresponding inverse transform is defined as N −1 1 X x[n] = X[k]ej(2πk/N )n N k=1

These can be computed in MATLAB using the fft and ifft commands, respectively. 11

The 1D Fourier transforms (cont’d...) Example - same as before 0.6 0.4 0.2 0

|X [k]|

0

1

x1[n]

1

−1 0

20

40

60

80

7000

n 1 0.5 0 −0.5 200

9000

−3

2 2

|X ([k]|

2

x (t)

x 10

8000 k

250

300 time, t (s)

350

1 0 6500 7000 7500 8000 8500 9000 9500 k

Note that the bottom right spectrum is also discrete, but plotted here as a solid line just for clarity. 12

Background - Correlation and Convolution DIP Textbook 3.4.2 Spatial Correlation and Convolution Correlation is the processing of filtering a mask over an image, exactly as explained in linear spatial filtering or neighborhood operations in lecture 2. Correlation operator ⋄, of a filter w(x, y) of size m × n with an image f (x, y) Pa Pb w(x, y) ⋄ f (x, y) = s=−a t=−b w(s, t)f (x + s, y + s)

Convolution is the same except the mask is reflected about both axes (exactly as the transpose in binary morphological image processing. Also described as the mask being rotated by 180 degrees) Convolution operator ∗, of a filter w(x, y) of size m × n with an image f (x, y) Pa Pb w(x, y) ∗ f (x, y) = s=−a t=−b w(s, t)f (x − s, y − s) 13

Background - Correlation and Convolution

14

Background - Correlation and Convolution

15

Filtering in the time-domain If we have a system which is linear and time-invariant (LTI) x(t)

y(t) = T fx(t)g

T fx(t)g

the output can be calculated using the convolution integral y(t) = x(t) ∗ h(t) =

Z

∞

x(τ )h(t − τ )dτ

−∞

or equivalently in the discrete-time case, by the convolution sum y[n] = x[n] ∗ h[n] =

∞ X

x[m]h[m − n].

m=−∞

The convolution sum can be calculated in MATLAB using the conv command.

16

Filtering in the frequency domain We also know that for LTI systems, the Fourier transform is very powerful, since the convolution in the time domain can be replaced with a multiplication in the frequency domain, i.e. F

y(t) = x(t) ∗ h(t) ←→ Y (ω) = X(ω) · H(ω) And correspondingly in the discrete domain (in principle) F

y[n] = x[n] ∗ h[n] ←→ Y [k] = X[k]H[k] Explore this in the development phase this week. There are some details we need to look into as well...

17

Length of filtered signals When we filter a signal x[n], n = 0, 1, . . . , N − 1 with a filter h[n], n = 0, 1, . . . , M − 1, the output signal y[n] will go from n = 0, 1, . . . , (N − 1) + (M − 1) − 1. In MATLAB notation: length(y)=length(x)+length(h)-1. The output signal is longer than the input! So, if we are doing the filtering in the frequency domain, using the DFT (fft in MATLAB), we need to make the frequency vectors longer as well. This is accomplished by zero-padding the transform. The filtering is done by typing: y = real(ifft(fft(x,M+N-1).*fft(h,M+N-1)))

18

Frequency selective filtering So, by the frequency domain interpretation of filtering (convolution) we see that we can easily design frequency selective filters. Recall the previous signal example 1 |X1(f)|

x1(t)

1 0

0.5

−1 0

2

4 time, t (s)

6

0 −50

0 frequency, f (Hz)

50

We realize now that we could remove the high-frequency component by multiplying the frequency-domain with a window that cuts away the undesired frequencies, i.e. a filter which spectrum is zero for frequencies above some threshold.

19

Frequency selective filtering ...the result 1 y (t)

0

1

1

x (t)

1

−1 0

0 −1

2

4 time, t (s)

6

0

2

4 time, t (s)

6

This is an example of a low-pass filter, since low frequencies pass through undisturbed while high frequencies are cancelled out.

20

Now, let’s go 2D! In two dimensions, working with discrete signals (digital images), the convolution theorem is generalized as follows: y[m, n] =

−1 M −1 N X X

x[k, l]h[m − k, n − l],

k=0 l=0

where x[m, n] is the original image and h[m, n] is the two-dimensional filter, sometimes called filter mask. In MATLAB this is done by conv2.

21

The 2D discrete Fourier transform The extension of the Fourier transform theory to the two-dimensional case is straightforward. Since we’re working with digital images, let’s focus only on the discrete transform. The 2D discrete Fourier transform is defined as: X[u, v] =

M −1 N −1 X X

x[m, n]e−j2π(um/M +vn/N )

m=0 n=0

And the corresponding inverse transform M −1 N −1 1 X X x[m, n] = X[u, v]ej2π(um/M +vn/N ) M N u=0 v=0

These are available in MATLAB as fft2 and ifft2. 22

2D filtering in the frequency domain Details... Just as in the 1D case, the output image of the filter will be larger, due to the size of the filter mask. Using conv2, i.e. filtering in the spatial domain, the results can be cut automatically using the ’same’ option to conv2 (see MATLAB help for details). Doing the filtering in the frequency domain, the transform has to be zero-padded to the proper size. After the filtering you can then cut the image to the original size. As in the 1D case, the 2D fourier transform and its inverse are infinitely periodic (in both dimensions), ie. they wrap around. Consider this when you observe effects at the edge of your image when converting back to spatial coordinates with the inverse fourier transform. You will practice these tasks in this weeks development phase. 23

2D filtering in the frequency domain As the 2D discrete fourier transform (DFT) is complex, it can be expressed in polar coordinates with a magnitude, and an anglular frequency (also known as the phase). The components of the spectrum determine the amplitudes of the sinusoids that combine to form the resulting image. At any given frequency in the DFT image, a large amplitude implies a greater prominance of a sinusoid of that frequency in the image. The phase is a measure of the displacement of the various sinusoids with respect to their origin. While the magnitude of the 2D DFT gives the amplitudes of the sinusoids and determines the intensities in the image, the phase specifies the angles (displacements) of the sinusoids and where objects are in the image. Investigate this in the development phase this week.

24

The concept of frequency, revisited When going to 2D in general, and digital images in particular, we need to think about what we mean with frequencies. For images, the transformed quantity X[u, v] is a function of two spatial frequency variables. We no longer talk about Hertz or radians/sec, but simply spatial frequencies. Low frequencies still represent slow variations in the image, and high frequencies represent rapid changes. A smooth image, e.g. a photograph often contains large areas of slowly varying intensity. These images have a dominant part of their energy concentrated to a few low-frequency components in the Fourier-domain. Edges in images imply rapid changes in the intensity, and will therefore be represented by the high-frequency components of the image’s spectrum. 25

2D filtering Example Consider the following filter mask: h[m, n] =



1/25, 0,

m, n = 0, 1, . . . , 4 otherwise

,

i.e. a 5 × 5 moving average filter (taking the average of 25 pixels). Let H[u, v] be the 2D discrete Fourier transform of the filter mask. Also, consider an inverse filter, G[u, v] = 1/(H[u, v] + µ). (Where µ is a small value to prevent division by zero)

26

2D filtering Example (cont’d...) Looking at the 1D frequency response of the uniform averaging filter we get the following

27

2D filtering Example (cont’d...) Looking at the frequency responses of the two filters we get the following |H(u,v)|

|G(u,v)| 1

−0.5

−0.5

4

0.8 3 0

v

v

0.6 0

2

0.4

1

0.2 0.5 −0.5

0 u

0.5

0.5 −0.5

0 u

0.5

Consider, which frequencies respond strongly in each filter? Which frequencies are ”passed” (not removed) by each filter? 28

2D filtering Example (cont’d...) Some comments: The first filter seems to be a low-pass filter, since low frequencies are allowed to pass through, but higher frequencies are damped. The second filter is the inverse of the low-pass filter, so it should be a high-pass filter. The low-pass filter should cause edges and other rapid changes in the image to be blurred, while smooth surfaces should be left more or less unaffected. The second filter should keep edges but smooth surfaces should be attenuated (reduced in intensity / magnitude). Let’s look at the result...

29

2D filtering Example (cont’d...) original

grayscale

low−pass filtered

high−pass filtered

30

2D filtering Some general comments Being able to analyze the spectral content of an image provides tools to construct filters. For example, if a distortion of an image can be modeled as a linear filter, then the inverse filter should be able to compensate for the distortion. One example could be motion blurring, caused by the motion of the object or the camera. In the presence of noise, inverse filtering becomes difficult. This is because the filter may contain zeros in the spectrum. If there is noise at those frequencies, the inverse filter will blow up the noise, drowning all other interesting information. Hence the µ (small constant) in the inverse example. You will experiment with this in this weeks development phase.

31

Summary Repetition of convolution and Fourier transforms in the 1D case. Generalization to the 2D world. Examples of filtering. Chapter 4 covers frequency domain filtering in detail Chapter 5 covers image restoration and topics involving removing blur and periodic noise from an image.

32

Tips for Development Phase 6 Be aware of what the range of intensities (values) are in your images. Be careful with image rescaling as it can hide the range of intensities from you and makes small magnitures look as significant as large magnitudes. You could plot them in 3D with the mesh command. Extra reading 4.6.3 Periodicity 4.6.5 Fourier Spectrum and Phase Angle 4.7 The Basics of Filtering in the Frequency Domain 4.7.2 Frequency domain filtering fundamentals 4.7.3 Summary of steps for filtering in the frequency domain 4.7.4 Correspondance between filtering in the spatial and frequency domain

33

Fourier Descriptors for Dev Phase 8 Consider we have the perimeter of a simple convex object and we want to characterise the shape of the object. Consider a centre point within that object, and then consider that we represent the perimeter of the object in polar coordinates (r, θ) If we plot a graph of r,theta we will have a periodic function We could apply a 1D discrete fourier transform to this plot to identify the major frequencies that characterise the shape of the object. More on this in Development Phase 8, and in DIP 11.2.3

34

Summing Up Consider the following three questions; What do i need to work on? What have I learnt today? What was the main point left unanswered? Write your answers in the journal with the lecture number at the top of the page.

35