Target direction indication and acoustic pulse analysis

[0001] The present invention relates in a first aspect to a method and apparatus for providing an indication of the direction of a target, and in a second aspect to a method and apparatus for estimating the direction of arrival of an acoustic pulse. Both aspects may be employed typically, although not exclusively, in a system for identifying the direction and/or range of a sniper and/or determining the trajectory of a bullet fired by the sniper.

[0002] U.S. Pat. No. 6,178,141 describes a system for determining the trajectory of a supersonic bullet by observing the shock wave from the bullet. Muzzle blast observations may also be used to estimate the exact sniper location along the trajectory. The sensors are helmet-mounted, and each helmet is supported by wireless communications, GPS and a backpack or otherwise mounted computer system that computes and displays the solutions.

[0003] A first aspect of the invention provides a method of providing an indication of the direction of a target, the method comprising the steps of detecting a signal originating from the target; processing the signal to estimate the direction of the target relative to a weapon; and providing an indication when the weapon is aligned with the estimated direction of the target.

[0004] The first aspect of the invention also provides apparatus for providing an indication of the direction of a target, the apparatus comprising a sensor for detecting a signal originating from the target; a processor for processing the signal to estimate the direction of the target relative to a weapon; and a direction indicator for providing an indication when the weapon is aligned with the estimated direction of the target.

[0005] The first aspect of the invention presents a new method and apparatus for presenting target direction information to a user. By providing a user with an indication when a weapon (typically held or controlled by the user) is aligned with the estimated direction of a target such as a sniper, the user can quickly and accurately return fire.

[0006] Preferably the apparatus comprises an orientation sensor for providing an output indicative of the direction of the weapon, and the direction indicator provides an indication when the output from the orientation sensor indicates that the weapon is aligned with the estimated direction of the target.

[0007] The signal originating from the target may be a visual signal (such as a visual muzzle flash) but most preferably the signal is an acoustic signal such as a muzzle blast or projectile shock wave.

[0008] The signal may be detected by helmet-mounted sensors as in U.S. Pat. No. 6,178,141. However in some situations it may be preferable for the user to be operating in a “helmets off” environment. Also, if the sensors are helmet-mounted then a complex system must be provided to estimate the target direction relative to the weapon. Therefore preferably the signal is detected by one or more sensors mounted on the weapon.

[0009] Preferably the indication is also provided by an indicator mounted on the weapon.

[0010] Preferably the sensor, indicator and processor are mounted in a common housing.

[0011] In a preferred embodiment an indication is provided when the weapon is aligned with a vertical and/or horizontal plane containing the estimated target direction.

[0012] The indication may be visual or audible. In a preferred example the colour of a visual indicator is changed when the weapon is aligned with the estimated direction of the target.

[0013] The signal is typically detected by a plurality of sensors which may be arranged in a planar array, or in a volumetric array. Preferably the spacing (or baseline) between the sensors is small, in order to reduce the size of the apparatus. Small spacing does however make it difficult to accurately estimate the target direction. Therefore preferably one of the following methods is used to address this problem:

• • [0014] i. using a conventional time-of-arrival method in which a threshold crossing point is detected for each sensor. In this case the sampling frequency can be increased to take account of the reduced sensor baseline; or
  • [0015] ii. cross-correlating one or more pairs of sensor signals in the time domain so as to generate one or more cross-correlation signals which are then used in the target direction estimation; or
  • [0016] iii. spectrally decomposing each sensor signal so as to generate one or more spectral components for each signal which are then used in the target direction estimation.

[0017] Method iii) is preferred, and is discussed in further detail below with reference to the second aspect of the invention.

[0018] The weapon may be a tank or helicopter, but preferably the weapon is a hand-held weapon such as a rifle.

[0019] Preferably the (or each) sensor is mounted on the sight of the weapon. This automatically aligns the sensors with the sight and enables more accurate targeting.

[0020] The direction and/or range indicators may also be mounted on the sight to provide a compact integrated system.

[0021] In U.S. Pat. No. 6,278,141, fine resolution time difference of arrivals, for the shock wave arrivals, are determined via cross-correlation between the channels. The data are interpolated eight times during the cross-correlation process.

[0022] A second aspect of the invention provides a method of estimating the direction of arrival of an acoustic pulse, the method comprising:

• • [0023] a. detecting the pulse with a plurality of sensors so as to generate a plurality of sensor signals;
  • [0024] b. spectrally decomposing each sensor signal so as to generate one or more spectral components for each signal; and
  • [0025] c. processing the spectral components to estimate the direction of arrival of the acoustic pulse.

[0026] The second aspect of the invention also provides apparatus for detecting the direction of arrival of an acoustic pulse by the above method, the apparatus comprising an array of sensors; and a processor programmed to perform steps b and c.

