Acoustical Direction Finding with Time-Modulated Arrays

Time-Modulated Linear Arrays (TMLAs) offer useful efficiency savings over conventional phased arrays when applied in parameter estimation applications. The present paper considers the application of TMLAs to acoustic systems and proposes an algorithm for efficiently deriving the arrival angle of a signal. The proposed technique is applied in the frequency domain, where the signal and harmonic content is captured. Using a weighted average method on harmonic amplitudes and their respective main beam angles, it is possible to determine an estimate for the signal’s direction of arrival. The method is demonstrated and evaluated using results from both numerical and practical implementations and performance data is provided. The use of Micro-Electromechanical Systems (MEMS) sensors allows time-modulation techniques to be applied at ultrasonic frequencies. Theoretical predictions for an array of five isotropic elements with half-wavelength spacing and 1000 data samples suggest an accuracy of ±1∘ within an angular range of approximately ±50∘. In experiments of a 40 kHz five-element microphone array, a Direction of Arrival (DoA) estimation within ±2.5∘ of the target signal is readily achieved inside a ±45∘ range using a single switched input stage and a simple hardware setup.


Introduction
There are numerous situations where the range and origin of a transmitting source need to be determined. Some of the most common methods that achieve this use a receiver that electronically or physically scans space to find the Direction of Arrival (DoA). Electronic systems usually incorporate an array of sensors [1][2][3].
Once the input data is gathered it can be post-processed for DoA applications by established techniques such as the "MUltiple SIgnal Classification" (MUSIC) [4] and "Estimation of Signal Parameters by Rotational Invariance Techniques" (ESPRIT) [5] algorithms. These algorithms make use of sub-space methods which involve a significant amount of computation to produce an accurate direction estimation for an incoming signal; furthermore, the data of each element in the array usually needs to be recorded concurrently. Systems with the specifications to perform these calculations can become expensive and therefore there is a need to reduce complexity for smaller, low-cost devices. One such device may take the form of a fully integrated Micro-Electromechanical Systems (MEMS) sensor with multiple sensor elements, and a low-cost processing element.
Ultrasonic systems in air for applications such as robotic or vehicle guidance [6,7] may also use an array of transducers, and use techniques such as Time Difference of Arrival (TDoA) [8] or phase steering [9] to obtain a DoA estimation. For TDoA techniques, the use of cross-correlation requires each sensor to be sampled; for phase steered arrays, each element's phase needs to be controlled to produce a beam in a particular direction, which is usually swept across a range of directions.
The present paper considers electrically scanned systems operating in the receive mode only which utilise Time-modulated Linear Arrays (TMLAs). TMLAs function by sampling the array in a switched fashion rather than by means of phase steering [10] and benefit from only requiring a single analogue input stage with the sensors being switched on periodically to drive that input. The switching waveforms and frequencies can be varied in order to change the properties of the array.
The use of time modulation was originally conceived by Shanks and Bickmore [10] as a method of increasing the capabilities of transducer arrays. In a time-modulated array (Figure 1), elements can be periodically switched on and off. Due to the nature of the switching, power in the frequency spectrum is distributed amongst the signal and harmonic components related to the switching frequencies. This effect can be controlled by changing the switching waveform [11].
... Time modulated arrays have received increased interest recently in direction-finding applications due to their inherent benefits, and although the MUSIC algorithm has been adapted for use with time-modulated arrays [12,13] other researchers have found that the study of the harmonic components is all that is required for DoA estimation and as a result fast algorithms that have low computational complexity have been developed [14,15].
He [15] has demonstrated that it is possible to extract an estimation for a signal's DoA using a two element time-modulated array by performing a two-point Discrete Fourier Transform (DFT). Whilst this method can significantly reduce the complexity of the algorithms, it can also be prone to errors due to multipath reflections and noise, especially in situations where a small number of samples are captured.
The present paper suggests an alternative computationally efficient method using TMLAs to estimate the DoA and demonstrates the method's effectiveness when used in a five-element broadband microphone array to locate a 40 kHz signal in air. The use of this method with larger numbers of sensors shows that DoA estimations can be accurately produced even in the presence of noise and multi-path reflections.

