Digital Signal Processing Lab 2

Filtering Applications and Correlation

Nikos Deligiannis,
Boris Joukovsky ([email protected]),
and Brent De Weerdt ([email protected])

1 Filtering in the frequency domain

Exercise 1. Filtering in the Fourier domain

Consider a Gaussian filter with dimensions 64 × 64 and σ = 0.5. You can implement it in Matlab
using the command fspecial.

filter=fspecial(’gaussian’,[64,64], 5.0)

i. Filter ‘gertrude.tif’ (after the conversion to double format) in the time domain using
conv2.

ii. Convert the filter and the image to the frequency domain and execute the filtering again,
corresponding to the pointwise multiplication of the spectrum of the signal and the filter. Pay
attention that you don’t introduce time domain aliasing (see below). Finish the operation by
performing the inverse transformation of the filtered signal.

iii. Compare the time that both operations take to complete (use tic and toc). Is the result the
same as the result through filtering in the time domain? Why is that? Is it useful to take a
detour to the frequency domain?

Time domain aliasing in 1D (similar for 2D): consider a signal x(n) with N samples and a FIR
filter h(n) with L samples. We know that through filtering, we get a an output y(n) with N + L − 1
meaningful samples. If we filter in the frequency domain, then Y (k) = DF T (y(n)) should also
be N + L − 1 long in order to match the result obtained in the time domain, which requires that
X(k) = DFT(x(n)) and H(k) = DFT(h(n)) are also N + L − 1 samples long. In other words, we
have to execute a N + L − 1 point DFT transformation for both signals to avoid time
domain aliasing. From last exercise section, this can be done in practice either by specifying the
number of samples in the fft command, or by padding the signal x(n) with L − 1 zeros.

1
Exercise 2. Overlap-save method

Figure 1: The overlap-save method

The overlap-save method is often used to filter very long signals or in real-time digital filtering ap-
plications. Figure 1 shows the working principle of this method. In this exercise, you will implement
the overlap-save method and apply it to an audio signal, which is given to you in the “song.mat”
file. You can play it using the following command (smile if you hear something familiar):

load song;
fs = 2e3;
sound(song, fs);

Follow the skeleton given in the file “overlap save.m” to finish the exercise. You will notice that
the filter is already provided in the matlab file. For each of the chunks, do not forget to perform
the filtering in the frequency domain.

2
2 Applications
Sometimes when you examine input data you may wish to smooth the data in order to see a trend
in the signal. Smoothing is how we discover important patterns in our data while leaving out
things that are unimportant (i.e. noise). Filters are data processing techniques that can smooth
out high-frequency fluctuations in data or remove periodic trends of a specific frequency from data.
The goal of smoothing is to produce slow changes in value so that it’s easier to see trends in our
data.
The following difference equation describes a filter that averages time-dependent data with
respect to the current time instance and the three previous time instances of data.
1 1 1
y(n) = x(n) + x(n − 1) + x(n − 2). (1)
3 3 3
The filter function
y = filter(b, a, x)
filters a vector of data x according to the following difference equation, which describes a tapped
delay-line filter:

a1 y(n) = b1 x(n) + b2 x(n − 1) + · · · + bNb x(n − Nb + 1)

− a2 y(n − 1) − · · · − aNa y(n − Na + 1), (2)

where a = [a1 a2 . . . aNa ] and b = [b1 b2 . . . bNb ] are the vectors of coefficients of the filter, Na is the
feedback filter order and Nb is the feedforward filter order. Thus, the output is a linear combination
of the current and previous elements of x and y.

Exercise 3. Smoothing a financial signal

The use of an algorithm to remove noise from a financial data set allows important patterns to stand
out. In our example we are given a dataset containing weekly closing stock prices from 6/1/2012
until 29/4/2016.

i. Load the dataset in matlab using the following command:

load(’stock data1.mat’)

ii. The dataset contains stock prices stored in the variable x and the corresponding dates stored
in the variable t. Plot your data to confirm their noisy nature.

iii. Create a moving average smoothing filter using the following commands:
a = 1;
avnum = 10;
b = 1/avnum * ones(1, avnum);

iv. Apply the filter to your data

y = filter(b,a,x);
and plot the filtered data. Experiment with different values of avnum.

3
Exercise 4. ECG Signals