[0027] The second aspect of the invention provides a processing method which enables the sensor signals to be sampled at a low rate (which may be as low as the Nyquist rate), whilst maintaining an accurate estimate. Spectrally decomposing the sensor signal reduces the effects of noise, compared with a conventional time-of-arrival method in which a threshold crossing point is detected for each sensor. The second aspect of the invention is particularly suited to use in a volumetric sensor with a small baseline.

[0028] In one example of the second aspect of the invention, the direction of arrival is calculated by performing a power versus bearing scan. This provides improved accuracy through “bearing coherent” averaging of noise. However for some applications this method may be overcomplicated. Therefore in a preferred embodiment the pulse is detected with a volumetric array of three or more pairs of sensors in step a. so as to generate three or more pairs of sensor signals; a plurality of pairs of spectral components are generated for each pair of sensor signals in step b; and step c. includes calculating a phase delay between each pair of spectral components, and determining the direction of arrival of the acoustic pulse in accordance with the calculated phase delays. The preferred embodiment is computationally simpler than a power versus bearing scan, and also provides improved accuracy through “bearing coherent” averaging of noise.

[0029] The method may be used in any suitable application, but in preferred examples it is used to estimate the direction of a target and/or to estimate the trajectory of a projectile.

[0047] FIG. 1 shows a rifle 1 with a sight 2, and an acoustic sensor device 3. The device 3 is shown in further detail in FIG. 2. The device 3 has a pair of mounting legs 4 which are screwed in place onto the sight 3 and a housing comprising upper and lower hemispheres 5,6. A digital range display 7 is mounted on the lower hemisphere 5, in this case showing a range of 800 m. A support plate 8 is screwed onto the lower hemisphere 5. Three light emitting diodes (LEDs) 9 are mounted in the lower hemisphere 5 and one in the upper hemisphere 6

[0048] FIG. 3 is a schematic view showing various component parts of the device 3, some of which are omitted from FIG. 1 for clarity. A volumetric array 20 of six miniature microphones is mounted to the support plate 8. The microphones are arranged in orthogonal positions around a sphere, with a diameter of between about 4 and 6 cm. The microphones are digitally sampled by an analogue to digital converter 21 and processed by an embedded digital signal processor (DSP) 22. A pseudo-omni output is formed for automatic detection of impulsive noise. The relative phases of all six microphones are then used to estimate the 3-dimensional direction of arrival of the detected sound. Systematic errors associated with the construction or mounting of the device (e.g. due to resonance, coupling or shading) can be precalibrated. Adaptive processing can be used to reject noise and interference.

[0049] The DSP 22 controls the LEDs 9, and range display 7, and also outputs data to a Bluetooth antenna or other wireless data link.

[0050] As shown in FIG. 3a, the Bluetooth antenna communicates with a personal digital assistant (PDA) 23 which is carried by the soldier. The PDA 23 also receives GPS data 24, and displays range and direction information in a map-based format. The PDA 23 is connected to a radio communication device 25 via an RS232 link.

[0051] The DSP 22 estimates the direction of the muzzle blast pulse, and the direction of the shock wave pulse. The range is then estimated by the DSP 22 using the method illustrated in FIG. 4. The shot is fired from the left hand corner, and is detected at the right hand corner. The muzzle blast travels directly between the target and the receiver. The shock wave from the passing bullet arrives from a different angle. By measuring the included angle, and the time delay between the two events, and making assumptions about the speed and deceleration of the bullet, a range estimate can be computed.

[0052] FIG. 3B is a schematic view of the components of an alternative acoustic sensor device. The device is similar to the device of FIG. 3a but includes an orientation sensor 23. The orientation sensor 23 contains an inclinometer (which measures the elevation angle of the device relative to the gravity vector) and a compass (which measures the azimuth angle of the device relative to the earth's magnetic field vector). Thus the sensor 23 provides an output indicative of the sensor's direction in space. The device of FIG. 3B operates as follows:

• • [0053] 1. The sensor array 20 outputs an estimated range and direction of a target to the DSP 22;
  • [0054] 2. The range display outputs the estimated range;
  • [0055] 3. The orientation sensor 23 senses the current orientation of the device, and outputs to the DSP;
  • [0056] 4. As the rifle is re-oriented, the DSP compares the estimated direction of the target with the current orientation of the device;
  • [0057] 5. If the device is aligned with a horizontal plane containing the estimated target direction, then the two LEDs 9 on the left and right-hand side of the range display shown in FIG. 2 change from red to green;
  • [0058] 6. If the device is aligned with a vertical plane containing the estimated target direction, then the two LEDs 9 above and below the range display shown in FIG. 2 change from red to green.

[0059] The comparison in step 4 remains valid whilst the translational position in space of the device is close to the position where the sensor array 20 measured the range and direction. This enables the device to give an immediate detection and localization of a target such as a sniper (for example to a particular window or building) after which usual tactics would apply.

[0060] The networked system shown in FIG. 3a enables the device to be coupled with GPS data 24 which can be used to fix the translational position of the device, and thus update the estimated range and direction of the target as the device is translated.

[0061] Extraction of the direction of the target involves two essential steps: (i) measure the propagation delay across the array, and (ii) infer the direction of arrival.

[0062] A feature of the sensor array is that it includes three orthogonal axes each with two microphones.

[0063] Therefore, direction can simply be extracted by computing the time delay component along each axis, and resolving them to form a direction vector. This direct form of direction estimation is simpler than say a full maximum-likelihood power vs bearing scan of the array, with almost as good results for high SNR pulses.

[0064] The preferred method of calculating direction will now be contrasted with other methods which may also be performed by the DSP 22. Three candidate methods are considered in detail:

• • [0065] a) Measure the time of arrival of the pulse at each microphone using a threshold crossing technique. Form the direction vector from the three component delays.
  • [0066] b) Cross-correlate each pair of microphones in the time domain, to estimate the delay based on the whole pulse shape. Form the direction vector from the three component delays.
  • [0067] c) The preferred method: decompose the pulse into its spectral components; compute the phase delay for each component along each axis; hence create a direction spectrum (ie. power-weighted direction vector vs frequency); and sum the weighted direction vectors, to obtain the overall direction vector.