Background
The normalized array factor of an N element, TMLA when receiving a signal from direction θ can be written as [16] AF(θ, t) where k is the wavenumber (rad/m), n is the element number, d is the spacing between elements in metres, and U n (t) is a time-switching function that turns each element either off or on in time defined by, where τ n,on and τ n,o f f are the switch on and off times respectively for the nth element in the array.
If each element is switched on then off uniformly in sequence, U n (t) becomes a periodic function and can be formulated as an infinite sum of m Fourier components each with coefficient c m,n [17] where the coefficients can be calculated as [11] c m,n = 1 T s The signal power is distributed among the harmonics of the switching pattern frequency f s = 1/T s , centred on a main frequency f 0 . Using the above equations, the array factor for a particular harmonic m can be calculated as, and if R is defined as the ratio of element spacing to signal wavelength λ, then the normalized array factor for each harmonic can be written as Figure 2 shows the normalized response of the dominant harmonics (i.e., f c ± 2 f m , f c ± f m , f c ) in a five-element array with each harmonic having a maximum amplitude at a different angle (i.e., ±53.1 • , ±23.6 • , 0 • ). From Equation (7), it can be shown that the direction in which the maximum response occurs for each harmonic can be calculated as: It can also be observed that in the case of a uniformly switched array (i.e., (τ n,o f f − τ n,on ) × N = T s , and each element is switched on in sequence), the point at which each harmonic's maximum amplitude occurs also marks the point where all other harmonics are at a minimum. The signal power associated with each frequency can be easily identified using a Discrete Fourier Transform (DFT). In practical terms this means that by analytical study of the amplitudes of each harmonic's DFT bin, a value for a given received signal's DoA can be derived. Array factors associated with each harmonic depend on the value of R as defined in Equation (6). As R becomes smaller, the main beam angles of each harmonic as calculated in Equation (8) will be more widely spaced. If R is greater than 0.5 (half-wavelength spacing), harmonics will obtain grating lobes due to the effect of spatial aliasing; in these cases, there are multiple angles in which harmonics become responsive, and the problem of direction-finding becomes more numerically complex. Throughout this paper, R is chosen to be 0.5 to ensure that each of the harmonics has a single main beam within the range of ±90 • .

Proposed Method of Determining the DoA
The approach taken in this paper is realised by finding the associated harmonic frequencies relating to both the signal frequency and the array switching frequency. For example, an array with five elements (N = 5) will distribute the power of an incoming signal into frequencies at f c − 2 f m , f c − f m , f c , f c + f m and f c + 2 f m . As shown in Equation (8), each of these harmonics has a known main beam angle θ m (i.e., −53.1 • , −23.6 • , 0 • , 23.6 • and 53.1 • ).
Consider a single frequency sinusoidal signal impinging on the array at an angle θ measured relative to the broadside direction. This signal will cause one or more frequency bins to become filled (as demonstrated by Figure 2) according to which harmonics are detected at the receiver. In the case where θ matches exactly a value of θ m , only one frequency is expected to be present in the output. It can then be assumed that if the value of the largest DFT bin is dominant, while the adjacent bins are small, then the signal angle is close to the main beam angle of that harmonic. Conversely, if the amplitudes of two adjacent bins are similar and also the largest measured, then the signal is approximately half way between the two main beams. The DoA can be estimated by looking at the average angle created by two beams, as weighted by the harmonic power level the array receives; estimation of the DoA can be expressed as: where X α and X α+1 are the two largest DFT bins that are adjacent to each other, and θ α and θ α+1 are the main beam angles associated with the corresponding harmonics.
The steps required to perform DoA estimation using the proposed method can be summarised as follows: 1. Collect a fixed number of samples from each microphone in turn (this can be repeated multiple times, however it is often adequate to do this just once). The harmonic frequency can be calculated as the reciprocal of the time taken to sample the set of all sensors. 2. Perform a DFT on the fundamental and expected harmonics of the switching frequency. 3. Select the DFT bins with the two largest amplitudes and identify their main beam angles. 4. Calculate a DoA using Equation (9) For example, if an unknown signal source produces the frequency spectrum shown in Figure 3 when sampled with a five-element array, then the power levels of harmonics 1 and 2 are taken as X α and X α+1 respectively and angles θ α and θ α+1 are found from Equation (8) Figure 4 shows the amplitudes of each harmonic plotted against the associated main beam angles. The source in this example is located at 25.0 • and the solid lines indicate the result for (N = 5). In this situation, the choices of X α and X α+1 are not clear. Harmonic 1 (m = 1) is the dominant frequency, but the two adjacent harmonics (m = 0, m = 2) are very similar in amplitude; both make possible candidates for the selection of X α or X α+1 . Using the method described, harmonics 0 and 1 will produce an estimation of 22.0 • whereas using harmonics 1 and 2 will produce an estimation of 25.0 • . Clearly, harmonic 2 is the correct value to use, however determining which of the two small harmonics can be challenging especially in the presence of noise. It can be concluded following this argument that the worst error case occurs when the incident angle of the signal is close to one of the harmonics' main beam angles.

