Description
Basics of Signals and
Systems
1 Part 1
In this assignment, 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 signalssuch 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).
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 integers). Since practical images are of finite size (such as 512×512, 1024×1024 etc.), we usually
have x[m,n] = 0 for m ∉ [0,Mx −1] or n ∉ [0,Nx −1] where the image size is Mx × Nx .
Recall that we store 1D signals in 1D arrays in Matlab. Since images are 2D signals, they
are stored in matrices in Matlab. An Mx×Nx image is stored in an Mx×Nx 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 namedReadMyImage.m and DisplayMyImage.m and the image named
Part5.bmp under the current directory of Matlab.
Issue the command
A=ReadMyImage(’Part5.bmp’);
×
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
You will see that a 532 800 matrix of type double will be created and stored in the
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×1[n]+α2×2[n] produces α1 y1[n]+α2 y2[n] for all x1[n], x2[n], α1, α2.
• x[n −n0] produces y[n −n0] for all x[n] and n0.
Similarly, a 2D DS system is called linear space invariant (LSI) if it satisfies:
• α1×1[m,n] + α2×2[m,n] produces α1 y1[m,n] + α2 y2[m,n] for all x1[m,n], x2[m,n],
α1, α2.
• x[m −m0,n −n0] produces y[m −m0,n −n0] for all x[m,n] and (m0,n0).
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
δ[n] =
(
1 if n = 0
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[n] =
X∞
k=−∞
x[k]δ[n −k] (1)
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:
y[n] =
X∞
k=−∞
x[k]h[n −k] (2)
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
δ[m,n] =
(
1 if m = 0,n = 0
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:
y[m,n] =
X∞
k=−∞
X∞
l=−∞
x[k,l]h[m −k,n −l]
= x[m,n]∗ ∗h[m,n]
=
X∞
k=−∞
X∞
l=−∞
h[k,l]x[m −k,n −l] (4)
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,Mh −1] where Mh ∈ 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,Mh −1] and n ∉ [0,Nh −1] where Mh,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 Mx × Nx ) is input to a 2D FIR DS LSI system whose impulse
response h[m,n] is of size Mh ×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 ≤ Mx −1 and 0 ≤ n ≤ Nx −1.
• h[m,n] can only be nonzero within 0 ≤ m ≤ Mh −1 and 0 ≤ n ≤ Nh −1.
Under these conditions Eq. 4 reduces to:
y[m,n] =
M
X
h−1
k=0
N
X
h−1
l=0
h[k,l]x[m −k,n −l] (5)
Based on the above equation, show that y[m,n] can only be nonzero within 0 ≤ m ≤ My −1
and 0 ≤ n ≤ Ny − 1. Determine My and Ny in terms of Mx , Nx , Mh 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 Mh×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 Mx×Nx 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 My × Ny 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 My × Ny
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
x =
2 1 1
−3 0 2
1 −1 2
and
h =
·
2 1
−1 0¸
5
6
4 Part 4 – Adding Noise and 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. In this part, first
of all we add noise to the picture and we will try to rescue thisimage from the noise without
disturbing the image itself as much as possible.
For the first step, write a simple MATLAB code to add Gaussian noise with different mean
and variances as given in the table below, and then comment about the differences you see
in the picture.
In order to add noise you can use the following command:
random(’norm’, mean , std ,size)
if you need further information about this command, look for it in MATLAB help.
Table 1: Noise parameter
σ µ
0.1 0
0.25 0
0.5 0
Since we know that the noise on the picture has high frequency content (just as many
other types of noise that we encounter in signal processing). 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:
where B is a free parameter that determinesthe bandwidth of the filter. We should select B
between 0 and 1.
Let D denote your ID number, and let D17 denote your ID number in modulo 17. That is
D ≡ D17 mod17
with 0 ≤ D17 ≤ 16. You can compute D17 using the rem command of Matlab. To learn how rem
works, type help rem in Matlab command window.
Now, let Mh = Nh = 20 + D17 and B = 0.5. 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 noise standard
deviation equal to 0.5 using this filter. Display the output image. Repeat the exercise with B
= 0.2 and B = 0.05 as well. Include all the output images to your report. Use subplot
command, and show the images on the same Matlab figure.
Comment on your results. Which value of B seems to be the most appropriate one? Include 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 h1 is equal to:
7
8
1
2
Prepare h1[m,n] and process your image with thisfilter. Let y1[m,n] denote the resulting
image. Display the image which is defined as s1[m, n] = y 2
[m, n]. Include this image to
your report. Which parts of the original image are emphasized? Include your answer to
your report.
Now let
Again, prepare h2[m, n] and process your image with this filter. Let y2[m, n] denote the
resulting image. Display the image which is defined as s2[m, n] = y 2
[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 h1[m, n]. Include your answer and
comments to your report.
6 Part 6 – A sample pattern recognition
application
A main problem of image (orsignal) 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 pattern 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.
9
To make the bright points more visible, compute the image |y [m, n]|
3
and display it. Include 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 without any confusion? Include your answer to your report.
Comment on the success of the method. Include your comments to the report.