[0068] The relative advantages and disadvantages of each technique are discussed below.

Method ‘a’

[0069] The array is small in terms of wavelengths (less than two wavelengths across) which means that if using a direct time-domain measurement method the signals must be sampled at a much higher rate (e.g. 400 kHz vs 20 kHz). This will increase power drain and increase the size of the electronics.

Method ‘b’

[0070] The use of correlation-based time measurement means that the necessary high sample rate can be obtained post-correlation, by interpolation. Method ‘b’ also has the advantage over method ‘a’ that using the whole pulse shape provides averaging of the effects of noise, whereas using a threshold-crossing method yields just one time sample per event. However, for sensible length interpolation filters the basic sample rate would have to exceed the Nyquist rate by a factor.

Method ‘c’

[0071] Method ‘c’ retains the advantage of method ‘b’ by using the whole pulse shape, but avoids the need to interpolate the correlation output to a very high effective sample rate, so allowing sampling at the Nyquist rate. The phase domain method also provides slightly improved accuracy through ‘bearing coherent’ averaging of noise, rather than just temporal averaging of noise.

[0072] Some more explanation of method c. will now be given. An alternative method of using phase in the direction estimation (employing a power vs bearing scan) will also be described.

[0073] The signal vector, across all six microphones, for a signal at elevation θ and azimuth θ is as follows:

[0074] $s\left(\theta ,\psi ,\lambda \right)=\left(\begin{array}{c}{e}^{j·2·\pi ·\frac{r}{\lambda }·\cos \left(\psi \right)·\cos \left(\theta \right)}\\ e^{j·2·\pi ·\frac{r}{\lambda }·\sin \left(\psi \right)·\cos \left(\theta \right)}\\ e^{-j·2·\pi ·\frac{r}{\lambda }·\cos \left(\psi \right)·\cos \left(\theta \right)}\\ e^{-j·2·\pi ·\frac{r}{\lambda }·\sin \left(\psi \right)·\cos \left(\theta \right)}\\ e^{j·2·\pi ·\frac{r}{\lambda }·\cos \left(\theta \right)}\\ e^{-j·2·\pi ·\frac{r}{\lambda }·\cos \left(\theta \right)}\end{array}\right)$

[0075] So the received signal, at true direction (θ₀ψ₀), in the presence of noise, is

[0076] $x\left(\lambda \right)=s\left({\theta }_o,{\psi }_o,\lambda \right)+\left(\begin{array}{c}{n}_0\\ n_1\\ n_2\\ n_3\\ n_4\\ n_5\end{array}\right)$

[0077] The optimal method of estimating the signal direction, assuming uncorrelated background noise would be to do a power vs bearing scan. ie. for beamweights w given by:
w(θ,ψ,λ)=( s(θ,ψ,λ))^T
we find the θ and ψ that maximise the following expression
P=w·X· w^T
where
X=x(λ)· x(λ)^T

[0078] However this is overcomplicated for the preferred application.

[0079] Let us consider a modified signal vector, in which opposite microphones along each axis are cross-correlated (or equivalently, the complex frequency-domain signals are multiplied in conjugate).

[0080] E.g. taking microphones 0 and 2 we have our correlated x-axis sensor:

