B ILKENT U NIVERSITY E LECTRICAL AND E LECTRONICS E NGINEERING D EPARTMENT

EEE - 321: Signals and Systems

Lab Assignment 4

Deadlines:
Section 1: 23.11.2020, 23:55
Section 2: 26.11.2020, 23:55

Please carefully study this assigment before coming to the laboratory. You may begin
working on it or even complete it if you wish, but you do not have to. There will be short
quizzes both at the beginning and end of the lab session; these may contain conceptual,
analytical, and Matlab-based questions. Within one week, complete the assignment in the
form of a report and turn it in to the assistant. Some of the exercises will be performed
by hand and others by using Matlab. What you should include in your report is indicated
within the exercises.
Note: Along with this pdf file, you will find a .rar file containing two Matlab functions and
several bmp pictures. Unzip the archive and place all files contained in it under the current
directory of Matlab.

1 Part 1
In this lab, you are going to do a bit of image processing. You will perform some simple
operations on gray scale (white and black) digital images and see the results.

In this course up to now, we have mainly dealt with one-dimensional signals such as x(t )
or x[n]. These signals are said to be one-dimensional because they are functions of one
independent variable (t or n).

1
 Images, on the other hand are signals of two independent variables. Therefore they are
two-dimensional signals.

A gray scale digital image can be represented with a 2D discrete function x[m, n] such
that x[m, n] denotes the light intensity of the (m, n)th pixel of the image (m and n are inte-
gers). Since practical images are of finite size (such as 512×512, 1024×1024 etc.), we usually
have x[m, n] = 0 for m ∉ [0, M x − 1] or n ∉ [0, N x − 1] where the image size is M x × N x .

Recall that we store 1D signals in 1D arrays in Matlab. Since images are 2D signals, they
are stored in matrices in Matlab. An M x ×N x image is stored in an M x ×N x matrix in Matlab.

In this first part, you will just practice reading images in Matlab and displaying them. You
will not write any code for this part but use the m-files provided for this assignment.

Place the m-files named ReadMyImage.m and DisplayMyImage.m and the image named
Part5.bmp under the current directory of Matlab.

Issue the command

A=ReadMyImage(’Part5.bmp’);

You will see that a 512 × 512 matrix of type double will be created and stored in the
workspace. This matrix contains the image. The original image is a color image, but it
is converted to a gray scale image by the ReadMyImage function.

Next, issue the command

DisplayMyImage(A);

A new figure window will open and the image will be displayed.

During the rest of this assignment, you will use these commands to read and display
various images.

For this part, just make sure that you run the functions above and see the image without
any problem. You do not need to provide anything for the report.

2 Part 2
Recall that 1D discrete time (DT) systems map a 1D input signal x[n] to another 1D signal
y[n] (that we name the output signal). (We say discrete time since the independent variable
usually denotes time.)

2
 Similarly, a 2D discrete space (DS) system maps a 2D input signal x[m, n] to a 2D output
signal y[m, n]. (Now we say discrete space since the independent variables usually denote
the space coordinates.) Recall that a 1D DT system is called linear time invariant (LTI) if it
satisfies the following two properties:

• α1 x 1 [n] + α2 x 2 [n] produces α1 y 1 [n] + α2 y 2 [n] for all x 1 [n], x 2 [n], α1 , α2 .

• x[n − n 0 ] produces y[n − n 0 ] for all x[n] and n 0 .

Similarly, a 2D DS system is called linear space invariant (LSI) if it satisfies:

• α1 x 1 [m, n] + α2 x 2 [m, n] produces α1 y 1 [m, n] + α2 y 2 [m, n] for all x 1 [m, n], x 2 [m, n],
α1 , α2 .

• x[m − m 0 , n − n 0 ] produces y[m − m 0 , n − n 0 ] for all x[m, n] and (m 0 , n 0 ).

Now let us develop the input output representation of 2D DS LSI systems by forming an
analogy with 1D DT LTI systems.

Recall that a 1D discrete impulse signal is defined as

