Beamforming enables you to find problems in your system before they
occur. As Martin explains， you can perform your own investigation
with a computer， a dual-channel A/D card， a pair of microphones，
and some software.
Tensions are high and a hostile submarine lurks in the waters just
off the coast. A patient exhibits troubling symptoms-do layers of
tissue mask a tumor？ A driver checks his mirrors and begins to
change lanes but never sees the truck hiding in his blind spot.
It‘s better to go looking for trouble than to have it find you. But
how do you do that？ Whether it's subhunting sonar， medical
ultrasound，or collision-avoidance radar， how do you uncloak danger？
One important arrow in your quiver is beamforming. Think of it as a
way to focus and measure the energy falling on an array of sensors
as you "look" in different directions. You can do your own handson
investigation with only a computer，a dual-channel A/D card， a pair
of microphones， and the software developed in this article.
Beamforming techniques come in two flavors： adaptive and fixed.
Adaptive methods are haute cuisine. Fixed methods are burgers with
fries. Each has its place， but blue-collar fare will be served up
here.
TIME DOMAIN BEAMFORMING
To understand how beamforming works， picture a line of identically
spaced microphones. Next， imagine a wavefront consisting of noise
plus the signal you want to detect， striking this uniform line
array. Finally，make two common simplifications，valid for many
applications. One，assume that the wavefront originated in the far
field so that it can be approximated as a flat surface propagating
toward the array and any curvature it possessed near the source can
be disregarded. Second， ignore the three-dimensional quality of
wave propagation and restrict your attention to the two-dimensional
plane. With this in mind， Figure 1 shows what it would look like if
you wanted to focus the energy arriving perpendicular （broadside）to
the array.
Each phone is excited by the same wavefront at the same time. To
get the total energy falling on the array at any instant， sum
together the data from each phone. Refer to Figure 2 if you want to
focus the energy from some other direction.
To sum the energy across a given wavefront， you need to insert
delays because the wavefront has to travel a different distance to
strike each phone. The extra distance the wavefront travels from
one phone to the next is d sin θ， where d is the uniform spacing
between the phones and θ，the look direction， is by convention
measured clockwise from broadside to the array. If the wavefront
travels through the medium（e.g.， air or water） at speed c， then the
difference between the time it strikes one phone and the next is （d
sin θ）/c. Generalizing， from Figure 2 the wavefront strikes phone p
at （pd sin θ）/c units of time before it strikes the reference
phone. So， for phone p and time t， if you call the phone‘s response
yp（t）， then the time-shifted summed output of the line array is
given by：
where P is the number of phones in the line array， numbered from 0
to P - 1.
Great， you're done. Just sweep around the array and find the
direction that gives the maximum response， right？ Actually， in
practice，this time domain， delay， and sum beamformer has drawbacks
because the data acquired is sampled， not continuous.
For an ADC clocking in data at some constant interval Δt， t can
assume only values that are integer multiples of Δ。 Of course， that
also means y is known only at integer multiples of the sample
interval as well. Equation 1 says you want the value read from each
phone （the y values） at time：
With Δt， d， and c fixed， and n and k restricted to be integers，this
limits the available choices for the look direction， θ。What should
you do？
PHASE SHIFT BEAMFORMING
If you‘re looking for periodic signals， like the vibration produced
by rotating machinery，you can use the Fourier transform and work
with Equation 1 in the frequency domain instead of the time domain.
（If you're not familiar with the Fourier transform，you can treat it
as a black box whose output provides the magnitude and phase of the
frequencies inherent in a time domain input.）
To begin， let the frequency-dependent function Yp（ω） be the Fourier
transform of the time-dependent function yp（t）。 Then， according to
the shifting property， the Fourier transform of：
is as follows：
Equation 1 transforms as：
The variables in Equation 5 are continuous. However， because your
data is sampled， in practice you will apply the discrete Fourier
transform with N bins， and for a representative bin n the discrete
form of Equation 5 is：
where fS is the sampling rate. By doing this you‘ve turned the time
delay beamformer into a phase shift narrowband beamformer and you
can point the steering vector wherever you want. The
Narrowband_Beamformer. sce file posted on the Circuit Cellar FTP
site shows how to apply Equation 6 using Scilab， a freeware
alternative to Matlab.
The five-channel， 5-kHz sampled wavefile （Narrowband.wav，posted on
the Circuit Cellar FTP site） contains noise plus synthesized data
for a 1，250-Hz narrowband signal traveling through a medium at
1，500 m/s （a typical sound propagation speed in ocean water），
falling on a five-element linear array with uniform spacing of 0.6
m between phones， at an angle of arrival of 30°。 （The phone spacing
for this example was chosen deliberately. Ideally， you‘d like the
spacing to be half a wavelength or less.） When that data is
processed using the code in the Narrowband_Beamformer. sce file
posted on the Circuit Cellar FTP site and the maximum magnitude of
the beamformer output （by frequency） is plotted for look directions
around a circle in 3° increments beginning at 0°， the result is
shown in Figure 3. The absolute maximum occurs when the array is
steered to point at the source-just what you want！ As a bonus， the
use of an array improves the signal-to-noise ratio because the
noise component tends to average out as more phones are added. But
there's more to the story. Note the shape of the beamformer output.
It‘s telling you that if you point a little bit away from the
source， your beamforming will leave the waves slightly out of
phase， but you'll still pick up a good deal of power.
But as you point farther and farther away from the source， the
waves become increasingly out of phase and add
destructively，causing the power received to drop rapidly. Keep
going and the waves come back into phase a bit， adding more
constructively before going back out of phase and adding
destructively again. This process repeats as you look around a
circle centered on the array.
The same idea holds if you fix the look direction and move the
source around the array instead. This is the basis for plots called
the beam pattern or radiation pattern of an array. The variables in
these plots are the number of elements in the array， the location
of the elements， the frequency of the wave，the angle of arrival of
the wavefront， and the look direction. For a simple uniform line
array， the beam pattern b（θ） can be expressed as：
θ
S is the steering angle and λ = c/f is the wavelength. Using
Equation 7，the beam pattern for the array in the
Narrowband_Beamformer.sce file on the Circuit Cellar FTP site for a
wave at 1，250 Hz is plotted in Figure 4.
点击查看Figure 4The large teardrop-shaped curve in the look direction is called the
mainlobe. The other curves are called sidelobes. When plotting beam
patterns for a line array， you need to deal with factors of 0/0.
L‘Hospital's rule（remember your calculus！） comes to your rescue
here. Using it you can show that this ratio equates to 1 for the
case at hand.
Why should you care about the beam pattern？ Because narrowband
tones may fall on your array from more than one direction. It can
happen if the wavefront is reflected and takes multiple paths to
your array. It can happen if there is more than one source present.
In any event， these undesired interferers can feed into the
sidelobes and adversely affect the output from the beamformer.
Depending on your application， you may be able to ignore the
problem caused by sidelobes. But if you can‘t，there are methods for
minimizing them called array shading that apply weights to the data
on a phone by phone basis.
Let me make a final point concerning Figure 3 and Figure 4
regarding the symmetry they possess. From the standpoint of a line
array for a given angle， the difference between the times a
wavefront strikes any particular pair of phones is the same
regardless of the side of the array from which the wavefront
arrives. For this reason， depending on your application，you might
want to consider using a different phone configuration or using two
nonparallel line arrays.
What if the wavefront consists of broadband data and you can‘t
resolve narrowband tones （e.g.， the data in the
BroadbandStochastic.wav file on the Circuit Cellar FTP site
consists of noise）？ The parameters are the same as those used with
the Narrowband.wav file on the Circuit Cellar FTP site， except that
the sample rate is 44.1 kHz and the phone separation is about 20.41
m. Making those adjustments to the code in the
Narrowband_Beamformer.sce file on the Circuit Cellar FTP site and
running it against this new file produces the output in Figure 5.
点击查看Figure 5
Although this wavefront is also falling on the array at an angle of
30°that's not what the phase shift beamformer is telling you. What
do you do？ The solution is to turn to the cross-correlation
beamformer.
CROSS-CORRELATING
Beamforming by cross-correlating works by comparing the wavefront
received at pairs of phones. The data from one phone is fixed in
time while you slide the data from the second phone past it. For
discrete data， the sliding is really stepping， where each step
represents one tick of the sample clock of your ADC （see Figure 6）。
At each step， calculate the similarity of the two data sets. The
goal is to find the number of steps you have to shift the second
data set to get the best agreement with the first data set. Once
you know the shift as a lead or a lag （often nominally called a lag
in either case）， it‘s easy to determine the angle of arrival of the
wavefront.
The first question is how do you measure similarity. That's where
cross-correlation comes into the picture. For two functions f（t）
and g（t），their cross-correlation is denoted by：
which is defined as：
which for a finite set of discretely sampled data becomes：
Be careful about blindly applying Equation 10. Suppose you
collected two sequences 2， 1， 1， 1 and 1， 1， 1， 2. Applying
Equation 10 gives Figure 7. You can see that because the data sets
are finite. The result at each stage is influenced by the number of
nonzero products in the sum〞more so， in fact，than by how well the
peaks align. To adjust the result， you can divide each sum by the
size of the interrogation area（see Figure 8）。 As expected， the
maximum occurs when the peaks are aligned.
There is one drawback to Equation 10. It is computationally
expensive for all but the smallest data sets. Once again，the
solution is to use the discrete Fourier transform and work in the
frequency domain. This is called fast correlation. Instead of
generating a set of values where each value is formed by summing
products （as in Figure 7）， fast correlation requires only the
discrete Fourier transform，a single pointwise multiply of the two
transformed data sets， and an inverse discrete Fourier transform.
It only sounds longer than the direct method！
In symbols， if x and y are the data sets to be correlated， and if F
and F
-1 are the Fourier transform and inverse Fourier transform operators
respectively， then：
where * means the complex conjugate.
There is a twist. The discrete Fourier transform assumes your data
is periodic. That means that when you perform fast correlation， the
result is the same as what you get by performing direct correlation
on the two periodic sequences in Figure 9a. The periodic extension
of the actual data bleeds into the sum of products and you＊ve got a
mess. This is called circular correlation （because you can
equivalently think of the actual samples of the first data set as
wrapping around in a circle as they march past the actual samples
in the second data set）。 If you want to match the results of direct
correlation （also called linear correlation），you need to zero-pad
your input data （see Figure 9b）。
But zero-padding is more work to code（in real-time C/C++ code，
anyway）。 In general， if one data set is N samples long and the
other set is M samples long， to correlate in the frequency domain
you must tack zeros to the end of both sets so that each set has
length N + M - 1. Then，following the inverse transform， you
need to select the proper subset. More importantly，padding makes
for larger， slowerto- execute Fourier transforms. So skip it！ You
can if the worst lag possible（which occurs when the wavefront
arrives from endfire） is small relative to the number of sample
counts by which the phones are separated. That＊s because over the
interval of interest，circularity will come into play for only a few
samples （see Figure 10）。 Add to that the fact that windowing
typically precedes the application of the Fourier transform，
tapering the ends of the sequences toward zero， and you can see
that the effects of circularity are further minimized.
Once you‘ve determined the lag， all that's left is to turn it into
a bearing. That‘s easy. As you can see in Figure 2，the time
difference of arrival could be expressed as：
Δt=（d sin θ）/c [12]
But note that：
Δt=L/f
s [13]
where fS is the sample frequency and L is the number of lags
required to align the data received at the relative phone with the
data received at the reference phone. （A delay is a positive shift
and an advance is a negative shift.） So， the bearing is given by：
The Broadband_Beamformer.sce file on the Circuit Cellar FTP site
pulls these ideas together with Equation 11 to calculate a
broadband bearing.
When this code is run against the BroadbandStochastic.wav file on
the Circuit Cellar FTP site， the result is shown in Figure 11. The
peak occurs at a shift of 300 sample counts， which translates into
a bearing close to the actual value of 30°。
点击查看Figure 11Typically， when beamforming via cross-correlation， you can use
multiple phone pairs-for N phones you can form N（N - 1）/2 phone
pairs-and add the results from each. This gives a better estimate
of direction of arrival. It is also especially valuable if the
phones are not collinear because it eliminates bearing ambiguity.
The true bearings will add constructively while the false bearings
will add randomly and contribute to the noise.
If there is a complaint to be had about the cross-correlation
beamformer，it is that it will be fooled by a noticeable narrowband
tone within a stochastic broadband signal. But there‘s a solution
for that too！
A slight modification to the crosscorrelation beamformer makes a
big difference in the result. After multiplying the Fourier
transform of the two input sequences， the phase transform，or PHAT，
divides the complex value in each frequency bin by its
magnitude.（This is equivalent to applying a filter to the original
time series data because multiplication in the frequency domain
corresponds to convolution in the time domain.） With the magnitudes
equalized， the correlation then relies on only the phase of the
input sequences. Interfering narrowband tones are neutralized. To
see what a dramatic difference this can make，refer to Figure 12
（cross-correlation without the PHAT） and Figure
13（cross-correlation with the PHAT）。The input file
（BroadbandStochastic- PlusTone.wav on the Circuit Cellar FTP site）
is the same as BroadbandStochastic. wav， except for an added
narrowband tone at a bearing of 64°。
点击查看Figure 12点击查看Figure 13Despite their low-brow status， fixed beamformers make their way
into commercial products， where they prove their worth. Here‘s
hoping they can point you in the right direction too， so the next
time you go looking for trouble， you find it！
Formally trained in mathematics，Martin Courtney has been writing
software for more years than he cares to admit. Although his
preference is signal processing and algorithm development， he has
written code for real-time embedded systems and graphically
intensive PC-based applications as well. You can reach Martin at
mcourtn1@san.rr.com.
PROJECT FILESTo download additional files, go to
ftp://ftp.circuitcellar.com/pub/Circuit_Cellar/2008/221.
RESOURCESW. S. Burdic, Underwater Acoustic System Analysis, Peninsula
Publishing,Los Altos, CA, 2002.
R. V. Churchill, Operational Mathematics,McGraw-Hill, New York,
NY,1972.
Scilab, www.scilab.org.
;){kind=link}