[0081] $\begin{array}{c}{a}_x=s_0·\stackrel{\_}{s_2}\\ =e^{j·2·\pi ·\frac{r}{\lambda }·\cos \left(\psi \right)·\cos \left(\theta \right)}·\stackrel{\_}{e^{-j·2·\pi ·\frac{r}{\lambda }·\cos \left(\psi \right)·\cos \left(\theta \right)}}\\ =e^{j·4·\pi ·\frac{r}{\lambda }·\cos \left(\psi \right)·\cos \left(\theta \right)}\end{array}$ $\mathrm{Therefore}$ $a\left(\theta ,\psi ,\lambda \right)=e^{j·\varphi \left(\theta ,\psi ,\lambda \right)}$ $\mathrm{where}$ $\varphi \left(\theta ,\psi ,\lambda \right)=\left(\begin{array}{c}4·\pi ·\frac{r}{\lambda }·\cos \left(\psi \right)·\cos \left(\theta \right)\\ 4·\pi ·\frac{r}{\lambda }·\sin \left(\psi \right)·\cos \left(\theta \right)\\ 4·\pi ·\frac{r}{\lambda }·\cos \left(\theta \right)\end{array}\right)$

[0082] But the direction vector of the true signal (in x y z coordinates) is as follows:

[0083] $d\left(\theta ,\psi \right)=\left(\begin{array}{c}\cos \left(\psi \right)·\cos \left(\theta \right)\\ \sin \left(\psi \right)·\cos \left(\theta \right)\\ \cos \left(\theta \right)\end{array}\right)$

[0084] Therefore we see the following relationship, from which signal direction is obtained directly:

[0085] $\varphi \left(\theta ,\psi ,\lambda \right)=4·\pi ·\frac{r}{\lambda }·d\left(\theta ,\psi \right)$

[0086] Therefore, in words, the phase of the cross-correlated microphone pairs, directly yields the signal direction vector, to within a frequency-dependent scale factor. Therefore by analysing the phase of the received pulse, we can estimate the bearing of each spectral component. We can then produce an overall bearing estimate from the power-weighted sum of the component direction vectors. This is not quite optimal (compared to a power bearing scan) but is simpler. Also it does not require any of the microphones to be matched in gain, and only requires them to be phase-matched in pairs, which simplifies construction or calibration of the array.

Example Calculations

[0087]

Maximum frequency	fmax := 10000		Hz
Speed of sound in air	c := 344		m/s
Acoustic radius of array	r := 0.025		m
Acoustic wavelength	λmin := c f ⁢ ⁢ max	λmin = 0.0344	m
Acoustic diameter in wavelengths	2 · r λ ⁢ ⁢ max = 1.4535
Time delay across array	τ := 2 · r c	τ · 106 = 145.35	μs
Time delay for 1 degree bearing	τ 1 := τ · sin ⁢ ⁢ ( 1 · π 180 )	τ1 · 106 = 2.536687	μs
Sample rate for Nyquist sampling	fsNyq := 2 · fmax	fsNyq · 10−3 = 20	kHz
Sample rate to resolve 1 degree change	fs 1 ⁢ ⁢ deg := 1 τ 1	fs1deg · 10−3 = 394.21	kHz

Three Algorithms Compared

[0088] The three algorithms are compared using a simple simulation.

[0089] A band-limited pulse is constructed.

[0090] An uncorrelated noise background is obtained by phase-randomising band-limited pulses, to ensure simulation consistency regardless of sample rate.

[0091] The scenario is run 100 times, with a different true pulse bearing each time, and the resulting accuracy is compared for the three methods.

[0092] Sample rates and data collection block lengths for the three methods:

[0093] $T_{\mathrm{block}}:=\frac{16}{f\max }$ $T_{\mathrm{block}}·{10}^3=1.6$ ${\mathrm{fs}}_a:=32·f\max$ ${\mathrm{fs}}_b:=4·f\max$ ${\mathrm{fs}}_c:=2·f\max$ ${\mathrm{fs}}_a·{10}^{-3}=320$ ${\mathrm{fs}}_b·{10}^{-3}=40$ ${\mathrm{fs}}_c·{10}^{-3}=20$ $N_a:=\mathrm{round}\left(T_{\mathrm{block}}·{\mathrm{fs}}_a\right)$ $N_b:=\mathrm{round}\left(T_{\mathrm{block}}·{\mathrm{fs}}_b\right)$ $N_c:=\mathrm{round}\left(T_{\mathrm{block}}·{\mathrm{fs}}_c\right)$ $N_a=512$ $N_b=64$ $N_c=32$

[0094] Nominal pulse, and (approximately) equal energy flat-spectrum pulse used to generate noise.

[0095] We assume that the pulse is windowed or gated, so as to exclude extraneous noise.