(
1 if n = 0
δ[n] =
0 otherwise

While developing the input-output relation for a 1D DT LTI system, we first wrote the
input signal x[n] as a superposition of shifted impulse signals as:
∞
X
x[n] = x[k]δ[n − k] (1)
k=−∞

where we interpreted x[k] as the coefficient of the impulse shifted by k units (that is, x[k]
is the coefficient of δ[n − k]). Then, we named the response that the system gives to δ[n] as
h[n] (impulse response). Using the time invariance property of the system, we recognized
that the response of the system to δ[n − k] should be h[n − k]. Then, using the linearity
property of the system, we wrote the input-output relation of the 1D DT LTI system as:
∞
X
y[n] = x[k]h[n − k] (2)
k=−∞

We named the above operation as the convolution of x[n] and h[n] and used the short hand
notation
y[n] = x[n] ∗ h[n] (3)
to denote it.

3
 Now, define the two dimensional discrete impulse signal as
(
1 if m = 0, n = 0
δ[m, n] =
0 otherwise

and going through the same steps, derive the input-output relation of a 2D DS LSI system,
or derive the formula for 2D convolution of x[m, n] and h[m, n] that we denote as:
x[m, n]∗∗h[m, n]. In other words, derive the 2D version of Eq. 2. Show—explain your steps
as is done above for the 1D case. Include your work to your report. You should obtain the
following result:
∞
X ∞
X
y[m, n] = x[k, l ]h[m − k, n − l ]
k=−∞ l =−∞
= x[m, n] ∗ ∗h[m, n]
∞
X ∞
X
= h[k, l ]x[m − k, n − l ] (4)
k=−∞ l =−∞

3 Part 3
Recall that a 1D DT LTI system is called FIR if the impulse response h[n] contains a finite
number of nonzero values such that h[n] = 0 for n ∉ [0, M h −1] where M h ∈ Z + . Similarly, a
2D DS LSI system is called FIR if the impulse response h[m, n] contains a finite number of
nonzero values such that h[m, n] = 0 for m ∉ [0, M h − 1] and n ∉ [0, Nh − 1] where M h , Nh ∈
Z +.

In this part, you will write a Matlab function that computes the output when a finite sized
input image x[m, n] (of size M x × N x ) is input to a 2D FIR DS LSI system whose impulse
response h[m, n] is of size M h × Nh . From another perspective, your code will compute the
2D convolution of two finite-size 2D signals x[m, n] and h[m, n]. We assume that

• x[m, n] can only be nonzero within 0 ≤ m ≤ M x − 1 and 0 ≤ n ≤ N x − 1.

• h[m, n] can only be nonzero within 0 ≤ m ≤ M h − 1 and 0 ≤ n ≤ Nh − 1.

Under these conditions Eq. 4 reduces to:

MX
h −1 NX
h −1
y[m, n] = h[k, l ]x[m − k, n − l ] (5)
k=0 l =0

Based on the above equation, show that y[m, n] can only be nonzero within 0 ≤ m ≤ M y −1
and 0 ≤ n ≤ N y − 1. Determine M y and N y in terms of M x , N x , M h and Nh . Include your
work to your report.

Next, write a Matlab function of the following form

4
 function [y]=DSLSI2D(h,x) where

• h of size M h ×Nh denotes the impulse response of the system, such that h(1,1)=h[0, 0],
h(1,2)=h[0, 1], ..., h(k,l)=h[k − 1, l − 1].

• x of size M x ×N x denotes the input signal x[m, n], such that x(1,1)=x[0, 0], x(1,2)=x[0, 1],
..., x(k,l)=x[k − 1, l − 1].

• y of size M y × N y denotes the output signal, such that y(1,1)=y[0, 0], y(1,2)=y[0, 1], ...,
x(k,l)=y[k − 1, l − 1].

Basically, you will implement Eq. 5. You should first create a zero matrix of size M y × N y
for y. Then, you can use a nested for loop as in the following code:
for k=0:Mh-1
for l=0:Nh-1
y(k+1:k+Mx,l+1:l+Nx)=y(k+1:k+Mx,l+1:l+Nx)+h(k+1,l+1)*x;
end
end

Note: Do not use any built in command of Matlab. Directly implement Eq. 5.

Note: Before writing your code, carry out the mini exercise below. You do not need to
provide anything for the report, but make sure that you fully understand the interpretation
of the following commands.

• Type A=zeros(3,3);. Then type A(1:2,1)=[1;2];. Examine A.

• Now type A(2:3,2:3)=[5 8;6 9];. Examine A.

• Now type A(3:3,1:1)=[3];. Examine A.

• Now type A(1:1,2:3)=[4 7];. Examine A.

• Let B=A(1:2,2:3);. Examine B.

• Let B=A(3:3,1:2);. Examine B.

Check your function: If you take

 
8 1 6
x = 3 5 7 
4 9 2
and · ¸
1 3
h=
4 2

5
you should get  
8 25 9 18
35 34 48 33
 
y=
16 47 67 20

16 44 26 4
Include your code to your report.

4 Part 4 - Image Denoising

Now it is time to see some practical image processing examples.

Read the picture named Part4.bmp in a matrix named x and display it. You will see a
picture that is highly corrupted by noise. In this part, we will try to rescue this image from
the noise without disturbing the image itself as much as possible.

Suppose we know that the noise on the picture has high frequency content (just as many
other types of noise that we encounter in electronics engineering). We also know that typi-
cal daily life images taken by a typical camera have low-frequency content. Therefore, why
should not we apply a low pass filter to the noisy image and try to eliminate the noise? A
typically used 2D FIR low pass filter has the following impulse response:
(
0 if m < 0 or n < 0 or m > M h − 1 or n > Nh − 1
h[m, n] = n ³
M h −1
´o n ³
Nh −1
´o
sinc B m − 2 sinc B n − 2 otherwise

where B is a free parameter that determines the bandwidth of the filter. We should select B
between 0 and 1.

Let D denote your ID number, and let D 17 denote your ID number in modulo 17. That is

D ≡ D 17 mod17

with 0 ≤ D 17 ≤ 17. You can compute D 7 using the rem command of Matlab. To learn how
rem works, type help rem in Matlab command window.

Now, let M h = Nh = 20 + D 17 and B = 0.8. Prepare the h matrix that represents h[m, n].
If you wish, you can use nested for loops over m and n, this time it is allowed. Recall that
Matlab already has a built in function named sinc to compute the values of sinc(.) function.
Include your code for preparing h to your report.

Using the code that you developed in Part 3, process the noisy image with this filter. Dis-
play the output image. Repeat the exercise with B = 0.5 and B = 0.2 as well. Include all the
output images to your report. Use subplot command, and show the images on the same
matlab figure.

6
 Comment on your results. Which value of B seems to be the most appropriate one? In-
clude your comments to your report.

As you see, it is possible to greatly clarify the information in a corrupted signal using only
the simple concepts of convolution, LSI systems, etc. In this part, we tried to remove the
noise from a corrupted image. Such problems are called image denoising problems. Many
algorithms have been developed by many researchers that try to clear the images from the
corrupting noise while giving minimum damage to the original image.

5 Part 5 - Edge Detection

Another widely studied interesting problem of image processing is the detection of edges
in an image, which is called the edge detection problem. In a typical image, edges are the
set of points in the close vicinity of which the pixel values change abruptly. Since changes
are sudden and great, edges are inherently associated with high frequencies. Therefore, we
can make use of high pass filters to detect edges.

Read the picture named Part5.bmp in a matrix named x and display it.

Consider a 2D FIR DS LSI system whose impulse response h 1 [m, n] is equal to:

1

 if m = 0, n = 0
h 1 [m, n] = −1 if m = 0, n = 1

0

otherwise

Prepare h 1 [m, n] and process your image with this filter. Let y 1 [m, n] denote the resulting
image. Display the image which is defined as s 1 [m, n] = y 12 [m, n]. Include this image to
your report. Which parts of the original image are emphasized? Include your answer to
your report.

Now let 
1

 if m = 0, n = 0
h 2 [m, n] = −1 if m = 1, n = 0 .

0

otherwise
As you see, h 2 [m, n] = h 1 [n, m]. Again, prepare h 2 [m, n] and process your image with
this filter. Let y 2 [m, n] denote the resulting image. Display the image which is defined as
s 2 [m, n] = y 22 [m, n]. Include this image to your report. Which parts of the original image are
emphasized now? Comment on the difference with the image you obtained with h 1 [m, n].
Include your answer and comments to your report.

7
 Finally, process the original image with the filter h 3 [m, n] = 0.5h 1 [m, n] + 0.5h 2 [m, n].
Let y 3 [m, n] denote the resulting image. Display the image which is defined as s 3 [m, n] =
y 32 [m, n]. Include this image to your report. Which parts of the original image are empha-
sized now? Comment on the difference with the images you obtained with h 1 [m, n] and
h 2 [m, n]. Include your answer and comments to your report.

6 Part 6 - A sample pattern recognition application

A main problem of image (or signal) processing is the detection of certain objects within a
picture.

Read and display the image named Part6x.bmp. You will see our national soccer team.
Our purpose in this part is to detect the faces in the image. Such problems are named pat-
tern recognition problems. Again, many algorithms have been developed to solve such
problems. Here, we will use one of the basic methods which is called matching filter
method. Read the image Part6h.bmp into Matlab and think of it as the impulse response
of a 2D DS LSI FIR system. Display the image to see how the impulse response looks like.
You will see an inverted face. That is, to detect the faces within the image, we are using a 2D
DS LSI system whose impulse response is an inverted face! Note that by inverted, we mean
that the face is rotated by 180 degrees. (When the rotation is 180 degrees, its direction -
clockwise or counterclockwise - is unimportant.)

Now, pass the input image through this system and find the output image. Display the
absolute value of the output image, that is, display |y[m, n]|. Include this image to your
report. Search for the points that look bright. Where do they occur? Do you always see a
face at the location that you think there is a bright point, or do you sometimes find bright
points where there is no face? Include your answers to your report.

To make the bright points more visible, compute the image |y[m, n]|3 and display it. In-
clude this image to your report. Also compute the image |y[m, n]|5 and display it. Include
this image to your report. Taking which power is sufficient for you to detect the faces with-
out any confusion? Include your answer to your report.

Comment on the success of the method. Include your comments to the report.