ee453 lab 3

EE 453: Lab 3
Applications of Windowing in Signal Processing
April 18, 2021
Overview
The purpose of this final lab is to see applications of windowing in signal processing. This lab will not feature Simulink, but instead focus on applications
using standard MATLAB scripting. The first application will be to antenna
array signal processing. This will be useful to anybody interested in radar or
ultrasound array image processing. The final application will be an application
of windowing in filtered back projection, the method used to invert the Radon
Transform, used routinely in medical image processing. These applications will
show some applications beyond filter design of both windowing and digital signal
processing.
1
1 System Setup
Before beginning the lab, you must make sure the system is setup as follows:
• Connect your headphones to the Headphone jack on the front of the computer
• Log on to the computer using your PSU access account. You must save
all files on your PSU W: drive. Please do not store any files
on the computer itself. It takes up hardrive space, and will be
deleted periodically. While not necessary to save the files, they will be
good to have as a reference for the future labs this semester.
• Download the .zip file from Canvas containing files needed to complete
the lab.
2
2 Antenna Array Processing
Antenna arrays are used everywhere in communication, defense, and imaging applications. One of the tremendous benefits of using arrays is that the
individual performance of any one sensor/antenna/element becomes less important in determining performance of the entire system. An array system could
be made up of many inexpensive, easy to obtain elements. The benefit of using
more specialized elements decreases with the amount of elements, as well as the
signal processing that conditions the signal.
3
The figure shows an incoming wave that is received at each sensor. The incoming signal is a plane wave, and so has constant phase along its phase front.
This means that when it impinges on the array, it must be projected down onto
the array geometry (a line in this case). Often in antenna array systems, dipoles
are used for the array elements. Anyone who has taken an antenna course will
know that the pattern of a dipole is approximately omnidirectional in the x − y
plane. This leads to the approximation of a Uniform Linear Array (ULA). Under this approximation, the elements are arranged on the same line, and each
element has an omnidirectional response. The omni-directional response means
that it doesn’t matter what the angle of arrival (AoA) of the incoming signal
is. The individual elements will give the same voltage (signal) no matter what
direction the wave comes from.
One of the most important theorems in electromagnetics is the Reciprocity
Theorem. This theorem says that if a plane wave at angle θ impinges on an
antenna, denote the signal measured at the input of the antenna as V . If this
same signal V is sent to the antenna, a wave at angle θ will be transmitted.
Equivalently, it doesn’t matter whether you are receiving or transmitting, the
signal will be the same. Therefore, whether we are trying to detect what direction a wave is coming from (like in cell phones) or are trying to send a wave in
a certain direction (beamforming), the signal processing is the same.
4
The signal at each antenna of a ULA can be written as vm = e
−jkdm sin(θ)
,
where k =
2π
λ
, where λ is the wavelength of the incoming plane wave. dm is
the distance from the first element of each element. This expression indicates
that an incoming plane wave has a linear phase progression as you go from
element to element. The final signal received is simply the superposition of all
the responses at each element,
V (θ) =
N
X−1
m=0
ame
−jkdmcos(θ)
am is an optional term that adds a different weight to each element. Notice the
similarity between the expression for total voltage, and the expression for the
Discrete Time Fourier Transform X(e
jω) = P∞
n=−∞ ane
−jωn. This similarity is
how we will use digital signal processing, specifically windowing techniques, to
control the performance of the antenna array.
By analogy, the array response can be thought of as the DFT of the weights
that are applied to each element. Different windows can be used as the weights,
with the array response defined by their Fourier Transform. There are two main
concerns when evaluating the performance of an array. The first is the width
of the main beam. This will determine the resolution in terms of angle of the
array. Less resolution means you may think two waves coming in from similar
angles are coming from the same direction. Equivalently, when transmitting in
a certain direction, there will be some spread of the beam around that direction. The second parameter of interest when designing an array is the side lobe
level. Having lower side lobe levels is important for controlling the signal to
noise ratio in your system. If an interferer has a very strong signal at the side
lobe level, it could drown out the main signal. The ideal weighting scheme will
cause the array have an infinitely small beamwidth, while having zero side lobes.
However, you have seen through windowing that this is impossible, and in fact
the smaller the main lobe the higher the side lobe levels. This lab will explore
the different signal processing windows and show their usefulness in array design.
5
Lab Steps
1. The magnitude response of the ULA is equivalent to the Fourier Transform
of the weights. Therefore, we will use the various windowing functions you
have seen in order to create beamformed images
2. Open up MATLAB (I recommend script, but you may simply use the
terminal). Use rectwin(N) in order to design a rectangular window of
length N. (This is simply a vector of ones, but the syntax is the same for
all windows).
3. Using the custom MATLAB function [azm, mag] = arrayResponse(w)
from Canvas, input the window weights as w to the function. The received parameters are azm, or the angles of the response in degrees, and
mag, the magnitude response at those angles.
4. Using the subplot command, first plot the response of the rectangular
window for length N = 21 using the normal plot(azm,mag) command.
5. In the next subplot plot the polar response of the rectangular window,
polarplot(theta,rho). Note that the azm returned from arrayResponse is
in degrees, while polarplot expects the angle to be in radians. polarplot
also expects rho to be in linear scale, while arrayResponse returns the
magnitude in dBV.
Call your TA over to verify that your program is working
6. Fill out the table in Question 1, using Rectangular, Chebyshev, BlackmanHarris, Hamming, and Bartlett windows. Answer Q1 on your worksheet
7. By evaluating the performance of each of these windows, you will have
noticed that the windows always point at 0 degrees. Before electronic
beamforming, the arrays were fit with mechanical motors that would rotate the array, so that 0 degrees was always in the direction of interest.
8. The benefit of moving to electronic beamforming, as explored here, is the
ability to electronically modify the weights, and therefore the direction at
which the array operates. By changing electronically, there is much less
wear and tear on the physical components, and the angles can be switched
at a much faster rate.
9. Question 2 on your worksheet, asks you to apply a DTFT property in order to modify the angle of the main beam. HINT: The direction of the
main beam is mainly governed by the phase of the weight vector
6
Answer Q2 ond Q3 on your worksheet
10. Using the subplot command, plot in both linear and polar the response of
a Chebyshev window, with 10 dB SLL, at between 15 and 20 degrees
Call your TA over to verify that your program is working
7
3 Filtered Back Projection
This final part of the labs will return to the Radon Transform of Lab 2, and
show another application of windowing. The Radon Transform is given by
R(L) = Z
L
f(x)dx
When biomedical sensing measurements are taken, the received signal is given by
R(L). Therefore, by getting measurements around all lines through the object,
the hope is to reconstruct f(x). f(x) describes the part of the body that is
obstructed from view, and could contain information on whether someone has
cancer, a stroke, or a concussion.
8
The figure shows the common way to represent Radon Transform measurements, where the two axes represent the two parameters that uniquely identify
the line (distance from the origin and angle with from the y-axis). The first step
in reconstructing f(x) is set g(r, θ) = R(L). Then take the Fourier Transform
of g(r, θ) with respect to the r variable:
G(ρ, θ) = Z ∞
−∞
g(r, θ)e
−j2πρrdr
One of the amazing facts about the Radon Transform is the Projection-slice
Theorem. This theorem states that:
G(ρ, θ) = F1D{g(r, θ)}
F(u, v) = F2D{f(x)}
G(ρ, θ) = F(ρcos(θ), ρsin(θ))
f(x) = F
−1
2D {F1D{g(r, θ)}}
Amazingly, the function f(x) can be found by taking the 1D Fourier transform
of the measurements g(r, θ), making a small change of variable: u = ρcos(θ)
and v = ρsin(θ), and taking the inverse 2D Fourier transform.
The last step in the formula can be written explicitly as:
f(x) = 1
4π
2
Z π
0
Z ∞
−∞
F1D{g(r, θ}e
jρ(xcos(θ)+ysin(θ)
|ρ|dρdθ
The 1
4π2 comes from the fact that the inverse Fourier transform has a factor of
1
2π
out front, and this a two dimensional inverse transform. The |ρ| comes from
standard calculus change of coordinates formula (Jacobian) for integrals from
cartesian to polar. The |ρ| actually causes challenges when trying to reconstruct
the image. Since G(ρ, θ) and |ρ| are multiplied in the frequency domain, they
are equivalently convolved in the time domain, g(r, θ) ∗ φ(r). Changing to the
time domain is convenient, since if the two signals can be convolved efficiently,
there is no need to do an expensive 2D inverse Fourier transform. If such a
function could be found, then the expression for reconstruction is:
f(x) = 1
2π
Z π
0
[g(r, θ) ∗ φ(r)]dθ
Unfortunately there is no such function that satisfies this condition
Fφ(r) 6= |ρ|
However, there are several discrete approximations to this high pass filter, and
that is where windowing will make an appearance. The most commonly used
9
filter, and the filter that is the standard in MATLAB’s inverse radon command,
is the Ram-Lak Filter:
φ(rn) =



1
4d2 n = 0
1
(nπd)
2 n = odd
0 n = even
In the expression, the function is sampled at r = nd. Since this filter is a high
pass filter, it can be extremely susceptible to noise. That is where the windowing
comes in. The window works by multiplying the filter. Using a window such as
Hann or cosine window can help make the algorithm robust to noise, while still
giving good reconstruction of f(x).
10
Lab Steps
1. To begin, start with the Shepp-Logan phantom, with size 128×128, P =
phantom(N).
2. Next, take the radon transoform using 180 angles. th = linspace(0,180,180),
R = radon(P,th)
3. Add white gaussian noise to the transformed signal, so that it has a 10
dB SNR. R = awgn(R,10,’measured’).
4. Now take the inverse Radon transform, without using any filtering I =
iradon(R,th,’linear’,’none’), and find the psnr comparison between the recovered image and the original psnr(I(1:N,1:N),P)
Call your TA over to verify that your program is working
5. Question 4 will ask you to fill out a table when there is no window used,
the Ram-Lak, Shepp-Logan, Hamming, and Hann windows are used.
Answer Q4 on your worksheet
6. Change the number of sensors used (angles that are measured) from 180
to 20. Is windowing more effective or less effective when fewer sensors are
used?
Answer Q5 on your worksheet
11