[0096] $f_o:=\frac{f\max }{2}$ $\mathrm{sync}\left(x\right):=\mathrm{if}\left(x=0,1,\frac{\sin \left(x\right)}{x}\right)$ $\mathrm{win}\left(t\right):={\sin \left(\pi ·\frac{t}{T_{\mathrm{block}}}\right)}^4$ $p\left(t\right):=\begin{array}{c}{\mathrm{sync}\left[2·\pi ·\frac{f_o}{2}·\left(t-\frac{T_{\mathrm{block}}}{2}\right)\right]}^2·\\ \sin \left[2·\pi ·f_o·\left(t-\frac{T_{\mathrm{block}}}{2}\right)+\frac{\pi }{4}\right]\end{array}$ $q\left(t\right):=\begin{array}{c}\mathrm{sync}\left[2·\pi ·f_o·\left(t-\frac{T_{\mathrm{block}}}{2}\right)\right]·\\ \cos \left[2·\pi ·f_o·\left(t-\frac{T_{\mathrm{block}}}{2}\right)\right]\end{array}$ $n:=0\dots N_a-1$ $x_n:=p\left(n·\frac{1}{{\mathrm{fs}}_a}\right)$ $X:=\mathrm{FFT}\left(x\right)$ $y_n:=q\left(n·\frac{1}{{\mathrm{fs}}_a}\right)$ $Y:=\mathrm{FFT}\left(y\right)$ ${\mathrm{SNR}}_{\mathrm{adjustdB}}:=20·\log \left(\frac{\mathrm{stdev}\left(x\right)}{\mathrm{stdev}\left(y\right)}\right)$ ${\mathrm{SNR}}_{\mathrm{adjustdB}}:=1.4118$

[0097] Common simulation parameters and functions

• [0098] :=1000
• [0099] k:=0 . . . K−1

[0100] ${\theta }_k:=2·\pi ·\frac{k}{K}$

• [0101] SNR_dB:=16

[0102] $\mathrm{getpulse}\left(N,\mathrm{fs},\mathrm{delay}\right):=|\begin{array}{c}\begin{array}{c}\mathrm{for}n\in 0\dots N-1\\ x_n\leftarrow p\left(n·\frac{1}{\mathrm{fs}}-\mathrm{delay}\right)·\mathrm{win}\left(\frac{n}{N}·T_{\mathrm{block}}\right)\end{array}\\ x\end{array}$ $\mathrm{getnoise}\left(N,\mathrm{fs},{\mathrm{SNR}}_{\mathrm{dB}}\right):=|\begin{array}{c}\mathrm{for}n\in 0\dots N-1\\ |\begin{array}{c}{x}_n\leftarrow q\left(n·\frac{1}{\mathrm{fs}}\right)\\ w_n\leftarrow \mathrm{win}\left(\frac{n}{N}·T_{\mathrm{block}}\right)\end{array}\\ X\leftarrow \mathrm{CFFT}\left(x\right)\\ \varphi \leftarrow \mathrm{runif}\left(N,-\pi ,\pi \right)\\ \mathrm{for}n\in 0\dots N-1\\ Y_n\leftarrow X_n·e^{j·{\varphi }_n}\\ y\leftarrow \mathrm{Re}\left(\mathrm{ICFFT}\left(Y\right)\right)·{10}^{-\frac{{\mathrm{SNR}}_{\mathrm{dB}}-{\mathrm{SNR}}_{\mathrm{adjustdB}}}{20}}\\ \overrightarrow{\left(y·w\right)}\end{array}$
Method ‘a’

[0103] Generate pulses for each of 100 bearings, for 4 microphones (0=north, 1=east, 2=south, 3=west)

[0104] ${\mathrm{pulse}}_{a0}^{\langle k\rangle }:=\mathrm{getpulse}\left(N_a,{\mathrm{fs}}_a,-\frac{r}{c}·\cos \left({\theta }_k\right)\right)$ ${\mathrm{pulse}}_{a1}^{\langle k\rangle }:=\mathrm{getpulse}\left(N_a,{\mathrm{fs}}_a,-\frac{r}{c}·\sin \left({\theta }_k\right)\right)$ ${\mathrm{pulse}}_{a2}^{\langle k\rangle }:=\mathrm{getpulse}\left(N_a,{\mathrm{fs}}_a,\frac{r}{c}·\cos \left({\theta }_k\right)\right)$ ${\mathrm{pulse}}_{a3}^{\langle k\rangle }:=\mathrm{getpulse}\left(N_a,{\mathrm{fs}}_a,\frac{r}{c}·\sin \left({\theta }_k\right)\right)$

[0105] Generate noise waveforms
noise_a0^<k>:=getnoise(N_a,fs_a,SNR_dB)
noise_a1^<k>:=getnoise(N_a,fs_a,SNR_dB)
noise_a2^<k>:=getnoise(N_a,fs_a,SNR_dB)
noise_a3^<k>:=getnoise(N_a,fs_a,SNR_dB)

[0106] Pulse detection threshold (say set at 4 sigma wrt noise, or ¼ of the peak, whichever is greater)

