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
0 comments