Electrocardiography (ECG or EKG) is the process of recording the electrical activity of the heart
over a period of time using electrodes placed on the skin. The QRS complex is a name for the
combination of three of the graphical deflections seen on a typical electrocardiogram as you can
see in Fig. 2. It is usually the central and most visually obvious part of the tracing. In adults, it
normally lasts 0.06 − 0.10 secs; in children and during physical activity, it may be shorter. The Q,
R, and S waves occur in rapid succession, do not all appear in all leads, and reflect a single event,
and thus are usually considered together.
Next, we will experiment with ECG signals downloaded from the PhysioNet database, which was
created by a team of investigators from Children’s Hospital Boston (CHB) and the Massachusetts
Institute of Technology (MIT) (https://www.physionet.org).

1. We will load an ECG signal recorded at 360 samples per second. Our signal is stored in
‘103.dat’. Load the signal in matlab using the following commands:
filename = ‘103.dat’;
fid=fopen(filename,‘r’);
fs = 360;
time=10;
f=fread(fid,2*fs*time,‘ubit12’);

2. Subsample the signal so it can be better visualized and plot it:

x ecg=f(1:2:length(f));
figure; plot(x ecg); title(‘input ECG signal’);

3. Smooth the signal using a moving average filter. Experiment with different values of averaging
numbers and produce a signal that removes the noise, while it keeps the QRS information of
the signal.

4. Try to locate the “R waves ” using the matlab command findpeaks. The command takes as
input the addressed signal and the parameters ‘MinPeakHeight’ and ‘MinPeakDistance’.
[~, locs Rwave] = findpeaks(ECG data,‘MinPeakHeight’, ...
min ampl, ‘MinPeakDistance’,min dist);

5. Using the R values try to locate the Q and S positions by finding peaks in the neighborhood
of an R value.

Figure 2: Electrocardiography

4
Exercise 5. Removing Salt and Pepper noise (binary noise)

When an image is suffering from salt and pepper noise the image is contaminated with black and
white pixels (see Fig. 3). This type of noise can be removed in a simple way by using median
filtering. Median filtering is a non-linear operation where the target pixel is the median of the
source pixel and its 8 neighbouring pixels. Figure 4 clarifies this.

Figure 3: Salt and Pepper noise.

z=median([a b c d a f a e a])
Original image Filtered image

a b c a b ... ? ? ? ? ? ...
d a f a a ? z ? ? ?
a e a g a ? ? ? ? ?
a f a e d ? ? ? ? ?
. .
. .
. .

Figure 4: Median filtering.

Load the greyscale image ‘lenagray.png’. Convert to double format. Add salt and pepper
noise with the imnoise command. Now use the medianfilter function (given) to apply me-
dian filtering. Look at the result. Also look at the implementation of the medianfilter function.

Exercise 6. Unsharp masking

Unsharp masking is an operation that is applied to strengthen the edges in an image. The operation
takes the original image and adds a signal which is proportional to the difference between the original
image and a lowpass filtered version of this image. This is equivalent to taking the original image
and adding a highpass filtered version of this image. The unsharp masking operation is described
as:
u(m, n) = v(m, n) + λg(m, n) (3)

5
where v is the original image, u the image after unsharp masking and g the highpass filtered image;
λ is always larger than 0.

i. Load the image ‘aerial.tif’ and convert it to double.

ii. Execute the process of unsharp masking by yourself. Use the sharp filter as a high-pass filter:
sharp = 1/9 *[-1 -1 -1; -1 8 -1; -1 -1 -1]
Apply it to the image using conv2 command:
conv2(image, sharp,‘same’)
and display the obtained highpass filtered image. Use equation (3) to complete the unsharp
masking operation. Compare the result with the original image.

iii. Use the matlab command fspecial to create an unsharp masking filter:
unsharp=fspecial(‘unsharp’);
Apply it to the image using conv2 command and display the result. Compare it with the
unsharp masking operation you implemented in step 2.

6
3 Correlation: Template Matching
Exercise 7. Template matching using correlation

The goal of this exercise is to compare two images based on a correlation coefficient. Consider
the reference image ‘lena.png’ as template image and ‘lenalips.png’ as target image. The
target image is “placed” over the template image and the correlation coefficient for each pixel in the
template image is computed to construct the correlation map. After sliding through all the pixels
in the template image, the maximum coefficient is obtained from the map. The pixel position with
maximum value is the starting point of the target image.
To calculate the correlation coefficient you can use the Matlab command corr2:
r = corr2(A,B)
which returns the correlation coefficient r between matrix A and matrix B; r is a scalar double.

50

100

5
150

10
200
15
250
20
300
25
350
30
400
35
10 20 30 40 50 60
450

500
50 100 150 200 250 300 350 400 450 500

Figure 5: Template image (left), target image (middle) and template matching (right).

Note: A colored RGB image with dimensions m, n is represented in matlab by a 3-dimensional

tensor M of size m × n × 3. You need to apply the correlation operation only to one “slice” M(:,
:, 1) of this tensor.