Frequency Domain Processing (Fourier Transform)

**Objectives:
**To learn about
the Discrete Fourier Transform (DFT)

To learn about the Fast Fourier Transform (FFT)

Please
DO NOT Print this lab as it will not appear correctly

Preparations:

As usual, create a MS Word document and submit your answers and the requested figures as one file for the in-lab activities and another one for the post-lab.

Fourier Transform

In this lab, we will see some of the applications of DFT and FFT. In class we talked about the 1-D and 2-D discrete Fourier Transforms. The 1-D Discrete Fourier Transform (DFT) is defined F(u) as:

The 1-D Inverse Discrete Fourier Transform (IDFT) is defined as:

In the first equation *i = 0, 1,..., M-1* represent the number
of samples of the continuous function and in the second equation those
corresponds to their corresponding Discrete Fourier Transform
coefficients.

Any *x, *(*x _{i}*) can be computed using:

where

where the second term (

The 2-D DFT is defined as:

The 2-D IDFT is defined as:

For the 2-D case, in these equation *i = 0, 1,..., N-1* and *j
= 0, 1, ..., M-1* represent the number of samples of the continuous
function in the *x *and *y* direction respectively and in
the second equation these values correspond to the Discrete Fourier
Transform of those points.

Similar to the 1-D case, any *x* or *y *can be computed
using the initial value and the number of grids that that point is away
from the origin. Thus,

,

from these two, we will find:

,

That being said, let's solve some examples.

**Fourier Transform
in MATLAB **

We want to first conduct some experiments with Fourier Transform.
If you run the following three commands you will produce a function
displayed on Figure (1).

f = zeros(512,512);

f(245:265,235:275)=1;

imshow(f)

Figure (1)

We can find the Fourier Transform of this function:

fft2 produces real and imaginary values of a 2-dimensional
function such as f. The
real values are as small as -5.3731e-007 and as large as
183. As we mentioned in class one way to avoid the problem
dealing with both the imaginary and real values is to compute the
spectrum of the result. The spectrum is defined as:

|F(u,v)|
= (R(u,v)^{2} + I(u,v)^{2} )^{1/2}

Where R represents the
real part and I denotes the
imaginary part. We can find the real and imaginary part of the result
using the following two commands respectively:

r = real(F) and i = imag(F).

As you may have noticed it is not possible to display negative
values. Also displaying large values would not produce a good grayscale
image. Sometimes to display the result of a Fourier Transform we
use a log transformation of
the positive values to convert the large values to a more reasonable
range. Note that log will
reduce the values between:

[10-100] to [1-2],

[100-1000] to [2-3],

and so forth.

>> F2 = log(abs(F));

>> imshow(F2)

This produces:

Figure (2)

You can see the periodicity behavior on this image. The four
bright corners come from periodicity property of Fourier
Transform. Usually, we wish to see the very bright spots (main
phase) at the center. One way to do this is to multiply every
element of the mesh by (-1)^{x+y}, i.e., by (-1) to the
power of the sum of indices on x and y axes. This will shift the
phase to the center and distribute other points around the
center. For example a pixel at row 2 and column 3 would be
multiplied by (-1)^{2+3}= -1. Following is a codes
segment that does this for us.

One way to do this is to use a for loop (you can create a temp.m
file and copy and paste these lines there to run):

%
temp.m

F = fft2(f);

F = fft2(f);

F2 =
log(abs(F));

imshow(F2)

for
i = 1:512

for j = 1:512

f(i,j) = (-1)^(i+j)*f(i,j);

end

end

F =
fft2(f);

F2 =
log(abs(F));

figure,imshow(F2)

This produces:

Figure (3)

There is also a command in MATLAB that has the same effect:

Fc = fftshift(F) % to shift
the values after computing the Fourier transform

The result of the shift should be the same as that we obtained by
multiplying pixels by (-1)^{x+y}.

Lab Activity (1)

Assuming you have gone through above example and have obtained Figure
(2) and (3). We now want to conduct an experiment. By
executing the following commands, you will move the rectangle on the
original image by 40 pixels to the left.

>> f = zeros(512,512);

>> f(245:265,195:235)=1;

>> imshow(f)

Figure (4)

How do you think the centered Fourier transform of the new rectangle
will look like? You can do this by modifying the temp.m file by replacing the line
that produces the rectangle with:

f(245:265,195:235)=1;

that moves the line to the left.

Run the script.

Question1:
What difference moving to the left by 40 pixels had on your
output?

In addition to moving to the left by 40 pixels, try 3 more cases by
moving the rectangle

1) move the original line up by
40 pixels,

2) move the original line up and
left by 40 pixels on each direction

Question2:
Explain the affect of moving in each direction on the output
Fourier transform. Include your answers and the resulting 4 images in
your in-lab file.

Lab Activity (2)

In this activity, we will compute a 1-D signal and will compute the spectrum and phase for that signal. We begin by creating 20 values between 0 to 720 degree:

t =
[0: 10 : 720]*pi/180; % 20 angles between 0 to
720 in Radian (2 full period)

s = sin(t);

plot(t,s);

Figure (5)

Fsin
= fft(s);

F2 = abs(Fsin);

figure, plot(t, F2);

Figure (6)

Fc =
fftshift(F2); % to center it

figure, plot(t, Fc);

Figure (7)

Just to highlight a few things, Let's repeat this for sin(2x).
Sin(2x) is the same wave function as sin(x) but the frequency is twice
larger.

s = sin(2*t);

figure, plot(t,s);