[0107] ${\mathrm{Thresh}}_a:=\max \left(4·\mathrm{stdev}\left({\mathrm{noise}}_{a0}\right),\frac{1}{4}·\max \left({\mathrm{pulse}}_{a0}\right)\right)$ ${\mathrm{Thresh}}_a=0.2392$ $:=0\dots N_a-1$

[0108] Compute times of arrival for each pulse

[0109] $\mathrm{getTOA}\left(x,T\right):=|\begin{array}{c}n\leftarrow 0\\ \mathrm{while}\left(x_n<T\right)·\left(n<N_a-1\right)\\ n\leftarrow n+1\\ n\end{array}$ ${\mathrm{TOA}}_{a{0}_k}:=\mathrm{getTOA}\left({\mathrm{pulse}}_{a0}^{\langle k\rangle }+{\mathrm{noise}}_{a0}^{\langle k\rangle },{\mathrm{Thresh}}_a\right)$ ${\mathrm{TOA}}_{a{1}_k}:=\mathrm{getTOA}\left({\mathrm{pulse}}_{a1}^{\langle k\rangle }+{\mathrm{noise}}_{a1}^{\langle k\rangle },{\mathrm{Thresh}}_a\right)$ ${\mathrm{TOA}}_{a{2}_k}:=\mathrm{getTOA}\left({\mathrm{pulse}}_{a2}^{\langle k\rangle }+{\mathrm{noise}}_{a2}^{\langle k\rangle },{\mathrm{Thresh}}_a\right)$ ${\mathrm{TOA}}_{a{3}_k}:=\mathrm{getTOA}\left({\mathrm{pulse}}_{a3}^{\langle k\rangle }+{\mathrm{noise}}_{a3}^{\langle k\rangle },{\mathrm{Thresh}}_a\right)$

[0110] Compute bearing error for each pulse
θerr_ak:=arg[[(TOA_a2k−TOA_a0k)+j·(TOA_a3k−TOA_a1k)]·e^−j·θk]
Method ‘b’

[0111] Generate pulses for each of 100 bearings, for 4 microphones (0=north, 1=east, 2=south, 3=west)

[0112] ${\mathrm{pulse}}_{b0}^{\langle k\rangle }:=\mathrm{getpulse}\left(N_b,{\mathrm{fs}}_b,-\frac{r}{c}·\cos \left({\theta }_k\right)\right)$ ${\mathrm{pulse}}_{b1}^{\langle k\rangle }:=\mathrm{getpulse}\left(N_b,{\mathrm{fs}}_b,-\frac{r}{c}·\sin \left({\theta }_k\right)\right)$ ${\mathrm{pulse}}_{b2}^{\langle k\rangle }:=\mathrm{getpulse}\left(N_b,{\mathrm{fs}}_b,\frac{r}{c}·\cos \left({\theta }_k\right)\right)$ ${\mathrm{pulse}}_{b3}^{\langle k\rangle }:=\mathrm{getpulse}\left(N_b,{\mathrm{fs}}_b,\frac{r}{c}·\sin \left({\theta }_k\right)\right)$

[0113] Generate noise waveforms
noise_b0^<k>:=getnoise(N_b,fs_b,SNR_dB)
noise_b1^<k>:=getnoise(N_b,fs_b,SNR_dB)
noise_b2^<k>:=getnoise(N_b,fs_b,SNR_dB)
noise_b3^<k>:=getnoise(N_b,fs_b,SNR_dB)
:=0 . . . N_b−1

[0114] Compute time delays for each axis by circular correlation and interpolation

[0115] $\mathrm{getDT}\left(x,y\right):=|\begin{array}{c}X\leftarrow \mathrm{CFFT}\left(x\right)\\ Y\leftarrow \mathrm{CFFT}\left(y\right)\\ Z\leftarrow \overrightarrow{\left(X·\stackrel{\_}{Y}\right)}\\ \mathrm{for}n\in 0\dots N_a-1\\ Z_{{\mathrm{intp}}_n}\leftarrow 0\\ \mathrm{for}n\in 1\dots \frac{N_b}{2}-1\\ |\begin{array}{c}{Z}_{{\mathrm{intp}}_n}\leftarrow Z_n\\ Z_{{\mathrm{intp}}_{N_a-n}}\leftarrow Z_{N_b-n}\end{array}\\ z_{\mathrm{intp}}\leftarrow \mathrm{ICFFT}\left(z_{\mathrm{intp}}\right)\\ \mathrm{for}n\in 0\dots N_a-1\\ |\begin{array}{c}{\mathrm{sortarray}}_{n,0}\leftarrow n\\ {\mathrm{sortarray}}_{n,1}\leftarrow z_{{\mathrm{intp}}_n}\end{array}\\ \mathrm{mod}\left({\mathrm{csort}\left(\mathrm{sortarray},1\right)}_{N_a-1,0}+\frac{N_a}{2},N_a\right)-\frac{N_a}{2}\end{array}$ ${\mathrm{DT}}_{b_{{20}_k}}:=\mathrm{getDT}\left({\mathrm{pulse}}_{b2}^{\langle k\rangle }+{\mathrm{noise}}_{b2}^{\langle k\rangle },{\mathrm{pulse}}_{b0}^{\langle k\rangle }+{\mathrm{noise}}_{b0}^{\langle k\rangle }\right)$ ${\mathrm{DT}}_{b_{{31}_k}}:=\mathrm{getDT}\left({\mathrm{pulse}}_{b3}^{\langle k\rangle }+{\mathrm{noise}}_{b3}^{\langle k\rangle },{\mathrm{pulse}}_{b1}^{\langle k\rangle }+{\mathrm{noise}}_{b1}^{\langle k\rangle }\right)$