Error Reduction Technique
This issue can be mitigated by processing the sampled data in a different way. If a number of samples has been taken from each element in turn but not repeated, it is possible to remove the samples recorded by the outer elements so that data from an N − 2 element array is obtained. This produces a different frequency spectrum due of the effect of removing elements, the modulation frequency becomes higher and the resultant harmonic beam directions change according to Equation (8). By selecting the number of elements used, errors caused by estimating the source direction using harmonic amplitudes close to the null regions can be avoided. Using the data provided, if harmonics 0 and 1 are taken as X α and X α+1 in the three-element case then the method produces an estimation of 24.6 • .

Numerical Results
In this section, results of a simulation developed in MATLAB R2015a are provided to illustrate the use of the proposed technique. A source of 40 kHz was assumed, impinging on a five-element array, where each element was separated by λ/2. A sampling rate of 500 kHz was chosen and 200 samples were taken from each simulated isotropic receiver in succession, meaning that the expected harmonics are separated by 500 Hz. The results are shown Figure 5   The abrupt changes in estimation error observed at approximately 23 • -25 • for the five-element array in Figure 5 are due to the algorithm selecting an incorrect smaller harmonic close to the null regions as explained in Section 2.2.1. In spite of there being no noise in this simulation, effects such as under-sampling and numerical errors still give rise to errors close to the null. If now we consider the estimation error for N = 3 (shown as the dotted line in Figure 5) and then the results are re-analysed using the same method with new values of harmonic data, the inaccuracies can be mitigated by taking a composite of the two traces. Since the value of N has changed, the null regions located using Equation (8) for N = 3 and N = 5 are non-coincident. It is known at which angles the inaccuracies will occur and hence by selecting the result of either N = 3 or N = 5, the best case accuracy overall can be achieved as shown in Figure 6. Multiple values of N were tested using a single source with additive white Gaussian distributed noise. The signal to noise ratio was fixed throughout all tests at 10 dB. Figure 7 shows that using a larger number of elements increases the accuracy of the method, as well as increasing the effective angular range of the system.  Using this method, it is also possible to detect multiple incoherent signals simultaneously if the frequency of the second signal and its harmonics generated by switching do not interfere with the DFT approximations performed on the first signal and its harmonics. In this case, the designer of the system must space the microphones at less than half a wavelength apart at the highest frequency to avoid spatial aliasing.