Fsin
= fft(s);

F2 = abs(Fsin);

figure, plot(t, F2);

What kind of difference do
you observe, if any?

I now will create a signal: sin(t) +
sin(2t) + sin(4t) and will see what the above analysis will do on this
one.

s =
sin(t) + sin(2*t) + sin(4*t);

figure, plot(t,s);

Fsin
= fft(s);

F2 = abs(Fsin);

figure, plot(t, F2);

What do you learn from this
spectrum?

Question3:

Consider Figure (5), what is the frequency (number of cycles per unit
time (sec) ) of the signal? What is the amplitude (the largest
absolute value) of this signal?

Question4:

What kind of relationship do you
observe between the sin signal and its Fourier Transform?

Question5:

On Figure (6), What will the
distance between the two maximums suggest? What is the relation
between these peaks and the behavior of the signal on Figure (5).

Lab
Activity (3) - If you
didn't finish this part in the lab, do it for post-lab

One of the applications of Fourier Transform is to locate features in
an image by performing convolution. For instance occurrences of a
letter in a text. The MATLAB Image Processing Toolbox has an
example in which the letter **a **has been extracted from a sample
text. The procedure for extracting that a letter is explained
below. Before you try this please copy three image files: text.jpg , e.jpg
and n.jpg
to your directory. Please do not chnage the size of the
images. Save and use them as they are.

These two images are the images of letters "n" and "e" with the same size as they appeared in the text.jpg image file. The idea (if it works) is to find where the location of "e" or "n" is in the image shown below.

text.jpg Image

The correlation of the image e.jpg and text.jpg can be computed by:

1) Read the text.jpg image and store it in an array
called text. Read the n.jpg and e.jpg images and store them in
corresponding n and e arrays.

text = imread('text.jpg');

n = imread('n.jpg');

e = imread('e.jpg');

2) use imshow to display the text image.

3) use figure and imshow to display one of the
letter images (let's say n).

4) Rotate the letter image by 180:

rn = rot90(n, 2); % this rotates n by 90
twice (180)

5) Compute the Fourier Transform of the rotated
image and bring the size to 256x256 by padding. Note that the text
image is 256x256.

rn = double(rn); % convert to double

fn = fft2(rn,
256, 256);

6) Compute the Fourier Transform of the text image:

ftext = fft2(double(text));

7) Multiply element-by-element the resulting two
transfroms.

invboth = ifft2( ftext .* fn);

8) Take the real part of the result.

realboth = real(invboth);

9) Find the MAX value in invboth. You need this
value for Thresholding. Suppose you found that the MAX was 2700000.
Then set the threshold value a little below this. Let's say 2500000.

maxt = max(max(realboth))

Suppose you got 2700000. This is just an example. Set a
threshold value based on the maximum value, for example:

thres = 2500000;

Note: I usually take the MAX and subtract 10% from it to get the
threshold. Thus, thres = maxt - 0.1*maxt

10) Display the result with pixel values above the
threshold ONLY.

figure, imshow( realboth > thres);

You should see a blank page with bright dots on it. You may need to
play a little bit with the threshold value to fine tune the image. The
bright dots locate where that letter has appeared in the text.

Try the second letter. HAVE FUN. Submit the list of commands
and the resulting images with your file.

**Post-lab - Due Wed
Nov. 1**

Problem1. In MATLAB we have a function imrotate that rotates an image by
a given angle. For this problem, you will try different rotation
angles between 0 to 180 to find out the effect of rotation on Fourier
Transform. We will use the rectangle image (Figure 1) that we
used in the lab to complete this lab.

f = zeros(512,512);

f(245:265,235:275)=1;

imshow(f)

Perhaps the best way to complete this post-lab is modify the temp.m
script.

Produce the centered Fourier Transform of the rectangle after

1.1) 45 degree
rotation,

1.2) 90 degree
rotation

1.3) 135 degree
rotation, and

1.4) 180 degree
rotation.

The syntax for imrotation
is:

B = IMROTATE(A, ANGLE, METHOD)

Where A is the given image, ANGLE is the counter-clockwise rotation
angle in degree, and METHOD is the interpolation method. As you
may have noticed when you rotate an image by an angle, let's say 45,
the edge of the shape may become the diagonal of some pixels. A
pixels may only have one gray value representing it. Thus, you
will settle for a value that best describes the pixel. There are
several METHODS to determine such a value for a questionable pixel.
Three are used by MATLAB:

1) 'nearest' (default) nearest neighbor
interpolation

2) 'bilinear' bilinear interpolation

3) 'bicubic' bicubic interpolation

For now you can use the second method bilinear interpolation.

Important: Also, when you rotate
the 512x512 image the size may no longer be 512x512 depending on the
rotation angle. For example for 45 degree you will not produce a
512x512 image. In such a case you need to cut the center 512x512
part of the image to use in your transformation. The best way to
make sure the image is correct before you apply the Fourier transform
is to display it first. For example, display the 45 degree
rotated image and then cut the 512x512 center part.

Example (if you rotate the 512x512
image by 45):

g = imrotate(f, 45, 'bilinear');

size(g) % This gives the size of g

725 725

If you run the imshow on the rotated image, you will get a warning:

>> imshow(g)

Warning: Image is too big to fit on screen; displaying at 56% scale.

In truesize>Resize1 at 308

In truesize at 44

In imshow at 161

A warning is displayed about the dimension.

After I cut the center 512x512 pixels, and used imshow to display
the result, I got this image. Note:
you want to cut the center 512x512 image.