[0116] Compute each bearing error for each pulse

[0117] $\theta {\mathrm{err}}_{b_k}:=\arg \left[\left({\mathrm{DT}}_{b{20}_k}+j·{\mathrm{DT}}_{b{31}_k}\right)·e^{-j·{\theta }_k}\right]$
Method ‘c’

[0118] Generate pulses for each of 100 bearings, for 4 microphones (0=north, 1=east, 2=south, 3=west)

[0119] ${\mathrm{pulse}}_{c0}^{\langle k\rangle }:=\mathrm{getpulse}\left(N_c,{\mathrm{fs}}_c,-\frac{r}{c}·\cos \left({\theta }_k\right)\right)$ ${\mathrm{pulse}}_{c1}^{\langle k\rangle }:=\mathrm{getpulse}\left(N_c,{\mathrm{fs}}_c,-\frac{r}{c}·\sin \left({\theta }_k\right)\right)$ ${\mathrm{pulse}}_{c2}^{\langle k\rangle }:=\mathrm{getpulse}\left(N_c,{\mathrm{fs}}_c,\frac{r}{c}·\cos \left({\theta }_k\right)\right)$ ${\mathrm{pulse}}_{c3}^{\langle k\rangle }:=\mathrm{getpulse}\left(N_c,{\mathrm{fs}}_c,\frac{r}{c}·\sin \left({\theta }_k\right)\right)$

[0120] Generate noise waveforms
noise_c0^<k>:=getnoise(N_c,fs_c,SNR_dB)
noise_c1^<k>:=getnoise(N_c,fs_c,SNR_dB)
noise_c2^<k>:=getnoise(N_c,fs_c,SNR_dB)
noise_c3^<k>:=getnoise(N_c,fs_c,SNR_dB)
:=0 . . . N_c−1

[0121] Compute spectra for each microphone
X_c0^<k>:=FFT(pulse_c0^<k>+noise_c0^<k>)
X_c1^<k>:=FFT(pulse_c1^<k>+noise_c1^<k>)
X_c2^<k>:=FFT(pulse_c2^<k>+noise_c2^<k>)
X_c3^<k>:=FFT(pulse_c3^<k>+noise_c3^<k>)

[0122] Compute bulk delays (by circular correlation, without interpolation) to a resolution of 1 sample (ie. half wavelength max) to eliminate phase ambiguity.

[0123] $\mathrm{getBulkDT}\left(x,y\right):=|\begin{array}{c}X\leftarrow \mathrm{CFFT}\left(x\right)\\ Y\leftarrow \mathrm{CFFT}\left(y\right)\\ Z\leftarrow \overrightarrow{\left(X·\stackrel{\_}{Y}\right)}\\ z\leftarrow \mathrm{ICFFT}\left(Z\right)\\ \mathrm{for}n\in 0\dots N_c-1\\ |\begin{array}{c}{\mathrm{sortarray}}_{n,0}\leftarrow n\\ {\mathrm{sortarray}}_{n,1}\leftarrow z_{{}_n}\end{array}\\ \mathrm{mod}\left({\mathrm{csort}\left(\mathrm{sortarray},1\right)}_{N_c-1,0}+\frac{N_c}{2},N_c\right)-\frac{N_c}{2}\end{array}$ ${\mathrm{BulkDT}}_{c_{{20}_k}}:=\mathrm{getBulkDT}\left({\mathrm{pulse}}_{c2}^{\langle k\rangle }+{\mathrm{noise}}_{c2}^{\langle k\rangle },{\mathrm{pulse}}_{c0}^{\langle k\rangle }+{\mathrm{noise}}_{c0}^{\langle k\rangle }\right)$ ${\mathrm{BulkDT}}_{c_{{31}_k}}:=\mathrm{getBulkDT}\left({\mathrm{pulse}}_{c3}^{\langle k\rangle }+{\mathrm{noise}}_{c3}^{\langle k\rangle },{\mathrm{pulse}}_{c1}^{\langle k\rangle }+{\mathrm{noise}}_{c1}^{\langle k\rangle }\right)$

[0124] Compute direction vector spectrum, accounting for bulk delay phase ambiguity.