Experimental Results
A linear five-element microphone array using MEMS microphones was set up (Figure 8) with an inter-element spacing of 4.25 mm. This array spacing was chosen to be half the wavelength of a signal transmitted from a narrowband 40 kHz piezoelectric acoustic source. This source was a generic device of type "400ST160" having a bandwidth of approximately 2 kHz and a transmit pressure of approximately 110 dB Sound Pressure Level (SPL). It was driven from the 50 Ω output of a function generator producing a 40 kHz sine wave. A National Instruments "myRIO" was used for the processing of the received signal. The Field-Programmable Gate Array (FPGA) on the myRIO was setup to produce a digital waveform to control a set of complementary metal-oxide-semiconductor (CMOS) analogue switches that select between the individual microphone elements. A single analogue signal was captured at a rate of 500 kHz; 200 samples were captured for each element before switching to the next element meaning that 1000 samples were taken overall. The real-time processor chip on the myRIO performed DFTs at 39 kHz, 39.5 kHz, 40 kHz, 40.5 kHz and 41 kHz during the acquisition stage. Where null regions were detected, the process switched to using data from the centre three microphones, and DFTs were performed at 39.17 kHz, 40 kHz, and 40.83 kHz. The method described in Section 2.2 was used in LabVIEW 2015 to predict the source's one-dimensional heading. d= λ/2 The transmitter was securely mounted and fixed in position, 0.5 m from the array centre and carefully aligned so that the array was illuminated in the main beam. The microphone array was mounted on a positioner and rotated around its axis as shown in Figure 1. For each bearing angle measured, 1000 samples were taken and the DFT amplitudes were analysed to estimate the source direction θ est and this estimation was compared with the known direction θ. In addition, performance data was produced using LabVIEW's "Tick Count" facilities. Figure 9 shows that the beam patterns of the constructed unit are in good agreement with those predicted theoretically. It can also be seen that there are a few distortions in the patterns of the outer-harmonics; there are many possible explanations for this but the theory assumes omnidirectional receivers and does not account for the presence of a surface behind the array or the complexities of real transducers (e.g., intra-element coupling). Nonetheless, Table 1 shows that the results represent good DoA performance and agree with the overall trends that would be predicted. It is clear that the system is generally more accurate when using the data from a larger number of elements, but the use of switching to using the data from three elements has proved beneficial in the cases where the signal direction was close to one of the harmonic's null regions (23.6 • ).
The mean average time taken to convert the set of eight DFTs to an angle estimation was 29 µs. The time to perform the method described in [15] was also measured to compare the proposed method's computational complexity; the mean average time taken to convert the two DFTs to an estimation in this case was 41 µs. This is expected as the handling of complex arithmetic is usually significantly slower in these systems.  Figure 9. The received harmonic content of an experimental five-element linear microphone array compared to the expected numerical results.

Conclusions
A direction-finding method with both low hardware and computational requirements has been presented and introduced for acoustic signals. The paper has presented both numerical simulations and a full implementation of the algorithm using MEMS sensors to achieve fast direction-finding with ultrasonic frequencies; the use of small MEMS sensors makes practical implementations of ultrasonic acoustic arrays possible. The technique was demonstrated using a prototype system containing an FPGA and a real-time processor, but the method can be implemented using any single device with one analogue input and digital lines to switch between sensors in a circuit.
Furthermore, since the proposed method only requires real-valued magnitudes of the frequency content, it is possible to use optimised algorithms, such as the Goertzel algorithm [19], to compute the amplitudes of the harmonic content for the array, minimising the processing time and complexity needed to obtain a DoA estimation.
It is acknowledged that more computationally demanding estimators have been proposed and that in the case of [15], these make use of highly sensitive functions that may require careful handling of complex arithmetic. However, the technique demonstrated in this paper has obtained good DoA resolution in comparison and can be performed using relatively simple systems that can be adapted to suit a variety of requirements. If the system to be implemented requires high angular accuracy in the presence of noise, then the use of a greater number of microphones is preferable. If the system is required to find positions of signals that are short in duration, then it is possible to reduce the number of samples per microphone whilst keeping consistency in the direction-finding resolution.