Frequency Resolution:

[0125] $N_c=32$ $\delta f_c:=\frac{{\mathrm{fs}}_c}{N_c}$ $\delta f_c=625$ ${\mathrm{freq}}_c\left(n\right):=\left(\mathrm{mod}\left(n+\frac{N_c}{2},N_c\right)-\frac{N_c}{2}\right)·\frac{{\mathrm{fs}}_c}{N_c}$ $\mathrm{arrgh}\left(x\right):=\mathrm{if}\left(x=0,0\arg \left(x\right)\right)$ $\begin{array}{c}\mathrm{getDvecspect}\\ \left(\begin{array}{c}X0,X1,X2,\\ X3,\mathrm{DT}20,\mathrm{DT}31\end{array}\right)\end{array}:=|\begin{array}{c}\mathrm{for}n\in 0\dots \frac{N_c}{2}-1\\ |\begin{array}{c}f\leftarrow {\mathrm{freq}}_c\left(n\right)\\ z02\leftarrow X{0}_n·\stackrel{\_}{X{2}_n}\\ z13\leftarrow X{1}_n·\stackrel{\_}{X{3}_n}\\ {\varphi }_{\mathrm{adj}20}\leftarrow 2·\pi ·\mathrm{DT}20·\frac{f}{{\mathrm{fs}}_c}\\ {\varphi }_{\mathrm{adj}31}\leftarrow 2·\pi ·\mathrm{DT}31·\frac{f}{{\mathrm{fs}}_c}\\ {\theta }_{{02}_n}\leftarrow \mathrm{arrgh}\left(z02·e^{-j·{\varphi }_{\mathrm{adj}20}}\right)+{\varphi }_{\mathrm{adj}20}\\ {\theta }_{{13}_n}\leftarrow \mathrm{arrgh}\left(z13·e^{-j·{\varphi }_{\mathrm{adj}31}}\right)+{\varphi }_{\mathrm{adj}31}\\ m_n\leftarrow \sqrt{z02·z13}\end{array}\\ \mathrm{for}n\in 0\dots \frac{N_c}{2}-1\\ {\mathrm{Dvec}}_n\leftarrow m_n·\left({\theta }_{{02}_n}+j·{\theta }_{{13}_n}\right)\\ \mathrm{Dvec}\end{array}$ ${\mathrm{Dvec}}^{\langle k\rangle }:=\mathrm{getDvecspect}\left(X_{c0}^{\langle k\rangle },X_{c1}^{\langle k\rangle },X_{c2}^{\langle k\rangle },X_{c3}^{\langle k\rangle },{\mathrm{BulkDT}}_{c{20}_k},{\mathrm{BulkDT}}_{c{31}_k}\right)$

[0126] Compute weighted sum direction vector

[0127] ${\mathrm{Dvecsum}}_k:=\sum _{n=0}^{\frac{N_c}{2}-1}{\mathrm{Dvec}}_{n,k}$

[0128] Compute bearing error for each pulse

[0129] $\theta {\mathrm{err}}_{c_k}:=\arg \left({\mathrm{Dvecsum}}_k·e^{-j·{\theta }_k}\right)$

[0130] Now compare the three results:

[0131] $\left(\begin{array}{ccc}\mathrm{stdev}\left(\theta {\mathrm{err}}_a\right)& \mathrm{median}\left(\overrightarrow{\theta {\mathrm{err}}_a}\right)& \max \left(\overrightarrow{\theta {\mathrm{err}}_a}\right)\\ \mathrm{stdev}\left(\theta {\mathrm{err}}_b\right)& \mathrm{median}\left(\overrightarrow{\theta {\mathrm{err}}_b}\right)& \max \left(\overrightarrow{\theta {\mathrm{err}}_b}\right)\\ \mathrm{stdev}\left(\theta {\mathrm{err}}_c\right)& \mathrm{median}\left(\overrightarrow{\theta {\mathrm{err}}_c}\right)& \max \left(\overrightarrow{\theta {\mathrm{err}}_c}\right)\end{array}\right)·\frac{180}{\pi }=\left(\begin{array}{ccc}0.6928& 0.4886& 2.2567\\ 0.5096& 0.3562& 1.4211\\ 0.3638& 0.2388& 1.2215\end{array}\right)$

[0132] So typically method c wins, whether compared in terms of standard deviation, median, or maximum error. However, the advantage is not overwhelming, so the main advantage is that we can use a lower sample rate (and hence use less standby power).

[0133] Note that this analysis does not account for frequency-dependent bearing distortion (e.g. short multipath reception), and so there may be further, possibly more significant, advantages in bearing coherent processing when more complex environments are considered.

[0134] Although the invention has been described above with reference to one or more preferred embodiments, it will be appreciated that various changes or modifications may be made without departing from the scope of the invention as defined in the appended claims.