AU5877099A

AU5877099A - Improvements relating to beamformers

Info

Publication number: AU5877099A
Application number: AU58770/99A
Authority: AU
Inventors: Keith John Jones
Original assignee: UK Secretary of State for Defence
Current assignee: Qinetiq Ltd
Priority date: 1998-09-15
Filing date: 1999-09-15
Publication date: 2000-04-03
Also published as: WO2000016308A1; EP1114412A1; GB9819973D0; JP2002525584A

Description

WO 00/16308 PCT/GB99/03075 IMPROVEMENTS RELATING TO BEAMFORMERS This invention relates to beamforming apparatus and methods of forming conventional beams for use, especially but not exclusively, with 5 an array of acoustic hydrophones. It is envisaged that the invention has application outside of the field of acoustic sensor arrays, for example for radar, communications, imaging, geographical exploration, astrophysical exploration and biomedical applications. 10 The object of a sonar system is to detect energy from one or more acoustic sources incident upon an array of sensors over a wide spatial field of view, whilst retaining the ability to reject signals from undesired directions in order to be able to determine the direction and range from which the incident energy is being received. 15 The requirements for a high sensitivity over a wide field of view with a high degree of directionality are to a degree mutually exclusive. In order to provide both a wide field of view and a good directional sensitivity, the outputs from the sensors in the array are processed using 20 an assembly known in the art as a conventional beamformer, referred to hereafter simply as a beamformer. Each sensor in the array typically has an omnidirectional sensitivity to incident radiation over a range of frequencies of interest. Each sensor 25 produces a time sampled output which is dependent upon the incident radiation and which is fed to the beamformer. Beamforming is then classically performed using a delay-and-sum process, whereby an appropriate time delay is applied to the sampled data from the sensors, and the time delayed samples then added together to produce an output 30 known as a conventional beam, referred to hereafter simply as a beam.

WO 00/16308 PCT/GB99/03075 2 Through the application of appropriate time delays to the sampled data, it is possible to steer the directional sensitivity of the array to ensure that acoustic wavefronts incident upon the array from the direction of the formed beam are sampled in-phase, as shown in Figure 1(a). When the 5 incident wavefronts are in-phase, the output sampled data is said to be coherent, and it is the coherence of the sampled data which produces the directional sensitivity. Typically, a multitude of different beams are formed covering a wide range of beam steering angles. Each beam is produced by applying different time delays to the sampled data. In one 10 possible system, the number of beams selected is chosen to ensure that a full field of view, between 00 and 1800 in azimuth is covered, with the difference in sensitivity over this range of angles dropping by no more than 3 dB, say, from the maximum value. This is shown in Figure 2 of the accompanying drawings. 15 Assuming that each sensor in the array produces a constant stream of discrete time sampled output data, with each output signal containing amplitude and phase information dependent upon the acoustic energy incident upon the sensor, the time-domain output of a conventional delay 20 and-sum beamformer can be written in terms of the sampled sensor data, x,(t), where m corresponds to the sensor address in the array, as: M b(t,e9, ,r) = wmxni(t -rm(O,#~,r)) where t is the time index, 0 is the beam steering angle in the azimuth plane, 0 is the beam steering angle in the elevation plane, r is the 25 focusing range and z-m is the time delay at the rn-th sensor for a given beam look direction and range. The term wm is a shading coefficient for the m-th sensor, used to control the spatial response or beam pattern for a given array.

WO 00/16308 PCT/GB99/03075 3 The value of the shading coefficient, win, will vary for each sensor in the array dependent upon its location within the array. It will also vary with the geometry of the array, and in a traditional beamformer its value is calculated individually for each sensor. For the beamforming of 5 complex-valued, band-shifted sensor data, the shading coefficients will also be steering-angle dependent. A change in the shape or size of the array during operation would require all the coefficients to be manually recalculated. This is costly and time consuming, especially where a hard wired beamformer is produced for which the coefficients cannot readily be 10 altered. A major problem in the past with the beamforming of a large array, especially where a large number of beams are required for good directional sensitivity, is that it is difficult to process the large number of 15 samples in real time. The arithmetic complexity of a conventional wideband beamformer is typically O(M 2 F) operations per second, where M is the number of sensors and F, is the sampling rate. For the addition of full narrowband operation, where spectrum analysis is applied to every wideband beam, this increases to: 20 O(M2 + M*10g 2 (K)]*F) where K is the length of the associated temporal discrete Fourier transform (DFT) algorithm and time-frequency transformations are assumed to be carried out using a fast Fourier transform (FFT) algorithm. 25 This problem is compounded when the data to be processed by the beamformer is in a digital form, because with the processing of such data, a malfunctioning sensor may default, with equal probability, to either the minimum or maximum level. For a system with a high dynamic range, a 30 single sensor faulting in the high state may completely corrupt the spatial WO 00/16308 PCT/GB99/03075 4 information produced by the beamformer to such an extent that no meaningful data can be obtained. This is especially true where the signal being measured is of a small amplitude. 5 In the past, it has not been possible to compensate for malfunctioning sensors in real time, and teams of analysts have been employed to "guess" the effect of a failed sensor on the beam pattern in order to explain anomalies in the formed beams after the measurements have been made. Typically, where a sensor array is fitted to a ship 10 during exercises, this means analysing the data once the ship has returned to port. The need for highly intensive processing of large volumes of data is exacerbated by the fact that the sampling rate from the sensors has, in 15 the past, needed to be made significantly higher than the minimum sampling rate suggested by Nyquists theorem. Such high sampling rates are required in order to facilitate accurate beam steering. By sampling at a shorter or finer time interval, the time delay applied to each sampled output in order to form the beams may be adjusted in finer steps at the 20 high sampling rate. Lower sampling rates approaching the Nyquist limit have, in the past, been used in combination with computationally-intensive finite impulse response (FIR) digital filters, which interpolate the sampled data 25 after sampling, and possess both unconditional stability and a desirable linear phase characteristic. Nevertheless, such an approach has required the provision of a single FIR filter for each channel of data to be processed. This approach is both costly and computationally intensive.

WO 00/16308 PCT/GB99/03075 5 It is therefore an object of the present invention to provide a method of beamforming or an apparatus for forming beams which ameliorates the problems of the prior art, especially in relation to beamforming arrays of sensors of arbitrary size and geometry and in 5 particular for large arrays. In accordance with a first aspect, the invention provides a beamformer for use in combination with an array of sensors, each sensor being adapted to produce a channel of output sequential digitised time samples in response to a measurand, comprising: 10 (a) one or more than one first primary combining units, each primary combining unit comprising a processing means adapted to process input channels of sampled data from a subset of one or more of the sensors, each primary combining unit being adapted to process the input channels of data to produce a plurality of intermediate output channels of 15 data, each intermediate output channel comprising sequential data samples corresponding to a partially formed beam for the sensors associated with that primary combining unit, and (b) at least one partial beam combining unit adapted to receive the intermediate output channels of data from the one or more primary 20 combining units, and further adapted to combine the intermediate channels of output data to produce one or more channels of secondary data, each channel of secondary data comprising a formed beam for the complete array of sensors. 25 The beamformer may be a wideband beamformer, where wideband refers to the fact that the temporal frequency range covered is from dc (i.e. zero frequency) up to the maximum temporal frequency of interest. Such a beamformer is also applicable to the processing of complex valued, band-shifted sensor data for which the real and imaginary WO 00/16308 PCT/GB99/03075 6 components of the beam data may be produced using the same delay-and sum beamformer routine. The array of sensors may comprise an array of pressure sensors, 5 such as acoustic hydrophones, whereby the measurand will comprise one or more acoustic waves incident upon the sensors of the array. Preferably, the one or more primary combining units are adapted to process the input channels to produce the intermediate output channels of data simultaneously in real time, i.e. on-line processing. They may 10 therefore operate using a parallel processing technique, and a parallel processor may be provided. Thus, the invention of the first aspect provides a beamformer in which the input channels of data from the array of sensors are processed 15 in parallel by dividing the sensors into subsets and processing each subset simultaneously using a primary combining unit before combining the intermediate partial beams to produce a set of formed beams. The beamformer may output one different output beam for each 20 look direction, or produce more than one formed output beam corresponding to each beam look angle, for example by producing several beams focused at different ranges in each direction. In a refinement, a second level of parallelism may be achieved by 25 providing a partial beam combining unit comprising: (a) a first intermediate combining means adapted to combine the intermediate output channels of data from more than one of the total number of primary combining units to form at least one first intermediate channel of data, WO 00/16308 PCT/GB99/03075 7 (b) a second intermediate combining means adapted to combine the intermediate output channels of data from a different set of more than one primary combining unit to form at least one second intermediate channel of data, and 5 (c) a final combining means adapted to combine the first and second intermediate output channels of data to produce the channels of secondary data. For example, for an array of 960 sensors, 16 primary combining units may be provided, each processing 60 input channels of sampled 10 data. The primary combining units may be combined in groups of 4 primary combining units to produce the intermediate formed beams using four intermediate combining means, a single final combining means being provided to combine the four sets of intermediate formed beams to produce the secondary channels of output data. In this example, each 15 primary combining unit is adapted to process the same number of input channels. Of course, individual primary combining units may process different numbers of input channels of sampled data. An advantage of the first aspect of the invention is that the parallel 20 processing of the input channels increases the speed at which the sampled data can be processed to produce the formed beams. This enables large arrays of sensors to be constructed and beamformed. Each individual primary combining unit may also, as a result, be less complicated than a single serial-design beamformer. 25 Each primary combining unit may comprise a means for storing in a memory the sampled data from the input channels associated with that primary combining unit. The input data samples may be stored in a space-time array in a memory. The space-time array may comprise m 30 rows of samples stored in a block of memory, where m is the number of WO 00/16308 PCT/GB99/03075 8 sensors associated with a primary combining unit and t is the temporal length of the stored sample data in each row. For example, sampling the input data at unitary intervals in time would result in t samples being stored in the array in each channel of temporal length t for each of the m 5 channels. The array may be updated at every sampling interval when a new set of m samples of data is produced by the sensors. For example, the sample data may be added to the left-hand end of each row in the store and the remaining data already in the store may ripple across to the right with the oldest set of m samples being ejected from the right-hand side of 10 the store. The temporal length t of the array may be chosen to exceed the maximum predicted time that an end-fire signal will take to cross the array of sensors, this time corresponding to the maximum time-delay 15 applied to the samples during beamforming in order to ensure coherent data samples. Of course, it will be readily appreciated that for an array of temporal length t, the number of samples stored in the array will depend upon the sampling rate of the data, a higher sampling rate resulting in more samples being obtained in a pre-set length of time t. 20 The primary combining unit may further comprise means for interpolating the data samples stored in the space-time store. The interpolation of the data enables the samples to be "upsampled" to a higher sampling rate. 25 The provision of an interpolation means is advantageous as it facilitates the high frequency sampling of the data required for accurate beamforming to be achieved from a low fundamental sampling rate. This minimises hardware requirements on the sensors, and where the sampled 30 data is stored as digital samples after having been processed using an WO 00/16308 PCT/GB99/03075 9 analogue-to-digital converter (ADC), the speed, and hence the cost, of the ADC is reduced. A low sampling rate at the sensors is also advantageous as it minimises the required bandwidth of the cables carrying samples from the sensors. 5 The interpolation factor or upsampling rate may be an integer value denoted by P. Thus, after upsampling the data, the space-time store will contain P times as many samples within the same temporal length t. Each primary combining unit may further comprise a beamforming 10 sub-unit adapted to produce the one or more partial beams using a delay and-sum process. The beamforming sub-unit may produce a number of different beams from the same set of data in the space-time store. This may be achieved by applying different time delays to the data samples. 15 In a most preferred arrangement, the time delays required to produce each different beam are stored in a memory table, each entry in the table being stored with a unique identifying address which associates the data at that address with a particular beam and row of samples in the space-time store. The table is preferably held in a fast access memory, 20 such as an on-chip memory. Thus, the beamforming sub unit may be adapted to look up the delay to be applied to the samples on each input channel using the appropriate address and, once the required delay values from the table 25 have been obtained, to use the delay values to select the required samples of data from the store. The selected samples may then be added together by the primary combining unit to produce the output sample value for the chosen partially formed beam.

WO 00/16308 PCT/GB99/03075 10 In effect, the look-up table may therefore tell the beamforming unit which of the samples (or interpolated samples) to select from the space time store to produce a beam. As an example, for a set of m rows of samples Xm(t), the look-up table data may tell the beamforming sub-unit to 5 combine the sample x,(t) from the first row, sample x 3 (t) from the second row, sample xs(t) from the third row... and so on, the m resultant samples being combined to produce a single output sample for one partial beam. This is repeated for each sampling interval to produce a channel of output samples. In this example sequence, the samples x 2 (t) from the first row, 10 x 4 (t) from the second row and x 6 (t) from the third row... and so on, would be added to produce the next output sample for the beam. After a preset number of samples on one output channel have been produced, a different set of time delays or addresses are then obtained 15 from the store and the process repeated to produce a second channel of data samples for a second beam. In a most notable arrangement, the interpolated input samples stored in the space-time store and the time delay information in the look 20 up table may be pre-loaded into fast local memory before performing the delay-and-sum beamforming process on a block of samples. Thus, all the samples required for use by the beamforming sub unit and the data telling the beamforming sub-unit which samples to select may be pre-loaded into a fast memory which can be quickly accessed by the beamforming sub 25 unit. The beamformed outputs produced by the beamforming sub unit may then be written directly to a second area of fast memory associated with the beamforming sub unit. This ensures that the computations required in forming the beams can be performed as quickly as possible and are not delayed whilst the relevant information is accessed. The fast WO 00/16308 PCT/GB99/03075 11 access memory is preferably on-chip memory associated with the processing unit which performs the delay-and-sum routine. The primary combining unit may further comprise means adapted 5 to apply a shading coefficient to each of the interpolated samples selected by the beamforming sub-unit when forming a set of samples to be added to produce a beam. The value of the shading coefficient required for each selected sample may be held in a second look-up table stored in fast memory, i.e. on-chip memory. 10 Preferably the primary combining unit is adapted to apply shading coefficients to the selected samples prior to the delay-and-sum process. As for the time delay look-up table, the shading coefficients may be stored in a table which can be accessed using the same address system as used to access the time delay values. Thus, the beamforming sub unit 15 may produce a single address corresponding to the number or type of beam to be formed, and based upon this single address the beamforming unit can access the correct set of information stored in the table. The beamformer may also include a sensor array configuration and 20 management (SACM) unit, comprising a processing means adapted to calculate the time delays or addresses required for one or more formed beams and the shading values to be applied for each sample. The SACM unit may be adapted to update the values stored in the look-up tables which are accessed by the primary combining unit. The SACM unit may 25 update the look-up tables automatically in response to one or more variable parameters or may be controlled manually, for instance by a programmer. By providing a number of identical primary combining units, and 30 providing a plurality of look-up tables containing the time delay and WO 00/16308 PCT/GB99/03075 12 shading information which are accessed automatically by the beamforming sub-unit within each primary combining unit, the SACM unit may set the data in the look-up tables to tailor the formed beams to any shape or size of array. For instance, by selecting the appropriate data for time delays 5 and shading coefficients the SACM unit may allow the output beams formed by the primary combining unit to correspond to the required beams from either a linear, cylindrical or planar array of sensors (or, in fact, for any array geometry). Because the primary combining units only need to access the information in the look-up tables and to receive 10 information about the sampling rate of the input data in order to select the samples for the delay-and-sum routine, they can therefore be "dumb" and have no specific features relating to the array of sensors which they are processing. Thus, each primary combining unit may be substantially identical. 15 The beamformer may be adapted for use in combination with means adapted to determine the relative position and/or geometry of one or more of the sensors in the array. The SACM unit may be adapted to update the time delay and/or shading coefficient values stored in the look-up tables to 20 compensate for changes in the array shape during use. This may be done in real time. The SACM unit may be adapted to calculate the relative positions of the blocks of sensors associated with each primary combining unit. 25 Thus, it may be assumed that the sensors associated with each primary combining unit are fixed relative to one another whilst the sensors from different primary combining units may move relative to the sensors of another primary combining unit. Such an assumption is feasible where the number of sensors associated with each primary combining unit is 30 small. For example, a small number of sensors is considered to be up to WO 00/16308 PCT/GB99/03075 13 approximately ten or twenty percent of the total number of sensors in the array. This assumption simplifies the problem of calculating the array shape and calculating the corresponding time delays and shading coefficients. This may be especially applicable to the processing 5 requirements of long towed arrays. Alternatively, the SACM unit may be adapted to calculate the positions of each sensor within the array, making no assumptions about the relative positions of the sensors within each primary combining unit. 10 A sensor data malfunction detection unit may be provided associated with the beamformer. The sensor malfunction detection unit may be adapted to detect malfunctioning samples and thus detect erroneous data samples from the sensors before they are processed by the 15 primary combining units. The detection unit may detect erroneous samples by analysing statistics of the samples stored in the space-time array. The unit may be adapted to analyse the samples of data in real time. 20 The detection unit may be adapted to produce a plurality of output error codes which may be stored in a look-up table. The one or more primary combining units may each be adapted to interrogate the error codes stored in the look-up table. Preferably, a single sensor malfunction detection unit is provided which is associated with a respective one of the 25 one or more primary combining units. Alternatively, a single sensor data malfunction detection unit may be provided which is associated with all of the primary combining units, or with a sub-group of the total number of primary combining units.

WO 00/16308 PCT/GB99/03075 14 The sensor malfunction detection unit may form a part of the SACM unit. The error codes may be stored in the same memory which contains the time delay and shading coefficients, which preferably is fast memory which can be accessed quickly by the primary combining units. 5 The term quickly means that the rate of access of the data by the primary combining unit is sufficiently fast that the data access rate does not form the dominant factor in determining the rate of performing the delay-and sum process. The memory may be on-chip with the microprocessor of the primary combining unit. 10 The sensor malfunction detection unit may be adapted to produce a number of different error codes, each error code corresponding to a different type of malfunction. For example, a first error code may be produced to indicate that a single malfunctioning sensor from a block of 15 interrogated input channels is present, a second different error code may be produced to indicate that a pair of consecutive malfunctioning sensors from the block have been detected, and a third error code may be produced to indicate that a block of three consecutive malfunctioning sensors has been detected. 20 A single error code may be produced for a number of samples in a row of the space-time store prior to interpolation. The block may comprise at least one sample or perhaps all the samples in a row of the store at any instant in time. 25 The primary combining units may be adapted to interrogate the error codes produced by the malfunctioning sensor detection unit and to replace data identified as being erroneous with estimated data samples or estimated interpolated data samples. The estimated correct data may be WO 00/16308 PCT/GB99/03075 15 calculated by analysing one or more of the samples surrounding the erroneous sample data in the space-time store. Preferably, the sensor malfunction detection unit analyses the 5 sampled data prior to the primary combining unit producing the interpolated data samples. Preferably, the beamformer sub-unit associated with each primary combining unit may be adapted to convert the time-domain sample data 10 (which may include estimated samples) in the space-time store into frequency-domain data prior to interpolation by the interpolation means. The interpolation may therefore be performed with frequency-domain techniques, and may use a suite of algorithms based upon an FFT routine for computational efficiency. 15 The output signal from the delay-and-sum process may be decimated to produce a beamformed output at a frequency that is different from the sampling rate of the samples from the sensors and different from the upsampled data rate. 20 The combination of decimation and upsampling is advantageous as it allows the formed beams to be output from the primary combining units (and hence ultimately from the beamformer) at a rate which could not be achieved by simply dividing a high initial sample rate by an integer 25 decimating factor, so that the beamformer may truly be regarded as a multi-rate digital signal processing system. SThus, it will be understood by the skilled man that the first aspect of the invention provides a beamformer which can be easily scaled for use 30 with different sized arrays by providing various combinations of primary WO 00/16308 PCT/GB99/03075 16 combining units which are each adapted to receive various numbers of input channels from an array. To form a complete set of output beams for the whole array, all the input channels should be processed by the primary combining units. However, beams can be produced using only a part of 5 the array of sensors. It also comprises means for the novel detection and compensation of malfunctioning sensors in real time, and provides a degree of parallelism in the processing of the input data to enable large arrays to be handled whilst producing a large number of formed beams in real time at relatively low processing rates. 10 The beamformer may be adapted to process either real-valued sensor data, or complex-valued sensor data. Real-valued data is typically encountered with a low-pass system such as encountered with passive sensor data. Complex-valued data is typically encountered with active sensor data. 15 The beamformer may be adapted to operate with base-banded complex valued sensor data. This is especially beneficial when the ratio of signal band width to the transmission or carrier frequency is small (i.e. typically less than 10%) as it allows the processing requirements to be reduced. 20 When the beamformer is adapted to process complex-valued sensor data, the data may first be basebanded (i.e. bandshifted, filtered and decimated). The shaded coefficients used by the beamformer may also be complex valued. 25 The beamformer may include two space time arrays, a first array being adapted to hold the in-phase sensor data, the second array being adapted to hold the quadrature sensor data. A single set of beam sie ing coefficients may then be provided along with two sets of beam s ing coefficient pairs for carrying out the combined operations of array WO 00/16308 PCT/GB99/03075 17 shading, liner interpolation and coherent summarisation of the data held in the space-time store. In one particularly noted arrangement, the beamformer may be 5 adapted to process both real valued and complex valued sensor data. Switching means may be provided whereby the mode of operation (real or complex) may be switched in real-time. The switch may be manually operated by an operator of the beamformer. 10 The first aspect of the invention may alternatively comprise a method of beamforming for use in combination with an array of sensors, each sensor being adapted to produce a channel of output sequential digitised time samples in response to a measurand, comprising the steps of: 15 (a) providing one or more than one first primary combining units, each primary combining unit comprising a processing means adapted to receive input channels of sampled data from a subset of the sensors, and in which each primary combining unit is associated with a different subset of the sensors in the array, each primary combining unit 20 being adapted to process the input channels of data to produce a plurality of intermediate output channels of data, each intermediate channel conprising sequential data samples corresponding to a partially formed beam for the sensors associated with that primary combining unit, and (b) providing at least one partial beam combining unit adapted to 25 receive the output channels of data from the one or more primary combining units, and adapted to combine the intermediate channels of output data to produce a plurality of channels of secondary data, each channel of secondary data comprising a formed beam for the complete array of sensors. 30 WO 00/16308 PCT/GB99/03075 18 In accordance with a second aspect, the invention provides a beamformer technique for use in combination with an array of sensors, each sensor being adapted to produce a channel of output data samples in response to a measurand, comprising the steps of: 5 (a) storing a block of data samples from mn sensors in a first store means; (b) analysing the data samples in the store means to produce statistics which are interpreted to determine whether the data from a sensor is erroneous; and 10 (c) providing a fault indicating means adapted to store an error code indicative of the function of each sensor in the array in a second store means. The sampled data may be stored in a memory and may be in the form of a space-time store comprising m rows of sampled data, each row 15 corresponding to a block of sampled data of temporal length t from a respective sensor. Step (b) may comprise the additional step of dividing the data stored in a row of the store into sub-blocks of data, and analysing each 20 sub-block to produce the statistics which are interpreted to determine whether the individual samples are erroneous. The fault indicating means may store the error codes in a sensor activity sequence, the activity sequence comprising a number of error 25 codes which are referenced by an address which corresponds to a unique sensor for a given block of data. The method may comprise a first or initial detection process, comprising the steps of measuring the median value of the integrated (i.e. 30 integration of the absolute sample values over time) temporal samples in a WO 00/16308 PCT/GB99/03075 19 sub-block, and comparing each individual integrated sample value in the sub-block to the median value to determine if the integrated samples and hence the original samples are erroneous. 5 This comparison is made possible by the fundamental fact that if at least fifty percent of the sensors are functioning, the median value of the integrated sample data will provide a reliable reference value from which comparisons may be made. In this way, an initial set of error codes corresponding to the first sub-block of data in each row may be produced 10 without the need for any historical data regarding the behaviour of the sensors. The method may comprise the further step of replacing samples of data which are indicated as erroneous by the fault detection means by estimated samples derived from the correct samples neighbouring the 15 faulty samples. For example, if the sensor activity sequence includes an error code indicating that one, two, or three malfunctioning adjacent samples (corresponding to malfunctioning adjacent sensors) are present, linear interpolation may be used to estimate the amplitude of the faulty samples, or extrapolation may be used to estimate the phase of the faulty 20 samples from correct samples on either side of the faulty samples. The beamformer may be adapted to form beams using the corrected data rather th-an the erroneous samples to produce more accurate results. Preferably, the method comprises the additional step of replacing a 25 block of sampled data by zero values where there are more than a predetermined maximum number of consecutive erroneous samples (corresponding to consecutive erroneous sensors), in this instance three, in a.sub-block.

WO 00/16308 PCT/GB99/03075 20 The method may comprise the additional step of updating the data samples contained in the store on each sampling interval by adding new samples to one end of the store and removing old samples from the other end of the store, so that the sampled data "ripples" across the store over 5 time. The method may comprise a second detection process used after an initial set of sensor activity sequence data has been derived using the first or initial detection process. 10 The second detection process, referred to as the standard detection process, may comprise, for each channel of sampled data, the steps of: (a) calculating the signal-to-noise ratio (SNR) for a number of sensors; 15 (b) removing one channel of data and calculating the SNR of the remaining sensors, after the channel has been removed; (c) comparing the two measured SNRs; and (d) determining whether or not a channel of samples is faulty based upon the results of the comparison. 20 The number of sensors for which the SNR is calculated may be all of the sensors in the array, or a sub-set of the total number of sensors. Preferably, the method comprises comparing the value of the SNR 25 for the removed channel with that of the remaining correctly functioning channels. In an alternative approach, instead of comparing the individual channel with the SNR of the remaining sensors in the array, the individual WO 00/16308 PCT/GB99/03075 21 channel may be compared with the SNR of the remaining correctly functioning sensor channels using a block of data from the store. The comparison of the data samples of one channel to the SNR 5 enables the functional status of a sensor to be determined because a malfunctioning sensor with a low SNR, for example, which falls below half the average of the remaining sensors, will effect a decrease in the SNR of the whole array. This decrease in SNR can be used to provide an indication that a sensor is faulty. Similarly, a malfunctioning sensor with 10 a high SNR might effect an increase in the SNR of the whole array. Preferably, the method may operate in real time to detect malfunctioning sensors and samples and calculate estimated data to replace the erroneous samples prior to input to a beamforming routine. It 15 may be automatic or semi-automatic. In an alternative approach, the second aspect may be considered to comprise a self-compensating beamformer apparatus for use in combination with an array of sensors, each sensor being adapted to 20 produce a channel of output data samples in response to a measurand, comprising: (a) a first store means adapted to store a block of data samples from m sensors; (b) a sample analysing means adapted to analyse the data 25 samples in the store means to produce statistics and interpreting means adapted to interpret the samples to determine whether the data from a sensor is erroneous; and (c) a fault indicating means adapted to store an error code indicative of the function of each sensor in the array in a second store 30 means.

WO 00/16308 PCT/GB99/03075 22 The self-compensating beamformer apparatus may further comprise a memory means adapted to define a space-time store, the store comprising m rows of the sampled data, each row corresponding to a block of 5 sampled data of temporal length t from a respective sensor. Means may be provided for dividing the data stored in a row of the store into sub-blocks of data, and analysing means adapted to analyse each sub-block to produce the statistics which are interpreted to determine 10 whether the individual samples are erroneous. In a further modification, the fault indicating means may define a sensor activity sequence in which the error codes are stored, the sensor activity sequence comprising a number of error codes which are 15 referenced by an address which corresponds to a unique channel block of data. It may further include measuring mean adapted to measure the median value of the integrated samples in a sub-block, and first 20 comparison means adapted to compare the individual integrated sample values in the sub-block to the median value to determine if the sample data in the sub-block is erroneous. The beamformer may further include means for generating 25 estimated samples from the correct samples neighbouring faulty samples, and means for replacing the faulty samples with the estimated samples. In another aspect the self-compensating beamformer may further comprise, associated with each channel: 30 WO 00/16308 PCT/GB99/03075 23 a) calculating means adapted to calculate the signal to noise ratio for the whole array of sensors; b) removal means adapted to remove one channel of data, the 5 calculating means being adapted to calculate the signal to move ratio of the remaining sensors; c) second comparison means adapted to compare the two measured signal to noise ratios; and 10 d) processing means adapted to determine whether or not a channel of samples is faulty based on the results of the comparison. In another arrangement the self-compensating beamformer may 15 further comprise a) an interpolating means adapted to interpolate the sampled data prior to beamforming; and 20 b) a decimating means adapted to apply an appropriate decimation factor to the output of the beamformer, in which the decimation factor is dependent upon the number-theoretic properties of the associated processing parameter so that the beamformed data is produced at an output rate which is different to the rate at which the data is sampled from 25 the sensors and different from the rate at which the data is interpolated. According to a third aspect, the invention provides a beamformer adapted to receive a number of channels of input samples from an array of sensors and adapted to produce one or more output beams by processing 30 the input samples, in which the beamformer comprises: WO 00/16308 PCT/GB99/03075 24 (a) a signal processing unit adapted to form a beam using a delay-and-sum routine to process the input samples; (b) a number of memory blocks; and (c) a data storing means, in which means are provided for 5 loading the input samples required to form a beam into a first block of memory prior to performing the delay-and-sum routine, with means being provided for storing the time delay values for each beam in a second block of memory. The delay-and-sum routine may be adapted to write the beam output 10 values to a third block of memory. Preferably, each block of memory comprises fast memory which can be quickly accessed by the delay-and-sum routine. Where the beamformer is implemented as one or more algorithms, which run on a 15 processor or chip, the memory may be on-chip memory such as cache memory. By providing sufficient fast memory for storing all the input time samples and time delay values prior to beamforming, and space for the 20 beamformed output values to be written in fast memory, the processing routine can operate faster than it would if it needed to wait to access data from a remote slow memory or wait for samples to be obtained and delivered to the processing routine. Thus, operation speed is optimised. 25 By "fast memory" it is meant any memory which can be quickly accessed by the processor operating the delay-and-sum routine. Preferably, the memory is built into the processor which performs the routine. Thus, the memory may comprise "on-chip" memory. It may be in the form of a memory cache. The access time of the memory is WO 00/16308 PCT/GB99/03075 25 preferably negligible compared to the time needed to perform the arithmetic processing needed to form the beams. The third aspect of the invention may be readily combined with the 5 first and second aspects of the invention to yield a high-speed parallel process. In such a case, the number of sensors allocated to each primary combining unit may be selected to optimise the use of the available fast memory on the processor. This may be dependent upon the amount of memory available, the number of sensors allocated to each primary 10 combining unit, and the temporal length of the store in which the sampled data is held. Alternatively, approaching the problem from the opposite direction, the amount of memory on the chip may be optimised to suit the number of sensors associated with each beamforming unit. 15 Where the problem of beamforming is broken down into several parallel operations using a number of primary combining units as described in relation to the first aspect of the invention, each primary combining unit would have the three blocks of memory to implement the invention of the third aspect. 20 The partial beams formed by the primary combining unit and stored in the third block of memory may be combined using a combination unit to produce formed beams for a whole array of sensors corresponding to a number of primary combining units. The resultant formed beams may 25 again be written directly into a fourth block of memory associated with the delay-and-sum process. Looked at another way, the third aspect of the invention may comprise a method of beamforming a number of channels of input samples WO 00/16308 PCT/GB99/03075 26 from an array of sensors to produce one or more output beams by processing the input samples, the method comprising: (a) providing a signal processing unit adapted to form a beam using a delay-and-sum routine to process the input samples, 5 (b) a number of memory blocks, (c) a data storing means, and (d) the additional steps of loading the input samples required to form a beam into a first block of memory prior to performing the delay and-sum routine, storing the time delay values for each beam in a second 10 block of memory prior to performing the delay-and-sum routine, and writing the beam output values to a third block of memory. According to a fourth aspect, the invention provides a beamformer for use in combination with an array of sensors, each sensor being adapted to produce a channel of output data samples in response to a 15 measurand comprising a processing unit adapted to process the channels of input data samples from a number of the sensors in the array to produce one or more formed beams, in which the beamformer further includes a SACM unit adapted to store information indicative of the time delays required for forming each beam in a store means, and in which the 20 processing unit is adapted to apply a delay-and-sum beamforming technique based upon the time delay values held in the store means. The SACM unit may be adapted to store information relating to the shading coefficients to be applied to each sample of data selected by the 25 processing unit during formation of the one or more beams. Thus, two independent units - a processing unit which forms the beams using a delay-and-sum routine and a SACM unit which stores information needed by the processing unit to form the beams - may be 30 provided which operate together to process the input samples.

WO 00/16308 PCT/GB99/03075 27 The SACM unit is so called since the values of the time delays and shading coefficients for the processing unit are dependent upon the configuration of the sensors connected to that beamformer from an entire 5 array. The SACM unit may further be adapted to store information relating to the operational status of one or more of the sensors, the processing unit being adapted to access the information relating to the 10 operational status from the SACM unit. For example, the SACM unit may be adapted to store one or more different fault codes, and a separate fault code may be stored which is associated with one or more samples obtained from the sensors. 15 Preferably, the SACM unit may be adapted to calculate the delay times to be used for each channel of data for a given beam based upon information about the type of array and/or the array size and/or spacing between the sensors and/or the focusing range and/or the data sampling rate. 20 The SACM unit may be user-programmable, and may include a graphical user interface (GUI) through which information about the shape, size, etc., of the sensor array may be input. For example, in order to construct a complete beamformer, a user may input data about the type of 25 array of sensors to the SACM unit. This may be through a graphical user interface. The user may also input a sampling rate corresponding to the rate at which samples are obtained from the sensors, and select a chosen output rate at which the formed beams will be output from the beamformer. The user may also input information about the spacing WO 00/16308 PCT/GB99/03075 28 between the sensors in the array, and for the case of a towed array, may input information about the position of the sensors. The SACM unit may be adapted to produce a plurality of tables of 5 time delays and/or shading coefficients and/or sensor status information, each table providing the information needed by a separate primary combining unit to produce a number of partial beams. Thus, a number of identical primary combining units may be operated in parallel. 10 The SACM unit further includes means to allow a user to input data to the SACM unit and/or change one or more parameters such as array size and sampling rate. The SACM unit may automatically recalculate the time delay and shading coefficient values if a change to the array parameters is made. 15 According to a fifth aspect, the invention provides a method of forming beams for use in combination with an array of sensors, each sensor being adapted to produce a channel of output data samples in response to a measurand in which each beam is produced as a stream of 20 beamformed output samples comprising the steps of interpolating the sampled data prior to beamforming and then applying an appropriate decimation factor to the output of the beamformer, in which the decimation factor is dependent upon the number-theoretic properties of the associated processing parameters so that beamformed data may be 25 produced at an output rate different to the rate at which the data is sampled from the sensors and different to the rate at which the data is interpolated. The method may therefore comprise the steps of forming the one or 30 more beams by interpolating the sampled data to an upsampled data rate, WO 00/16308 PCT/GB99/03075 29 producing an intermediate beamformed output by applying a delay-and sum routine to the interpolated data, and decimating the intermediate beam formed output to produce a beamformed output at the required rate. 5 The decimating factor may be selected to be exactly divisible (i.e. leaving no remainder) into the number resulting from the product of the interpolation factor and the number of valid samples in a temporal block of data to be beamformed. The interpolation factor may be chosen to be an odd number, whist the number of valid samples in the block may be 10 chosen to be an even number, so that the interpolation factor and the number of valid samples are relatively prime as required. The method may further comprise the step of linear interpolation of the interpolated data to further increase the accuracy of the beam output data. 15 Preferably, the decimation factor can be user adjusted to match the output frequency of the beamformer to the input frequency of any subsequent processing electronics for optimum data transfer rates. For example, when the beamformer comprises a number of primary combining 20 units operating in parallel to produce a number of partial beams, the data rate out of the primary combining units may be chosen to match the optimum data rate for the subsequent electronics used to combine the partial beams to produce complete beams. 25 The temporal interpolation may be carried out via frequency domain techniques using algoritms based upon the FFT routine. In an alternative approach, the fifth aspect may comprise a beamformer for use in combination with an array of sensors, each sensor 30 being adapted to produce a channel of output data samples in response to WO 00/16308 PCT/GB99/03075 30 a measurand, in which each beam is produced as a stream of beamformed output samples comprising: (a) an interpolating means adapted to interpolate the sampled data prior to beamforming; 5 (b) a decimating means adapted to apply an appropriate decimation factor to the output of the beamformer, in which the decimation factor is dependent upon the number-theoretic properties of the associated processing parameters so that the beamformed data may be produced at an output rate different to the rate at which the data is 10 sampled from the sensors and different to the rate at which the data is interpolated. According to a sixth aspect, the invention provides a method of beamforming for use in combination with a number of sensors, each sensor being adapted to produce a channel of data samples in response to 15 a measurand, the method comprising the steps of: (a) storing a block of sampled data form each of the sensors in a store means, whereby the temporal length of the store means is greater than or equal to the length of one block added to the time taken for an end fire signal incident upon the sensors to traverse the array of sensors, 20 (b) upsampling the block of sampled data in the store means and storing the upsampled data in the store means, (c) providing a second store means adapted to contain time delay values indicative of the sampled data from the store required to form more than one beam, 25 (d) producing a first block of data for a first formed beam for a block of stored data of temporal length equal to that of the temporal block t, based upon the time delay values of that beam, and (e) subsequently producing a second block of data for a second formed beam from the same block of stored data. 30 WO 00/16308 PCT/GB99/03075 31 In the prior art, it has been known to form a beamformed output data value for a first beam corresponding to a single moment in time and repeating until a set of beamformed output data for all the required beams has been produced before moving onto the next instant in time (i.e. the 5 next sampling period) and forming the next set of beamformed output data values. This necessitates changing the time delay values accessed by the beamformer m times (where m is the number of beams to be formed) for each instantaneous set of output samples before moving onto the next incident in time and forming a further set of beamformed output values. 10 The method of the sixth aspect of the invention provides that a complete temporal block of beamformed output values are produced for a single beam using one set of time delay data for all the samples in a block of samples (i.e. for every time sample in a block) before moving on to 15 apply the next set of time delays to the block of data to produce a second beam, and so on until all m beams are formed. In this way, in order to form n beams from m channels of sensor data, the beamformer need only access the time delay values m times. In the prior method, the beamformer must access different time delay data mxn times to form the 20 set of beams. The method may comprise providing n sets of time delay values stored under sequential addresses in an array, forming a block of data for a first beam using a first set of time delay values in the array 25 corresponding to a first address, and then incrementing the address used to access the time delay values to access the values needed to form the block of data for the next beam to be formed. Preferably, the sample data from the sensors may be stored in the 30 store means in the form of a space-time array having in rows of t samples, WO 00/16308 PCT/GB99/03075 32 where m is the number of input channels and t is the temporal length of the block. Each time sample in a row, Xm(t), may be provided with a unique address, the beamformer accessing the sample by selecting the appropriate address. 5 For example, for a block of data for each channel of length 256 samples, with an overlap of twenty-five percent between blocks and an output rate from the beamformer equal to the sampling rate of the sensor data, the beam steering addressees for the time delayed data samples need 10 only be accessed from the second store at a rate of once every 192 beamformed samples. In the prior art, the beam steering addresses would need to be accessed at the sampling rate. In an alternative approach, the invention may comprise a 15 beamforming apparatus for use in combination with a number of sensors, each sensor being adapted to produce a channel of data samples in response to a measurand comprising: (a) a store means adapted to store a block of sampled data from each of the sensors, whereby the temporal length of the store means is 20 greater than or equal to the length of one block added to the time taken for a signal incident upon the sensors to traverse the array of sensors, (b) an upsampling means adapted to upsample the block of sampled data in the store means and store the upsampled data in the store means, 25 (c) a second store means adapted to contain time delay values indicative of the sampled data from the store required to form more than one beam, (d) a means adapted to produce a first block of data for a first formed beam for a block of stored data of temporal length equal to that of 30 the temporal block t, based upon the time delay values of that beam, and WO 00/16308 PCT/GB99/03075 33 (e) a means adapted to subsequently produce a second block of data for a second formed beam from the same block of stored data. According to a seventh aspect, the invention provides a method of 5 interpolating sampled data from an array of sensors, each sensor being adapted to generate a channel of discrete samples in the time domain, comprising the steps of: (a) digitising the sampled data from the sensors for each channel, 10 (b) storing the sampled data in a store means comprising m rows of data having a first temporal length, each row comprising sequential data from a single channel, (c) dividing the data in each row of the store means into a number of blocks of data of a second shorter temporal length, 15 (d) converting the time-domain data in a block into the frequency domain, (e) applying an inverse DFT algorithm to a zero-padded version of the frequency-domain data, with the processing ignoring the zero valued data samples, to produce a block of interpolated time-domain data, 20 and (f) storing the interpolated data in the block in the store means, anid in which the blocks of data in the store means overlap in the time domain. 25 The store means may comprise a space-time store, comprising m rows of samples, where t is the temporal length of each row of the space time store and m is the number of sensors producing data, each row being divided into a number of blocks which overlap and therefore share common samples at the end. Thus, each sensor produces a stream of time 30 samples which fill each row of the store.

WO 00/16308 PCT/GB99/03075 34 Once the rows of the store have been filled, the time samples may ripple across the store with new samples entering at one end and the old samples being ejected from the other end, interpolation being formed on a 5 complete new block of data only once a complete data block of new samples has been added to the store. The overlapping of the data blocks, with the interpolation being performed on a single data block at a time, is advantageous, as most of 10 the errors introduced by a DFT-based technique for band-limited temporal interpolation arise from processing the ends of a block of data. The blocks may overlap by any amount, but an overlap of about twenty-five percent is preferred, an equal amount of overlap being provided at each end of the data block. 15 The blocks may comprise, in one arrangement, 256 data samples. A twenty-five percent overlap may be provided, which leaves 192 samples for use by a DFT algorithm to perform the interpolation. 20 The block of data may be interpolated using a 32-bit floating-point arithmetic processor running a suite of FFT algorithms to perform the DFT-based interpolation routine. The block of data may be multiplied by an appropriate window 25 function such as Hanning, Triangular, Hamming or Kaiser window to further increase the accuracy of the interpolation, reducing Gibbs-type effects. After the interpolated data has been returned from the frequency domain to the time domain, it may be scaled using a simple function to remove the effect of such a window. 30 WO 00/16308 PCT/GB99/03075 35 The frequency-domain interpolation algorithms may comprise a suite of compact FFT algorithms. The interpolation may employ a dual-channel processing scheme 5 whereby two channels of sensor data are processed simultaneously to increase the computational efficiency of the processing algorithms. The sampled data may be interpolated by zero-padding of the frequency-domain data and applying an inverse DFT, or IDFT routine. 10 The FFT algorithm suite may comprise an inverse Sino spectral interpolation algorithm based upon a row-column implementation of the IDFT process. Thus, the seventh aspect provides a method of interpolating sample 15 data from an array prior to beamforming using high accuracy frequency domain techniques. Looked at another way, the seventh aspect of the invention may comprise apparatus adapted to interpolate sampled data from an array of 20 sensors, each sensor being adapted to generate a channel of discrete samples in the time domain, comprising: (a) a digitising means adapted to digitise the sampled data from the sensors for each channel, (b) a store means adapted to store m rows of sampled data 25 having a first temporal length, each row comprising sequential data from a single channel, in which the stored data in each row of the store means into a number of blocks of data of a second shorter temporal length, (c) a time-domain to frequency-domain conversion means adapted to convert the time-domain data in a block into the frequency 30 domain, and WO 00/16308 PCT/GB99/03075 36 (d) a processing means adapted to apply an inverse DFT algorithm to a zero-padded version of the frequency-domain data, with the processing ignoring the zero-valued data samples, to produce a block of interpolated time-domain data, and store the interpolated data in the block 5 in the store means, in which the blocks of data in the store means overlap in the time domain. There will now be described, by way of example only, an embodiment of the present invention with reference to the accompanying 10 drawings of which: Figure 1 illustrates the principal of beamforming by applying time delay values for (a) a far-field beam incident at an acute angle to the array, and (b) a beam from a near-field source; 15 Figure 2 illustrates an ideal set of formed beams covering a wide range of beam look angles; Figure 3 is a schematic overview of a complete acoustic array system using a beamformer embodying several aspects of the invention; 20 Figure 4 illustrates the relationship between the functional building blacks of a complete beamformer for processing 960 acoustic sensors using 16 primary combining units; 25 Figure 5 is an enlarged view of the functional relationship between a block of four primary combining units, its associated intermediate combining unit and the SACM unit; WO 00/16308 PCT/GB99/03075 37 Figure 6 is an enlarged view of the interrelationship between the four intermediate combining units and the final combining unit which produces the final beamformed output data; 5 Figure 7 is a flow diagram illustrating the passage of the sampled data from the sensors through a primary combining unit and the role of the space-time store as a storage block for a time history of sampled data accessible to the primary combining unit; 10 Figure 8(a) shows the passage of the sampled data through the primary combining unit to produce interpolated data prior to application of a delay-and-sum process to produce the partial beam data; Figure 8(b) shows the subsequent processing of the interpolated data within the beamformer using time delay and shading coefficient values 15 produced by the SACM unit; Figure 9 illustrates the overlap of blocks of data for processing by the interpolation sub-unit of the primary combining unit; 20 Figure 10 shows the application of a Hanning window to the data in a block prior to interpolation; Figure 11 shows the arrangement of sampled data within the space-time array and the subdivision of rows of data into sub-blocks for processing 25 by the malfunctioning sensor detection unit; Figure 12 shows a typical theoretical beamformed output obtained from a prototype beamformer after compensation for five malfunctioning sensors which have failed in the "low" state; 30 WO 00/16308 PCT/GB99/03075 38 Figure 13 shows a similar set of traces to Figure 12 obtained where five sensors which have failed in the "high" state are included in the array; Figure 14 is a block diagram showing the row-column implementation of 5 the prime factor inverse DFT algorithm; Figure 15 is a block diagram showing the modified row-column implementation of the prime factor inverse DFT algorithm; 10 Figure 16 is a flowgraph of the inverse Sino spectral interpolation algorithm; Figure 17 illustrates a single dual-node processing unit; and 15 Figure 18 illustrates the Sino computing engine embodied in the beamformer. A complete beamformer has been implemented which employs many of the features of the present invention. An overview schematic of 20 the complete beamformer is shown in Figure 3 of the accompanying drawings. It comprises an array of m acoustic sensors 1 which each produce a channel of time sampled data. The m channels of data are fed into a beamformer unit 2 (described in detail hereinafter) which produces a number of beamformed outputs that are fed to a display 3 for 25 interpretation by an operator. Whilst the present invention is applicable to many types and shapes of sensor arrays, the present example is described in relation to a linear array of acoustic sensors or hydrophones for detecting incident acoustic 30 signals when towed behind a tow vessel. Of course, it will be understood WO 00/16308 PCT/GB99/03075 39 that the invention is not to be limited to towed arrays or linear arrays or acoustic sensors and many applications in other fields are envisaged. Conventional beamforming with a towed array comprising a line of 5 hydrophones has traditionally been carried out by assuming the array to be linear and horizontal, and then to produce beams by applying processing techniques to the sampled data obtained from the sensors in the array such as the well known delay-and-sum approach, to yield a computationally acceptable set of formed beams. 10 For a large towed linear array, of the order of 1000 metres long, the assumption that the array is of fixed uniform shape is not valid. For example, manoeuvres of the tow vessel can cause large deformations of the array shape, producing erroneous beamformed outputs unless 15 knowledge about the changing array shape is built into the processing. A further assumption that has been made in the past is that the wave fronts impinging upon the array are planar, i.e. in the far field. However, for a sufficiently large array of sensors, with the sensors spaced 20 by a minimum distance of 2/2 (to avoid spatial aliasing), where 2 is the wavelength corresponding to the highest frequency of interest, and assuming an array of 960 sensors, the total length of the array will be 806.4 metres. In this case, the near field extends to approximately 200 kilometres. As such, the assumption that all targets will be in the far 25 field will be clearly flawed. The present invention improves on prior art sensing arrays by providing means for focusing the formed beams so that objects in the near field may then be detected and range information produced. 30 WO 00/16308 PCT/GB99/03075 40 A further constraint on the prior art design of beamformers is the effect of malfunctioning sensors on the formed beams. This is especially significant in digital systems, where the output from a sensor may fail in either the high or low state. 5 The beamformer developed and described herein provides a solution to the problem of beamforming a large array of sensors, whilst compensating for changes in array size, attitude, the effect of malfunctioning hydrophones, and providing the possibility of focusing the 10 formed beams to within the near field of the array. The system employs digital design techniques and a degree of parallelism in the formation of beams to produce a fast compact system which can be tailored for use with varying sizes and shapes of array. 15 Specifically, the beamformer developed and described herein by way of example is adapted to form beams from an array of 960 hydrophones arranged in the shape of a simple linear array. With a design frequency of 880 hz, a wavelength of approximately 1.68 metres, 20 this results in a total array aperture D of 806.4 metres, the maximum duration for a signal travelling across the array being therefore 0.55 seconds. The processing system has a maximum temporal frequency of interest f, of 933 hz, and data from the sensors is sampled at a rate of 2800 hz, i.e. only 1.5 times the Nyquist rate in order to meet the anti 25 aliasing requirements of the system. The full array therefore produces 2,688,000 samples of data every second. Such a large amount of data obviously places a strain on any processing system, and this problem is addressed using parallel processing techniques.

WO 00/16308 PCT/GB99/03075 41 An 18-bit digital signal processing (DSP) system providing 18 bits of dynamic range is used, offering a maximum dynamic range of 102 decibels by storing data in a 24-bit 2s-complement form. 5 In order to reduce considerably the amount of processing needed so as to cope with the large amount of sampled data obtained from the hydrophones, the beamformer comprises a number of identical sub-units 4 which each process a proportion of the total number of hydrophones. From 960 hydrophone channels of input data, each sub-unit 4 processes 10 240 channels of hydrophone data. Thus, four sub-units 4 are needed as shown in Figure 4 and Figure 5. The four sub-units produce partial beam outputs which are combined by a combined processing unit 5 to produce the beamformed outputs. To further enhance the processing speed of the beamformer, each 15 sub-unit is divided into four primary combining units 6 and a sub combining unit 7, each primary combining unit 6 receiving the sampled input data from 60 hydrophones in the array. As shown, each primary combining unit 6 receives inputs from a block of adjacent hydrophones in the array and, using a suite of algorithms, processes the incoming data. 20 However, the hydrophones may in fact be divided randomly between the primary combining units 6 and need not be entirely adjacent. The dual level of parallelism provided by the primary combining units which feed into the sub-combining units and then the combining 25 units means that many of the arithmetic data processing steps involved in forming a beam from the input samples can be performed simultaneously, thus increasing processing speed. Increased speed allows larger arrays and higher numbers of beams to be handled for a given processor speed.

WO 00/16308 PCT/GB99/03075 42 As mentioned hereinbefore, each sub-unit 4 is identical, and each primary combining unit 6 is also identical. Beamforming is performed using a delay-and-sum technique within each primary combining unit, a set of partial beams from each primary combining unit being combined by 5 the sub-units to produce an intermediate set of partial formed beams which in turn are combined in a final sub-aperture combination unit to produce a single set of formed beams for the whole array of sensors. In addition to the primary combining units and the sub-combining 10 units, a SACM unit 8 is provided. The SACM unit 8 controls the information needed by each primary combining unit 6 to perform the correct delay-and-sum-routine. Each primary combining unit 6 obtains the required time delay values by reading values stored in a block of memory, and the time delayed values assigned to each primary combining 15 unit 6 will be different from primary combining unit 6 to primary combining unit 6 depending on the position of the hydrophones associated with that primary combining unit 6 within the whole array 1. It will thus be readily apparent to the skilled man that the primary combining units 6 are essentially dumb, performing beamforming based upon time-delay 20 values assigned to the memory 9 by the SACM unit 8. The SACM unit 8 also controls the values of the shading coefficients required by the delay an-d-sum beamforming process performed by each primary combining unit 6. The role of the shading coefficients is to reduce side lobes in the formed beams, and their values are written to the memory 9 by the SACM 25 unit 8 which is then accessed by each primary combining unit whilst forming a set of beams. A further role served by the SACM unit 8 is to analyse the samples generated by the sensors prior to beamforming, in order to assess the 30 accuracy and reliability of each of the samples. This is performed by WO 00/16308 PCT/GB99/03075 43 analysing statistics of the sampled data over time. The SACM unit 8 stores a number of fault codes corresponding to each sample or set of samples in the form of a sensor activity sequence 10, and the primary combining units 6 are adapted to access the fault codes from the SACM 5 unit 8 prior to beamforming. The primary combining units 6 can then alter the sampled data to compensate for erroneous samples. Thus, in its simplest form, the invention enables an array of virtually any size or shape to be beamformed using a system built from a 10 combination of just three simple building blocks: (a) a SACM unit 8; (b) a primary combining unit 6; and (c) a combining unit 5 and/or sub-combining unit 7 which reconstructs a complete set of beams from a set of partial beams formed 15 by each primary combining unit 6. This arrangement provides several degrees of flexibility. Firstly, the number of input channels of data connected to each primary combining unit (i.e. the number of sensors processed by each unit) can be tailored to suit the hardware and software upon which the beamformer is created. In 20 particular, the size of the primary combining unit can be tailored to suit the amount of available memory on a processor board, or to suit the desired sampling rate for a range of scanned frequencies using a hardware/software platform of fixed clock speed. 25 For a given size of array, increasing the number of primary combining units increases the amount of parallel processing and generally speeds up the rate at which beams can be formed. Thus, higher sampling rates or a larger frequency range of interest can be accommodated.

WO 00/16308 PCT/GB99/03075 44 The primary combining units may be identical which reduces the amount of design involved in producing a complete beamformer. Only a single primary combining unit need be designed which accesses a number of channels of input data, accesses time delays for beamforming from a 5 memory, and outputs a number of partially formed beams. The time delays are written to the memory by the SACM unit 8, and thus the behaviour of the beamformer and the type and number of beams to be formed from the input data can be altered by reprogramming the SACM unit 8 without altering the design of the primary combining unit. 10 By sharing a single SACM unit between all the primary combining units in the beamformer, a single set of changes to the SACM unit may produce all the changes required for each primary combining unit. For example, the SACM unit may be programmable through a graphical user 15 interface, allowing an operator to enter data about the number and type of beams to be formed and the array shape and size. A description of the operation of the various units of the beamformer is given in more detail hereinbelow. 20 The Space-Time Store Each sensor produces a channel of output data samples which contain phase and amplitude information about the acoustic signals incident upon the sensors. The sampled data from the sensors is held in a 25 space-time store 10 in memory which can be quickly accessed by a primary combining unit. An individual space-time store 10 of data is established for each primary combining unit. The store comprises m rows of sampled data, each row being of 30 temporal length t, where m is the number of sensors associated with the WO 00/16308 PCT/GB99/03075 45 primary combining unit. In this case, m = 60. The data in each row of the space-time store is divided into blocks of data of temporal length tb as shown in Figure 9. A number of blocks are formed, with the samples in each block overlapping by twenty-five percent. 5 The blocks are selected to have a length of 256 samples, with a twenty-five percent overlap. The purpose of dividing the data into blocks is two-fold. Firstly, the primary combining unit interpolates the sampled data block by block. Secondly, beamformed output data for each channel 10 is produced one block at a time. Overlapping the blocks increases the accuracy of the frequency-domain interpolation techniques employed. The size of the store is carefully selected to match the amount of memory available to the primary combining unit. The temporal length of 15 each row in the space-time store is chosen to be greater than or equal to the length of a temporal block tb, plus the maximum time that an acoustic signal will take to pass from one end of the array of sensors to the other. This ensures that when the primary combining unit forms the beams for a block of time samples, all the necessary delayed time samples are present 20 in the space-time array. Design Philosophy of the Primary Combining Unit. The role of the primary combining unit 6 (which can be 25 implemented in hardware or software) is to process a number of channels of sampled sensor data ( i.e. rows of data held in the space-time array) to produce a set of partially formed beams. Because the primary combining units operate in parallel, an increase in the number of primary combining units for a given size of array increases the degree of parallelism and 30 hence the processing throughput.

WO 00/16308 PCT/GB99/03075 46 To form accurate beams via conventional delay-and-sum beam forming techniques, it is necessary that the effects of time-delay quantisation are kept to a minimum. This is most easily achieved by using 5 sampled sensor data which has been sampled at a sufficiently high rate. Prior art beamformers have used sampling rates of 10 to 15 times that of the Nyquist rate, although for exact beam formation this should be greater than or equal to N, where N is the number of independent beams to be formed. Thus, for an array of 960 sensors, the sampling rate should 10 ideally be 960xNyquist. Such high initial sampling rates, however, are undesirable as they place a severe burden upon the ADC stage and upon the transmission bandwidth of the cables used. To avoid the problems presented by high initial sampling rates, the 15 required sampling rate for the beamforming sub-unit is obtained by upsampling the sampled sensor data within the primary combining unit. Thus, a low initial sampling rate can be used which need only meet the anti-aliasing requirements. For example, with a maximum frequency of interest of f,=933 Hz,, and a sampling rate of F=2800 Hz 20 (i.e.1.5xNyquist), the data can be upsampled by a factor of 7 in order to provide interpolated data to the beamformer corresponding to an acceptable rate of 10.5xNyquist. Such a sampling rate is generally considered sufficient to maintain the on-axis beamformer response to within an acceptable margin of error of the theoretical level. 25 As will be appreciated, the effect of interpolating the sampled data prior to beamforming is to facilitate increased accuracy for beam steering. Thus, associated with each partial beam is a set of integer time delay values, k,, for i= 1...m (where m is the number of sensors associated with 30 a primary combining unit). The primary combining unit accesses the WO 00/16308 PCT/GB99/03075 47 values of k, from the SACM unit, the values of k, being held in a block of fast access memory which may be quickly accessed by the primary combining unit. It is therefore relatively trivial to change the behaviour of the beamformer by varying the values of k, held in memory using the 5 SACM unit. A further benefit of upsampling the sampled data prior to beamforming is that it enables a lower initial sampling rate to be employed. For example, in a typical existing beamforming system, the 10 input data is sampled at a rate 10 to 15 times higher than the Nyquist rate, although for exact beam formation this figure should be greater than or equal to N, where N is the number of independent beams to be formed. Such a high sampling rate places a great strain on the hardware associated with the sensors, especially in terms of ADC and the transmission cable 15 bandwidth. In contrast to the requirements for a high input sampling rate, the output of the beamformer need only be produced at a rate in line with the Nyquist criterion. Thus, the output of the beamforming sub-unit of the 20 primary combining unit is decimated to a lower rate than the upsampled rate by dividing by an integer factor to yield an output rate which may not necessarily be a factor of the initial sampling frequency. Consider, for example, that the incoming hydrophone data, at a 25 sampling rate F, is upsampled by a rate P and then decimated by a factor Q, giving an output frequency for the partially formed beams of: Output frequency = (P x F) / Q where P and Q are integers.

WO 00/16308 PCT/GB99/03075 48 It is therefore possible to produce an output frequency which more closely matches the frequency requirements of the next stage of the beamforming system, in this case the combining units. 5 After interpolation, beamforming is performed by the primary combining unit implementing a delay-and-sum technique using the up sampled data. In the first stage of this technique, the primary combining unit applies a respective time delay k, to the samples on each channel stored in the memory. Different time delays are required for each beam 10 look direction. The interpolation is performed using a suite of computationally-efficient, frequency-domain algorithms. In the second stage of the beamforming process, after interpolation and selection of the required samples, each sample selected to form a 15 beam is multiplied by a shading coefficient, w,, and then added together. The effect of the shading coefficient is to reduce the magnitude of the unwanted sidelobes of the formed beam. For a dynamic array, which may be changing shape during use, a different set of shading coefficients will be required for each beam direction. Also, for the beamforming of 20 complex-valued, band-shifted sensor data, a different set of shading coefficients will be required for each beam direction. The shading coefficients are stored in memory which can be accessed by the primary combining unit when forming beams. As for the time delays, the values of the shading coefficients, w;, are adjustable using the SACM unit. 25 Thus, identical primary combining units can be used to build up a wide variety of beamformers for different array types and sampling rates, etc. Overview of Operation of the Primary Combining Unit.

WO 00/16308 PCT/GB99/03075 49 An overview schematic of the operation of the primary combining unit as a part of the complete beamformer is shown in Figure 7. A more detailed schematic of the primary combining unit is shown in Figures 8(a) and 8(b). 5 Referring to Figure 7, data sampled from the hydrophones is first digitised by an ADC 11 before being stored in the space-time store 10 which can be accessed by the primary combining unit. The sampled data is then processed by an optional malfunctioning hydrophone detection unit 10 12 (described hereinafter) before being passed directly to the primary combining unit 6 for processing. The malfunctioning hydrophone detection unit 12 produces a number of error codes which are stored by the SACM unit 8 in memory 9 which can be accessed by the primary combining unit. The primary combining unit comprises an interpolation 15 sub-unit which interpolates the sample data in the space-time store, a malfunctioning sensor compensation unit which replaces faulty samples with estimated data, and a delay-and-sum sub-unit which selects the required samples for producing a number of blocks of beamformed output data. 20 As shown in Figure 8(a), the primary combining unit 6 processes the data in the store in 256 sample overlapping blocks. Each set of blocks 11 is windowed and the data at the end of each block is ignored. The blocks of sample data are then converted to the frequency domain using a 25 digital signal processor adapted to perform an FFT routine. The frequency-domain data produced, which may be left in scrambled (i.e. bit reversed) form, is then passed to a band-limited interpolation sub-unit 14 which forms a part of the primary combining unit. This applies a compensating algorithm to the data to compensate for erroneous samples 30 before zero-padding the compensated spectrum and applying an inverse WO 00/16308 PCT/GB99/03075 50 DFT algorithm using the FFT routine. This produces a set of blocks of interpolated and corrected data which is written back to the space-time array 10. A scaling factor is applied to remove the effects of the window. This is repeated until all the blocks have been interpolated. The routine is 5 next performed when a fresh set of 256 samples has been written to the space-time store. As can be seen in Figure 8(b), the primary combining unit further comprises a delay-and-sum sub-unit 15 adapted to access the time-delay 10 values held in the space-time store 10 to obtain all the samples 16 needed to form a first beam for a single time instant. Shading coefficients are also applied at this stage using a shading coefficient multiplication unit 17. The samples are then added 18 to 15 produce a single data output value 19 for the first of the partial beams formed by the primary combining unit. This is repeated for each time interval in a block until a block of output values, of length 192 samples, has been produced for the first beam 20. 20 The delay-and-sum sub-unit then accesses the data for a second beam and generates a second block 20 of samples 19 for that beam and so on until a set of blocks of samples 21 has been produced for all the partial beams. The process then restarts for the next temporal block of 256 samples in the space-time store. In practice, as data ripples across the 25 array, the unit effectively processes the same (i.e. the latest) block of 256 samples once they have been interpolated. Detailed Description of the Interpolation Sub-Unit WO 00/16308 PCT/GB99/03075 51 Interpolation occurs prior to beamforming using the beamforming sub-unit. An approach has been developed for generating the upsampled data from the sampled hydrophone data which relies upon converting the sampled data into digital time-domain data and then converting into the 5 frequency domain using a FFT routine, zero-padding the spectrum, and then returning to the time domain with an inverse Sino spectral interpolation algorithm. The compensation can be built into the interpolation algorithm, producing a fast and computationally-efficient solution. 10 Inverse Sino Spectral Interpolation Algorithm To interpolate the sampled data, frequency-domain techniques may be adopted using a suite of algorithms based upon the FFT algorithm. A 15 suite of such algorithms has been developed which is computationally more efficient than the application of an FIR filter for the 960 sensor, 16 primary combining unit array. The suite of algorithms is based on the use of the DFT algorithm 20 which is given in 1-D form by the expression N-1 X[k] = Z x[n]Wnk n=0o where WN = exp{-i2r / N} the primitive N-th complex root of unity, with the inverse DFT, or IDFT, 25 being given by 1 N-I x[n]= - X[kW nk N k=o WO 00/16308 PCT/GB99/03075 52 The DFT-based temporal interpolation technique involves first transforming blocks of time-domain data to the frequency domain, via a DFT, zero-padding the centre of the resulting spectrum, then returning to the time domain via an IDFT. The DFT and IDFT routines may be 5 efficiently carried out with FFTs, whilst the amount of zero-padding determines the order of interpolation achieved. Note that the zero-valued samples are placed in the centre of the spectrum in order to preserve the symmetry of the spectral waveform. The main limitations of such a scheme are that conventional implementations suffer from high arithmetic 10 complexity, arising from the size of the IDFT, and low interpolation accuracy, arising from the presence of Gibbs-type effects in both the time and frequency domains. There is thus a need, in the high-performance DSP field, for a solution which is able to overcome both of these limitations. 15 In the suite of proposed algorithms used in the beamformer, the IDFT is modified to produce interpolated time-domain samples between the original N time-domain samples by means of a process referred to as zero-padding. Suppose that N zero-valued samples are inserted into the 20 centre of the DFT spectrum, denoted by X[O],X[1], ... , X[N-1], yielding a new spectrum Y[O],Y[1],..., Y[2N-1], i.e. an interpolating factor of two is to-be used. Then the IDFT of the new spectral sequence is given by 1 2N-1 x[n] = Y[kk2 2N k=O which reduces to 1 N-1 25 x[n] = X[k Nk=o where the range of summation has now been adjusted to reflect the fact that the middle N samples are zero-valued.

WO 00/16308 PCT/GB99/03075 53 Putting n = 2m in the preceding equation gives 1 N-I x[m]= X[k]WNk Nk=o for m=0,1,...,N-1, which represents the even-valued indices of x[n], so 5 that the 2N-point IDFT of the zero-padded spectrum reduces to that of the N-point IDFT at the even-valued indices. The odd values of the index n represent the interpolated IDFT values between the original N-point IDFT values. This result generalises in an obvious fashion to temporal interpolation by a factor P, where P _2, simply by zero-padding the 10 spectral sequence with a sufficient number of zero-valued samples, prior to taking the IDFT. Thus, for high interpolation factors, there is a requirement to carry out a long IDFT following the zero-padding of the spectral sequence, 15 which involves a high computational load, whilst the application of DFT techniques to both the time-domain and frequency-domain data sets introduces Gibbs-type effects in both the time and frequency domains, due to the discontinuities arising at the ends of both the time-domain and frequency-domain data blocks. 20 A high-accuracy solution to the problem of DFT-based temporal interpolation is provided which involves the appropriate overlapping and weighting of the time-domain data blocks. The solution offers reduced complexity through the use of a 2-D formulation of the IDFT via the 25 Chinese Remainder Theorem (CRT). The ubiquitous Chinese Remainder Theorem plays an important role in many modern signal processing algorithms as it is able to generate 1-D to m-D and m-D to 1-D index mappings for the DFT/IDFT algorithms, which are both unique and cyclic in every dimension.

WO 00/16308 PCT/GB99/03075 54 Thus, application of the Chinese Remainder Theorem to DFT/IDFT indexing ensures : (a) the original DFT/IDFT can be recovered in its correct form, 5 and (b) the correct periodicities are obtained for each of the small DFTs/IDFTs that constitute the m-D form, and given the importance of the m-D formulation of the DFT/IDFT for the results to follow, it is worth seeing at this point how, via such index 10 mappings, the m-D formulation can be achieved. Suppose firstly that the length N of a 1-D IDFT can be written as N = N 1 x N 2 where N 1 and N 2 are relatively prime factors. Then by applying the 15 relatively prime modulo (RPM) and CRT index mappings : n [N2n 1 + N 1 n 2 ] mod N & k [N 2 t 2 k 1 + N 1 t 1 k 2 ] mod N to the input data and output data, respectively, for n 1 = 0,1,...,N 1 -1 and n 2 = 0,1,...,N 2 -1, where t 1 and t 2 are given by 20 N2t 2 1 modN 1 & N 1 t 1 1 modN 2 the IDFT can be written in a 2-D form as 1 NI-I[1 N2-1 '] -ntkI x[n, n 2 ] - X[k, k21 W2k2 "jkNJ N k,=0 N2 k 2 =0 with WN, and WN 2 the primitive N 1 -th and N 2 -th complex roots of unity respectively. 25 The above equation thus reduces the N-point IDFT, for the two factor case of interest, to a true 2-D IDFT, without twiddle factors, consisting of just two partial-IDFT processes. Generalisation to the m-D case is a straightforward task, with the resulting expression being more WO 00/16308 PCT/GB99/03075 55 commonly known as the prime factor algorithm (PFA). The same result can be achieved by reversing the order of the RPM and CRT index mappings or even by applying the CRT index mapping to both the input and the output data, although it is not necessary to do so. 5 Returning to the 2-D IDFT formulation of the preceding equation, it is now seen how this can be computed in a highly-parallel fashion to yield high-throughput solutions to the original IDFT algorithm, and hence to the inverse spectral interpolation (i.e. band-limited temporal 10 interpolation) algorithm. Firstly, by putting 1 "2-1 Y[kl,n 2 ]= X[k,,k]W 2k

N

2 k 2 =O the equation can be written as 1 N-I 15 x[nln 2 ]= Y[k,,n 2 Wnk I k 1 =0 This shows that the first partial-IDFT can be computed by means of

N

1 independent N2-point IDFTs, referred to hereafter as the row-IDFT process. Similarly, the second partial-IDFT can be computed by means of 20 N 2 independent N 1 -point IDFTs, referred to hereafter as the column-IDFT process. This decomposition of the IDFT into row-IDFT and column IDFT processes is known as the row-column method, with the 2-factor decomposition being as shown in Figure 14. Given the independence of the small IDFTs in the row-IDFT and 25 column-IDFT processes, high-throughput solutions can be simply obtained for the 2-factor decomposition via parallel computation of the N2-point IDFTs in the row-IDFT process, followed by parallel computation of the WO 00/16308 PCT/GB99/03075 56

N

1 -point IDFTs in the column-IDFT process. However, before the column-IDFT process can be carried out, the matrix containing the row IDFT output data must first be transposed. An alternative approach to implementing the row-column method is now discussed, which is 5 particularly appropriate to the inverse spectral interpolation problem, eliminating the need for matrix transposition of the row-IDFT output. It can be seen that each Nl-point column-IDFT is carried out upon a set of N 1 outputs from the row-IDFTs, where each Ni-point output set 10 corresponds to the same N2-point IDFT coefficient. Therefore, the first column-IDFT can be processed as soon as the individual N2-point zero sample components are output from the row-IDFT processor. The remaining Nl-point column-IDFTs can be similarly processed as soon as the individual components of the corresponding IDFT coefficient set are 15 computed. Thus, by processing the row-IDFTs coefficient-by-coefficient, rather than IDFT-by-IDFT, as is conventionally done, the transposition of the row-IDFT matrix output can be eliminated, as seen in Figure 15, with 20 the output from each row-IDFT coefficient generator feeding directly into a .corresponding column-IDFT. This idea is now applied to the inverse spectral interpolation problem where particular simplifications resulting from the choice and ordering of the index mappings lead to an attractive, computationally-efficient solution. 25 The inverse spectral interpolation algorithm described in this and the following sections is based on the technique discussed hereinbefore, where the number of input samples to the interpolator is denoted by K and the interpolation factor by P. The number of interpolated time-domain WO 00/16308 PCT/GB99/03075 57 outputs, N, produced by the interpolator, referred to hereafter as the inverse Sino spectral interpolation (ISSI) algorithm, is given by N=PxK where P is taken to be an odd integer and K is an integral power-of-two, 5 i.e. P and K are relatively prime. By choosing the interpolation parameters in this way, it is possible to carry out the N-point IDFT required by the interpolation algorithm as a row-IDFT process, comprising the computation of K independent P-point 10 IDFTs, followed by a column-IDFT process, comprising the computation of P independent K-point IDFTs. If, in addition, the input index mapping is taken as the RPM mapping and the row-IDFTs are processed coefficient-by-coefficient, rather than IDFT-by-IDFT, as discussed above, then the resulting algorithm can be represented by the flowgraph of Figure 15 16, where each P-point row-IDFT possesses at most one non-zero input sample. Thus, the row-IDFT process can be decomposed into a single zero sample process and 2(P-1) dual-node processes, as carried out by the 20 dual-node processor of Figure 17, where each dual-node processor computes two P-point IDFT coefficients from the same non-zero input sample via the relationship: W: = (W"-P)* for p = 0,1,...,P-1, where Wp is the primitive P-th complex root of unity 25 and where "*" stands for complex conjugation. The output from the zero-sample and dual-node processes can be fed directly into the K-point IDFTs of the column-IDFT process, which can, if required, be carried out by any conventional radix-2 FFT 30 algorithm. A further simplification can be achieved, however, by selecting WO 00/16308 PCT/GB99/03075 58 a K-point radix-2 FFT algorithm with naturally ordered input and bit reversed output, i.e. a decimation-in-frequency (DIF) algorithm, as the bit-reversal index mapping can be removed from the FFT algorithm and combined with that of the CRT index mapping, as required by the PFA 5 for reordering of the output data, to yield a composite index mapping for the reordering of the interpolated time-domain output. Thus, by regarding the zero-sample process as a special case of the dual-node process, the N-point IDFT reduces to a very simple, elegant 10 computational structure, comprising 12(P + 1) identical, independent processes, with each process being carried out by the Sino computing engine shown in Figure 18, which comprises one dual-node processor and two scrambled output, K-point FFT modules. 15 Finally, the inverse Sino spectral interpolation algorithm provides a very flexible approach to the temporal interpolation problem, as the same computational building block, namely the Sino computing engine, can be used to construct temporal interpolators of arbitrary precision, with only the input and output index mappings needing to be recomputed, together 20 with the Fourier matrix elements required by each dual-node processor. Also, the algorithm does not require the spectrum to be physically padded out with zero-valued samples, so that the memory requirement, for both input data and input index mapping, is kept to a minimum, whilst the memory requirement for the output data buffer is proportional to the 25 number of interpolated time-domain outputs to be produced, i.e. upon the interpolation accuracy required. Performance Analysis of Inverse Sino Spectral Interpolation Algorithm WO 00/16308 PCT/GB99/03075 59 To properly assess the performance of the inverse Sino spectral interpolation algorithm, it is necessary to consider the arithmetic complexity, communication complexity and structure of the algorithm for the particular implementation concerned, whether it be on a von Neumann 5 uni-processor machine, or a high-throughput, parallel implementation using an array of processing elements to pipeline the computation, as for example with a systolic array. The performance, in terms of throughput, of a particular 10 implementation can be expressed in terms of its time-complexity, which comprises components corresponding to : (a) its arithmetic complexity, i.e. total number of arithmetic operations, and (b) its communication complexity, i.e. total number and length of 15 data routing operations, with each being expressed as a number of time units, where one arithmetic time unit corresponds to the time taken to carry out a basic arithmetic operation, and one routing time unit corresponds to the time taken to move a sample of data one distance unit, where one distance unit 20 is the distance between consecutive processing elements. For the arithmetic complexity of the ISSI algorithm, which would be the prime performance metric for a uni-processor implementation, an appropriate measure is given by the total number of complex 25 multiplication/accumulation operations required, as this operation forms the major computational load of any technique based upon the FFT algorithm. One arithmetic time unit therefore corresponds to the time taken to carry out one complex multiplication/accumulation operation.

WO 00/16308 PCT/GB99/03075 60 For a data sequence of length K and an interpolating factor P, the arithmetic complexity, SA, of the inverse Sino spectral interpolation algorithm can be expressed as SA = 12(P-1)K + 12PKxlog 2 {K} ~ 12PKxlog 2 {K} 5 arithmetic time units. For the direct DFT solution, using an FFT algorithm but without exploiting the zero-valued samples, the arithmetic complexity, DA, involves 10 DA = 1/2PKxlog 2 {PK} = 1/2PKxlog 2 {P} + 12PKxlog 2 {K} arithmetic time units, which implies that the inverse Sino spectral interpolation algorithm involves S12PKxlog 2

{P

} 15 fewer time units than that of the direct DFT/FFT solution. Thus, from these figures, it is evident that in terms of the arithmetic complexity alone, the inverse Sino spectral interpolation algorithm is an attractive solution for performing temporal interpolation, 20 when compared to the direct DFT/FFT solution, with its relative computational efficiency improving as the interpolation factor P is increased. It possesses, however, the additional attraction of having an elegant, regular structure, making it an ideal candidate for very-large scale-integration (VLSI) implementation, as now discussed. 25 It is also important to consider the communication complexity and structure of the inverse Sino spectral interpolation algorithm.

WO 00/16308 PCT/GB99/03075 61 The inverse Sino spectral interpolation algorithm decomposes into an input index mapping, a number of independent and identical processing modules, and a composite output index mapping. Such a decomposition lends itself naturally to a parallel implementation, where the parallelism 5 can be achieved in one of two ways : (a) assigning up to Q + 1 computing engines, where Q = 1/2 (P 1), for carrying out the processing, these operating in parallel; and (b) building parallelism into the design of the computing engine, with the two column-IDFTs, for example, being carried out in parallel, 10 with maximum throughput being achieved by assigning Q + 1 parallel implementations of the Sino computing engine to the processing, although in practise a trade-off would have to be made between hardware size, cost and performance. 15 Given the independence of the Sino processing modules, the communication complexity of the Sino implementation is limited to that of an individual computing engine, which is dependent upon how each of the K-point, scrambled-output FFT modules is to be implemented. For a parallel, radix-2 solution, comprising 12K processing elements, for 20 example, where the processing element is a complex butterfly operator, the maximum distance that data samples must travel to reach their required processing element is 1 A4K distance units, whereas for the conventional FFT solution, comprising 12PK processing elements, the maximum distance is given by 14PK distance units. 25 Since the time-complexity associated with the routing of data over d distance units is -V routing time units, the time-complexity of the routing operations for the algorithm, PC , is given by PC ~ O(PK) routing time units, whilst that for the inverse Sino 30 spectral interpolation algorithm, SC , is given by WO 00/16308 PCT/GB99/03075 62 SC - O(-K) routing time units, with the result that the communication complexity of a parallel implementation of the direct DFT/FFT algorithm, i.e. one involving 12PK complex butterfly units, is approximately VI times greater than that of a parallel implementation of 5 the inverse Sino spectral interpolation algorithm, i.e. one involving 12K complex butterfly units for each K-point FFT. For the case where K is sufficiently small, i.e. K _ 32, each K-point IDFT may be carried out in a fully-pipelined, systolic fashion by means of 10 a linear array of K processing elements, where each processing element carries out the same operation of complex multiplication/accumulation. The row-IDFT output from the dual-node processor can be pipelined into the two length-K linear arrays at a constant rate, to give a time complexity for the Sino computing engine, and hence for the inverse Sino 15 spectral interpolation algorithm, of O(K), with each computing engine producing all 2K outputs in just K time units. Thus, the inverse Sino spectral interpolation algorithm possesses a number of properties, not found in more conventional solutions, which 20 make it an ideal candidate for VLSI implementation. The regularity of its structure, i.e. the repetitive use of a single building block - the Sino computing engine, makes it very easy to design temporal interpolators of arbitrary precision. This, combined with the relatively simple communication requirements, which are dependent purely upon the 25 implementational complexity of the computing engine itself, lead to a cost-effective solution when implemented with VLSI technology. In addition to the properties of the algorithms set out above, a dual channel processing technique is employed. The technique involves 30 processing two channels of sensor data simultaneously, starting with the WO 00/16308 PCT/GB99/03075 63 forward FFTs, for the time-to-frequency transformation, then the wideband compensation, finishing with the long inverse FFTs, for the frequency-to-time transformation. The dual-channel processing is achieved by making use of the real and imaginary parts of the processing 5 algorithm via straightforward packing and unpacking of the data. One channel is stored as the real data and another is stored as the imaginary data. Thus both channels can be processed simultaneously using the algorithms. 10 Suppose that the time-series data sequences for consecutive sensor channels, as required for input to the dual-channel processing, are given by {x[n]} and {y[n]} with time-frequency relationships given by: {x[n]} ---------- DFT --------- > { XR[k]+iXI[k]} {y[n]j ---------- DFT ---------- > {YR[k]+iYI[k]} 15 for n=0,1...,. k-1 and k=0,1...,. K-1. If the real and imaginary components of the complex-valued time series data sequence are written as: zr[n] = x[n] and zi[n] = y[n] with time-frequency relationship given by: {zr[n] + izi[n]l --------- DFT ----------- > {ZR[k]+ iZI[k]} 20 then the real and imaginary components of the DFT spectra of the real valued time series data sequences, {x[n]} and {y[n]} , can be written as: XR[k]= (zR[k] + ZR[K - k]) x[k]= (ZI[k] - ZI[K - k]) YR[k] = (ZI[k] + ZI[K- k]) 25 YR[k] I(zR[K- k] - ZR[k]) 252 WO 00/16308 PCT/GB99/03075 64 where XR(K - k] = XRk] X4K - k]= -xk] YR[K - k]= YR[k] 5 YI[K - k] -YI[k] with ZR[0] = ZR[K] and ZI[0] =ZI[K], so that unpacking a single spectral component involves two real multiplications and two real additions. The frequency-domain sensor data is preferably in unpacked form 10 in order to apply the compensation techniques for malfunctioning samples. The data can be repacked via the relationships: ZR[k] = XR[k]- YI[k] & Zlfk]= xT[k] + YR[k] so that repacking requires only two real multiplications and two real additions. 15 The skilled man will appreciate that the interpolation routines described hereinabove have applications beyond the field of acoustic sensors and beamformers. We may wish to seek protection for any of the novel and inventive features of the algorithms as separate concepts, 20 including the dual- processing nature of the algorithms using the real and imaginary components to process two channels simultaneously. Detailed Description of the Delay-and-Sum Sub-Unit. 25 Once the sampled data from the sensors has been upsampled, using an appropriate interpolation unit in the primary combining unit, and compensated for the presence of any malfunctioning hydrophones and erroneous samples if required, it remains only to apply appropriate time delays and shading coefficients, for each hydrophone channel and beam to WO 00/16308 PCT/GB99/03075 65 be steered, prior to the summation of the data. To achieve this, a delay and-sum sub-unit accepts interpolated hydrophone data, for each successive hydrophone channel, and applies the appropriate time delay prior to application of the beam steering coefficients. The summation 5 process then combines the samples for all 60 channels to yield the appropriate partial beam. The shading coefficient and time delay values are obtained by the primary combining unit from the appropriate block of memory associated with the primary combining unit. The time delays and coefficients are written into the memory by the SACM unit which may 10 update the values at any stage, for instance to compensate for changes to array shape. The space-time store which holds the sampled data from the hydrophones prior to interpolation holds the running time history of the 15 array as conventional beams are formed by generating address sets across the store and summing the corresponding sets of addressed data. The stored time history is of sufficient length to contain all the hydrophone data necessary for forming a beam when the incoming signal is at end fire, i.e. when the signal has to traverse the entire array. For a 960 20 hydrophone array, it is necessary that 9 blocks of compensated, interpolated hydrophone data are maintained in the store. Each block comprises 256 time samples with a twenty-five percent overlap before upsampling, and 7 times this amount of samples after upsampling, a new block of 7x192=1344 interpolated samples becoming available every 25 1/15th second approximately. For the most efficient processing of the data, three units of random access memory (RAM) are provided which are available for high speed access by the delay-and-sum sub-unit, denoted by S 1 , S, and S 3

:-

WO 00/16308 PCT/GB99/03075 66 (a) S, is an area of memory for storing the integer values k, representing the time delays for formation of the beams by the delay-and sum sub-unit. The integer values k, may most conveniently be stored in the form of an array of m rows of N integers, where m is the number of 5 channels processed by the primary combining unit, and N is the number of beams to be formed; (b) S 2 is an area of memory in which the shading coefficients are stored. As for the time delays, a separate set of shading coefficients is required for each beam to be formed, and may be stored in an array of 10 dimensions m by N. Thus, incrementing by N automatically steps through the correct shading coefficients; (c) S 3 is an area of memory in which the results of the beamforming process performed by the delay-and-sum sub-unit are stored. 15 In use, the delay-and-sum sub-unit accesses the time delay values and the shading coefficients from the memory blocks S, and S 2 using the correct address n. Beamforming is performed on all the samples in a block in the space-time store using a single set of time delay values and shading coefficients corresponding to a single beam. The process is then 20 repeated using a set of time delay values and shading coefficients for the next beam by incrementing the indices by N. Thus, the time delay values anrrd shading coefficients are only accessed by the delay-and-sum sub-unit every time a complete block of samples for a single beam has been beamformed, i.e. every 192 samples. 25 After the beamforming operation has been performed, the third set of memory S 3 contains a set of data values for each of the partially formed beams. The partially formed beam values may be stored in the form of an array of N rows (output channels) of samples, where each row of samples 30 is of duration t. Thus, when N partial beams are formed, N channels are WO 00/16308 PCT/GB99/03075 67 output into the memory. The ability to access the samples in the space time store, the required time delay and shading coefficient values and write the results to memory directly without waiting for information to be transferred from a remote unit maximises the speed of the processing. 5 For the array described herein, S, need only comprise up to 240 bytes, S 2 240 bytes and S 3 4320 bytes, where S 3 contains the results for every beam being formed. The delay-and-sum sub-unit applies shading coefficients to reduce the unwanted sidelobes. A nominal set of shading coefficients (i.e. for an 10 ideal equi-spaced linear array of hydrophones), Wm, are given by the exact Blackman shading function: (2n(m-1) / m-- 1 wm = a o - a 1 co2Z +a 2 co 4.zr aM -L( M for m = 1, 2... M, with coefficients given by: ao = 7938/18608, a, = 9240/18608 and a 2 = 1430/18608. 15 However, for accuracy, it must be assumed that the shading coefficients are not ideal, but are shape dependent which means that a different set of coefficients are needed to be applied to the hydrophone data for each beam to be steered. 20 In addition to the above delay-and-sum beamforming process, further improvements in accuracy can be attained by applying a linear interpolator to the output of the inverse Sino spectral interpolation routine, so that the contribution to the unshaded beamformer output for a 25 given hydrophone channel and a given beam steering look angle can be expressed as a weighted combination of successive hydrophone samples x[i] and x[i+ 11] as: a.x[i] + .x[i+1] where a, / e [0,11.

WO 00/16308 PCT/GB99/0307 5 68 The interpolation weights for each hydrophone channel and beam steering angle may further be combined with the shading coefficients, stored in memory

S

2 , and accessed by the primary combining unit using 5 the same unique address identifier as per the shading coefficients and time delays. Again, the updated coefficients, referred to as the beam steering coefficients, may be calculated and stored by the SACM unit. The Malfunctioning Sensor Detection Unit 10 The purpose of the malfunctioning sensor detection unit is to analyse the sampled data held in the space-time store to detect erroneous samples. These erroneous samples may arise from faulty sensors or due to faults in the ADC electronics associated with each sensor. The unit 15 generates a set of error codes which are.stored in a memory block. Typically, a failed sensor will produce a digitised output which is either in the fully low (all bits set to zero) or fully high (all bits set to one) state, which could swamp the partially formed beams from a primary 20 combining unit. Such an error situation, if ignored, will result in possibly severe degradation in the performance of the subsequent signal processing, thereby negating the processing gains achieved by the introduction of 25 sophisticated signal processing algorithms. The effect of the malfunctioning sensors and erroneous samples on the generation of beam patterns is to reduce the directivity index (or, more generally, the gain) of the array and to increase the sidelobe 30 responses of the array to directional interference. For a towed array, the WO 00/16308 PCT/GB99/030 7 5 69 directional interference may be from tow-ship generated noise, nearby shipping or for the case of active processing, strong acoustic reverberation caused by discrete bottom scatterers. However, before a malfunctioning sensor can be compensated, it must first be detected and 5 classified according to the size or multiplicity of the block of malfunctioning sensors to which it belongs. The beamforming system therefore further includes a suite of sensor classification and compensation algorithms, where the 10 compensation algorithms are applied according to whether the malfunctioning sensor belongs to one of the following four sets; (1) isolated malfunctioning sensors, (2) pairs of malfunctioning sensors, (3) blocks of three consecutive malfunctioning sensors, or 15 (4) blocks of more than three consecutive malfunctioning sensors. The automatic detection scheme is based upon the effect to the array gain when the contributions of individual sensors are excluded from 20 the array gain equation. Following the automatic detection routine, a compensation processing scheme based upon amplitude interpolation and phase extrapolation techniques is then applied, with appropriate modification (i.e. taking into account the multiplicity of the block of malfunctioning sensors to which it belongs) to compensate for each type 25 of malfunctioning sensor. The resulting compensation techniques are applied in the frequency domain, and can be applied to both narrowband and wideband signal processing systems. For the case of narrowband processing, where a 30 frequency-domain beamformer is to be assumed, the compensation WO 00/16308 PCT/GB99/030 7 5 70 techniques need only be applied to those temporal frequencies of interest. For the case of wideband processing, however, as in the beamformer described herein, wherein the beamforming is carried out by frequency domain techniques but is in essence a time-domain beamfirmer, the 5 compensation techniques will need to be applied to every positive temporal frequency, so that the reconstructed time-domain signal may be suitably obtained. Although the results are discussed for the general case of an 10 arbitrary array of sensors, the application of the automatic detection and compensation techniques is not limited to that of the primary combining unit structure of the parallel beamformer array described herein. The signal wave front impinging upon the array of sensors is assumed to be planar, i.e. a far-field source, and the signal and noise fields are assumed 15 to be independent and additive. The sensors are assumed to be acoustic sensors (i.e. hydrophones). In order to fully understand the operation of the malfunctioning hydrophone detection unit, it is necessary to make the following 20 assumptions: (a) that the primary combining unit is approximately linear, and (b) that the target is in the far field of the primary combining unit, i.e. at least three quarters of a kilometre from the local origin of the acoustic sub-unit for the size of array considered herein. 25 By means of these assumptions, the signal wavefront impinging upon the primary combining unit can be regarded as being planar, and linear interpolation (for amplitudes) and extrapolation (for phase) can thus be justifiably applied by the primary combining unit to generate estimated 30 sample values.

WO 00/16308 PCT/GB99/03075 71 The malfunctioning sensor detection unit works by analysing statistics about the sampled data derived from the space-time store. As stated hereinbefore, the space-time store comprises a set of time samples 5 taken from each channel fed to a primary combining unit. The samples in each row of the space-time store are divided into blocks of samples of temporal length tb, comprising 256 samples in this case. The space-time array is refreshed every time that a block of samples have been obtained, and the samples held in the store between each refresh are analysed for 10 statistics which are used to detect malfunctioning channels of data. The detection of malfunctioning sensors and erroneous samples in the space-time store comprises two phases. Firstly, an initial start-up routine (for obtaining an initial set of error codes when no reliable historical data is available) is applied. Subsequently, a standard detection 15 routine is used once the initial error codes have been generated. The operation of each of the two units is set out herein. Consider a space-time store comprising 60 blocks (one block per sensor channel) of 256 samples of the most recent samples, as shown in 20 Figure 11. The space-time store actually contains space for more than one block per channel to account for the maximum time delay needed for beamforming, but on start-up only a single block of uncompensated data will be present, the remaining blocks initially being empty. Assuming that the store is replenished with an overlap of twenty-five percent between 25 blocks (to minimise errors in interpolation), the remaining three quarters of the block are divided into six sub-blocks of 32 samples. The initial detection scheme comprises the steps of: (a) integrating the samples in a sub-block (in magnitude form) for each channel; WO 00/16308 PCT/GB99/03075 72 (b) calculating the median value of the 60 resultant integrated values; if the median differs from the mean of the previous time block, as calculated in step (f), then replacing the median with the mean; (c) calculating upper and lower detection thresholds from the 5 median value; (d) comparing the integrated value for each sensor with the upper and lower values of the detection threshold and producing a set of error codes for the samples in the block based upon the results of the comparisons. 10 The error codes may be stored in a memory as a sensor activity sequence. For the remaining five sub-blocks in the 256 sample block on each 15 channel, generation of the error codes for the sensor activity sequence involves the steps of: (e) integrating the 32 samples (in magnitude form) for each sensor in the primary combining unit; and for each sensor; (f) calculating the mean of the remaining integrated values - i.e. 20 excluding the contribution from the sensor under test; if the mean differs significantly from the median, as calculated in step (b), then replacing the mean with the median; and (g) calculating the upper and lower detection thresholds from the mean (or median); and 25 (h) testing the integrated value for the sensor under test against the upper and lower detection thresholds, with the test results for each sub-block being used to update the error codes stored in the sensor activity sequence for that block.

WO 00/16308 PCT/GB99/03075 73 As the sensor activity sequence is updated every time a block of data in the space-time store of length 256 samples is produced, a malfunctioning sensor can be detected and flagged up in the activity 5 sequence every time 192 samples (due to the twenty-five percent overlap) are beamformed. Since the beamformer produces formed beams in batches by analysing in blocks of 256 samples, the values of the error codes held in the sensor activity sequence are effectively updated every time the block of 256 samples is passed to the beamformer. 10 The initial malfunctioning sensor detection algorithm set out above relies upon the fact that if at least fifty percent of the hydrophones are functioning correctly, the median value can be used as a reference for detecting faulty sensors producing erroneous samples. Once an initial picture has been developed using the initial routine, and a first set of 15 error codes are produced and written to the memory, a more accurate scheme is employed whereby single sensors are analysed by comparing the mean of the samples to the remaining sensors in a block. This is known as the standard malfunctioning sensor detection routine. 20 The sensor activity sequence containing the error codes comprises a block of digital data, with one bit being allocated per channel in order to allow a unique error code representing blocks of 1, 2 or 3 malfunctioning sensors to be indicated. A suite of compensating algorithms are (optionally) incorporated into the primary combining unit which modify 25 the analysed blocks of sampled data held in the space-time store dependent upon the values stored in the sensor activity sequence so as to substitute erroneous samples with estimated sample values. These values are derived by the algorithms based upon the outputs of samples surrounding the erroneous samples. In order to compensate for three types of error, WO 00/16308 PCT/GB99/03075 74 whereby 1, 2 or 3 adjacent sensors are in error in each block of samples, three algorithms are therefore used. It is possible to omit the malfunctioning sensor detection stage and 5 simply pass the digitised sample data straight to the interpolation stage of the primary combining unit to increase processing speed. Where a malfunction detection unit is provided, means may be provided to allow the operation of the malfunction detection unit to be manually or automatically switched on or off as required. Thus, where data from a 10 primary combining unit is known to be highly reliable, the detection routine may be switched off, for example manually switched off by the person using the system, for that primary combining unit, in order to reduce the number of processing steps required. 15 Automatic Detection of Malfunctioning Sensors So far, only brief details of the method for detecting erroneous samples from the statistics of the samples in the store have been provided. The error detection scheme based upon statistical analysis can be realised 20 as follows. The detection scheme comprises a standard and an initial detection routine. (a) The Standard Detection Routine 25 For an array comprising M sensors, with 2/2 spacing, the array output power, i.e. the output of a conventional beamformer, denoted Ps, can be written as: P =M 2 S+MN where S is the signal power at any one sensor, and N is the noise power. 30 The signal-to-noise ratio (SNRA) at the array output can then be written as: WO 00/16308 PCT/GB99/03075 75 SNR A = P - P ) / P where (P, - P.) is the array signal power and P. is the array noise power. If the array gain, AG, is defined as the ratio of the array signal-to-noise ratio, SNRA, to that of a single sensor, SNR,, then: 5 AG = SNRA / SNRH - = M i.e. the array gain is equal to the number of sensors M in the array. Suppose that the phase of the signal is assumed to be perfectly coherent across the array, but that the amplitude is perfectly random, 10 obeying a Gaussian distribution. The noise, as already stated, is assumed to be of uniform amplitude but totally incoherent across the array. Then after conventional beamforming, the signal-to-noise ratio for an array of M-1 correctly functioning sensors is given by: SNRM_ = (M-1)S / N 15 where the noise power is as previously defined, but where the signal power is now defined as the average power of the M-1 correctly functioning sensors. Suppose, in addition, that the remaining sensor is assumed to 20 possess the same noise power, but a signal power of arbitrary value, So. Then if the array output consists of a contribution from all M sensors in the array, the array output SNRM can be written as: SNRM ([(M- 1) 2 S + 2(M- 1)(SSo) 1 2 + So])/MN which reduces to: 25 SNR M =(MS)/N as required when So = S.

WO 00/16308 PCT/GB99/03075 76 The question now to be addressed, is at what level does the signal power, So, result in the condition: SM, 1_ SM , i.e. at what level of So can the corresponding sensor be considered to be malfunctioning. To answer this, let: 5 a = SNRM_1 / SNRM & So = ,2S The parameter ,8 is a detection coefficient that relates the signal pressure amplitude of the odd sensor to that of the average amplitude of the remaining M-1 correctly functioning sensors. Then: 10 f8=(M-1){[M / (M-1)a]2 -1} i.e. f can be expressed in terms of the SNR parameter a and the total number of correctly functioning sensors in the array. Through manipulation of the above equations, it is clear that when 15 a = 1, i.e. when SNRu.

1 = SNRM, the value of the lower detection coefficient, f8, is such that 84 1/2 as M increases. Thus, if the amplitude level of the odd sensor is being set so low as to measure less than 1/2 that of the average M-1 correctly functioning sensors, its effect is to decrease the array output SNR, so that the corresponding sensor may be considered 20 to be malfunctioning. This typifies the behaviour of conventional analogue systems, where the data generated from malfunctioning sensors generally defaults automatically to zero. A more serious problem may arise with the processing of digital 25 data, however, where the amplitude may default, with equal probability, to either the minimum or the maximum level. For example, when the signal power So, is such that a 1/2, i.e. when SNRM = 2xSNRM.I, the value of the upper detection coefficient, ,8, is 0.414M. Thus, if the amplitude level of the odd sensor is greater than 0.414M times that of the WO 00/16308 PCT/GB99/03075 77 average of the remaining M-1 correctly functioning sensors, its effect is to double the array output SNR, so that the corresponding sensor may also be considered to be malfunctioning. 5 To start up the standard detection procedure, just discussed, there is a problem in that there is no reliable knowledge as to the initial working status of the individual sensors, so that it is not possible to compare the amplitude of a given sensor against that of the average of the remaining sensors, as they may not all be correctly functioning. One very 10 large sensor amplitude, for example, could result in a grossly incorrect value for the average amplitude level. Thus, an initial start-up detection procedure is used for initial direction of malfunctioning sensors. This does not depend upon the use of a computed average, which may or may not be contaminated by the presence of malfunctioning sensors. 15 (b) The Initial Detection Scheme As described hereinbefore, the approach employed is to use the median value of the sensor amplitudes for the purposes of sensor 20 comparison, since provided at least half of the sensors are correctly functioning - which would certainly be expected to be the case, barring a catastrophic failure - the median value will at least correspond to the amplitude of a good, correctly functioning sensor. The upper and lower amplitude detection coefficients could be chosen to be less stringent in 25 this initial detection routine to allow for the possible discrepancy in performance resulting from the choice of the median amplitude instead of an average for the array. Any malfunctioning sensors which get through the initial detection stage undetected, because of the reduced thresholds, will be detected and flagged up in the standard detection procedure which 30 follows it.

WO 00/16308 PCT/GB99/03075 78 Note that for both the initial and standard detection procedures, given an appropriate choice of upper and lower detection coefficients for carrying out the amplitude comparisons, the individual sensors are flagged 5 up as being either correctly functioning, for which: flag = 0 or malfunctioning, for which: flag = 1 these flags being used to form a sensor activity sequence for the 10 corresponding primary combining unit. This sequence is in turn used to drive a set of simple compensation algorithms. On entry to the standard detection procedure, all those sensors flagged up as being correctly functioning are used in the derivation of the average sensor pressure amplitude for use in the amplitude comparison routine. 15 The fact that the sensor data being processed is digital is of the utmost importance, since it means that the amplitude of the data from each malfunctioning sensor may default, with equal probability, to either the minimum or the maximum level. This is in contrast to an analogue 20 system, where such data will generally default to the minimum level. To see why the sensors behave in this way, one need only look at the workings of the internal electronics of a digital processing system. For example: (a) an open gate input prevents a signal on that input from 25 getting to the gate circuitry and thus to the output; in a CMOS (complementary metal oxide semiconductor) gate an open gate input has an unpredictable effect; (b) an open gate output prevents a signal on any of the inputs from getting to the output, causing a break between the gate and the WO 00/16308 PCT/GB99/03075 79 output pin of the circuit, and resulting in no voltage on the output pin; and (c) an input or output shorted to the supply voltage will be stuck in the "high" state, i.e. at the maximum level, whereas an input or output 5 level shorted to the ground will be stuck in the "low" state, i.e. at the minimum level. Thus, intermittent openings or shorts in the sensor or array circuitry may well cause unpredictable behaviour in the sensor output, 10 depending upon where the broken connection occurs within the circuit. Apart from obvious physical causes such as these, problems might also arise as a result of telemetry errors, with "bad" sensor data being produced either intermittently or persisting over a prolonged period. However such data is produced, whether through a physical fault or 15 otherwise, given that spatial data and spatial-frequency data (i.e. beam patterns), like Fourier transform pairs, are broad in one domain when sharp or delta-like in the other, it needs just one sensor within the array to be malfunctioning in this fashion, for the spatial information provided by the conventional beamformer to be corrupted or possibly even 20 destroyed, and target tracks subsequently lost. Finally, it is possible for sensors to be incorrectly flagged as malfunctioning as a result of short-term signal activity, i.e. as a result of a short-term signal loss or a short-term signal burst. To cater for this 25 situation the standard detection procedure can be carried out on a continuous basis in order to provide the suite of compensation algorithms with the necessary sensor activity sequence. In addition, the initial detection procedure can be reapplied periodically to ensure that those sensors incorrectly classified can have their activity flag correctly reset. 30 If this check is not carried out, and blocks of sensors are misclassified in WO 00/16308 PCT/GB99/03075 80 this way, then array gain from the conventional beamformer may be unnecessarily lost. Compensation of Malfunctioning Sensors 5 Having identified erroneous samples, the block of samples may then processed to replace erroneous values with estimated data. The compensation techniques that have been developed are carried out in the frequency domain, i.e. following time-frequency transformations of the 10 time-domain sensor data via the FFT routine. The resulting sensor data is therefore to be of complex form, processing both amplitude and phase information. The complex sample, Xm, corresponding to the m-th sensor in the primary combining unit, for an arbitrary temporal frequency (i.e. for an arbitrary temporal FFT bin), is represented here in the form: 15 Xm = A exp(iDm) where Am is the amplitude: Am = IXmj and ~m the phase as denoted by: -, = arctan(x m / XRm) 20 where XRm is the real component of the complex sample and XIm the imaginary component. The estimated value of the sensor output for the case of a malfunctioning sensor is to be denoted by: Xm\ A i X. = A, exp(iew) to distinguish it from the output of a correctly functioning sensor. 25 Case 1 - Compensating Isolated Malfunctioning Sensors WO 00/16308 PCT/GB99/03075 81 For the case of an isolated malfunctioning sensor with index m, the amplitude of its output can be estimated by averaging those of its correctly functioning neighbours, A,,, _, and Ar+,: 5 = 2(A,. I + A,+, ). To obtain an estimate of the phase, the phase difference: (,n-1I rm-2 is first calculated from the two correctly functioning sensors to the left of the malfunctioning sensor. This phase difference is then added to the 10 value of the phase at the (m-1)-th sensor to give: Om = (D.-I + ((D.-2 - (D.-I) as required. Alternatively, this phase estimate could have been derived from the two correctly functioning sensors to the right of the malfunctioning one, to give: 15 =,- = (m+I - ((mn+2- (m+l) Thus, for the case of an isolated malfunctioning sensor, the amplitude estimate is obtained via linear interpolation of its nearest neighbours, whilst the phase is obtained by extrapolation from the two neighbours to its left or to its right, thereby avoiding the ambiguity 20 problem that conventionally arises in the computation of the phase estimate. Case 2 - Pairs of Malfunctioning Sensors 25 For the case of a pair of malfunctioning sensors, with indices m and m + 1, the amplitudes can be estimated by: ^1 A,,, = (2A,- + Am+2 and3 and WO 00/16308 PCT/GB99/03075 82 . +\ 1A Am+i (Am , +2Am+ 2 ) 3 via linear interpolation between the nearest functioning neighbours. The phase estimates, in turn, are given by: m m-1 +m-1 m-2 5 and Dm+I = (m+2 - (L'm+ 3 - 'Dm+2 with (Dm being obtained via extrapolation from the two functioning neighbours to the left, and (m,.+I being obtained via extrapolation from the two functioning neighbours to the right. 10 Case 3 - Blocks of Consecutive Malfunctioning Sensors Generalising, for the case of a block of more than two consecutive malfunctioning sensors, with the first and last indices given by ml and 15 m2, the amplitudes can be estimated by: Am, = (LAm, + Am2+1 /11 Am+, = - ((L-1)Am,_ +2Am2+1) 1 Am2-1 = (2Am 1 - + (L- 1)Am2+ 1 20 A.2 = (A., - + LAm2+1 ) L +1 via linear interpolation between the nearest functioning neighbours, where L = m2 - ml + 1, 25 WO 00/16308 PCT/GB99/03075 83 i.e. L corresponds to the number of malfunctioning sensors in the block. To obtain the corresponding phase estimates, the compensation algorithm, as for the case of a pair of malfunctioning sensors, works form both ends of the block. 5 Thus, the phase values of the two end sensors in the block: ml ml-1 m -1 m -2 and (m2 (Dm2+l - ((m2+2 - (Dm2+1) 10 are estimated first, followed by the neighbouring malfunctioning sensors: (1i+i + A -4mi-i) and m2-1 m2 - ((Dm2+1 m2 etc, exploiting the previous pair of phase estimates each time, until all the 15 malfunctioning sensors in the block have been processed. Discussion of Compensation Routines In compensating for malfunctioning sensors, the rule to be adopted 20 with the application of the techniques described herein is that all the isolated malfunctioning sensors are processed first, followed by all the pairs of malfunctioning sensors, followed by all the blocks of three malfunctioning sensors. This enables the phase estimates for the isolated malfunctioning sensors to be used for the processing of the pairs of 25 malfunctioning sensors, and the subsequent phase estimates for the pairs of malfunctioning sensors to be used for the processing of the blocks of three malfunctioning sensors. For the particular case of the isolated malfunctioning sensors, the rule will be that if correctly functioning neighbours are available to the left of the malfunctioning sensor, then they WO 00/16308 PCT/GB99/03075 84 be used first; otherwise, those to the right will be used. For blocks of more than three malfunctioning sensors, however, the only compensation to be applied will comprise zeroisation of the data. 5 Performance Considerations of Compensating Algorithm Suite At those temporal frequencies where the sensor spectra are dominated by spatially coherent signals, this compensation technique, at least for cases 1, 2 and 3 (i.e. for isolated malfunctioning sensors, pairs 10 of malfunctioning sensors and sets of three consecutive malfunctioning sensors), will provide a good approximation to that obtained if the sensors were actually functioning. Thus, as a result, the directivity index of the array will be maintained to a level close to the theoretical value, and the beam sidelobe responses to the directional narrowband signal will be 15 significantly reduced. At those temporal frequencies where no strong directional signal is present, the estimated amplitude spectrum will approximate that of neighbouring sensors, whilst the phase spectrum, relative to neighbouring sensors will be random. 20 The resulting lack of correlation between the phase of the estimated spectrum and that of adjacent functioning sensors reflects the noise-like nature of the acoustic field at these temporal frequencies, whilst the amplitude will also, on average, be the same as if the malfunctioning sensor was actually fully operational. 25 The techniques discussed for malfunctioning sensor compensation have been presented in the context of narrowband processing, where just a fraction of the total temporal frequency bandwidth - possibly a single temporal FFT bin - is of interest. However, the results can be extended 30 to the case of a wideband processing system, i.e. to one where every WO 00/16308 PCT/GB99/03075 85 temporal frequency is of interest. This is achieved by applying the compensation techniques to every positive temporal frequency bin covering the passband and transition band regions, resulting from application of the FFT routine, for every functioning sensor necessary for 5 carrying out the amplitude interpolation and phase extrapolation techniques. If the conventional wideband beamforming is to be carried out in the time domain, the negative frequency components are set up via the 10 complex conjugate symmetry property of the spectral components of each channel, comprising either the real or imaginary components of the complex input data, and the resulting two-sided compensated spectrum passed through an inverse temporal FFT routine, to yield the required reconstructed, wideband signals at the malfunctioning hydrophones. 15 The ability to identify erroneous samples allows the compensation of data to be performed in real time which has not previously been possible. Arithmetic Complexity Requirements for Treatment of Erroneous 20 Data To estimate the arithmetic complexity involved in correcting the blocks of sampled data, it is assumed that at any one time, a total of 5% of the sensors in an acoustic sub-module will on average be 25 malfunctioning, this involving a total of approximately 20% of the sensors for each of the primary combining units requiring transformations both to and from the frequency domain.

WO 00/16308 PCT/GB99/03075 86 The time-to-frequency and frequency-to-time transformations are carried out by means of a 256-point FFT algorithm, this involving 1740 real multiplications and 5008 real additions for the case of complex data. 5 The total FFT requirement for carrying out the compensation of malfunctioning hydrophones within a single primary combining unit of 60 hydrophones is therefore given by: M, = 12 x 1740 x 2800 / 256 = 228,375 10 real multiplications per second, and Aff, = 12 x 5008 x 2800 /256 = 657,300 15 real additions per second, where the sampling rate, as defined, is 2800 Hz. Suppose now, that the compensation of a single malfunctioning sensor, for a single temporal frequency, involves a total of 7 real 20 multiplications, 5 real additions, 1 square-root operation, 2 arc-tangent operations and 2 trigonometric operations - these figures have been achieved in practice in implementing the compensating routines set out hereinbefore. The arc-tangent and trigonometric operations are efficiently carried out via look-up table techniques, and can therefore be excluded 25 from the arithmetic complexity figures. The complexity of the square root operation, on the other hand, can be realistically equated to that of four real multiplications, so that the total arithmetic complexity requirement for carrying out the compensation of a single malfunctioning sensor, for a single temporal frequency, can be approximated by just 11 30 real multiplications and 5 real additions.

WO 00/16308 PCT/GB99/03075 87 Thus, the total arithmetic complexity requirement for carrying out the amplitude interpolation and phase extrapolation of the sensor data, for all positive temporal frequency bins, within a single acoustic sub-module, 5 is given by: Mcom = 3 x 128 x 11 x 2800 / 256 = 46,200 real multiplications per second, and Aco,, , = 3 x 128 x 5 x 2800 / 256 = 21,000 real additions per second, where it is assumed that 5% of the sensors in 10 the acoustic module are malfunctioning. Clearly, when more than 5% of the sensors are malfunctioning, a proportionately greater number of arithmetic operations will be involved, with an acceptable worst-case figure of 50%, for example, involving up to ten times the numbers of operations given above. As already stated, if more than 50% of the 15 sensors are malfunctioning, then only the first 50% will be fully compensated for, the remaining simply having their data zeroised. Array-Shape Compensation 20 Much of the potential performance loss resulting from deformation of the array shape during use (i.e. as when the array is towed behind a ship which makes a manoeuvre) can be recovered using a module which determines the position of each sensor in the array. Based upon this positional information, the required time delays for each sensor can be 25 recalculated by the SACM unit and the look-up tables stored in memory which are accessed by the primary combining units can be updated accordingly with new values. Exact measurements of the sensor positions can-be made for each sensor channel, or the positions can be estimated using measurements of only selected sensors. 30 WO 00/16308 PCT/GB99/03075 88 Ideally, a single heading and depth sensor should be provided for each acoustic sensor or hydrophone. This allows compensation for changes in spacing between the sensors and changes in the elevation of the sensors. However, such an approach is expensive and so, in an 5 alternative scheme, a polynomial interpolation curve such as a spline can be used to estimate hydrophone positions to a good accuracy with a single heading and depth sensor being provided for each acoustic module. To further improve the accuracy in such an approach, the time-history of the array shape can be used in order to refine the estimation based on a 10 hydrodynamic model of the media through which the array is towed. The ability to compensate for array shape changes in real time is further enhanced by making the assumption that each of the sets of sensors associated with each primary combining unit is approximately 15 linear (i.e. the sub-units of sensors are spatially fixed). Thus, it is assumed that only the blocks of sensors associated with each acoustic module moves relative to one another with the sensors in each block being fixed relative to one another. This considerably reduces the number of individual calculations required, making real-time array shape 20 compensation a reality. Performance Analysis of Beamformer Unit. To assess the performance of the above described beamformer against mathematically generated, 18-bit synthetic digitised sample data, 25 suppose firstly that the sampling rate for a 32-sensor linear array is given by 4096 Hz, and that the maximum temporal frequency of interest is given by 1024 Hz, i.e. one quarter of the sampling rate. As described hereinbefore, a 256-sample input data block is adopted, with a full Hanning window function, an overlap of twenty-five percent between 30 blocks and an interpolation factor of 5. The worst case root-mean-square WO 00/16308 PCT/GB99/03075 89 (RMS) interpolation error is given as -103 dB which corresponds to the quantisation noise level of an 18-bit signal processing system. The level of the interpolation error falls off at a rate of 6 dB per octave with decreasing temporal frequency. 5 A set of simulated mathematical data corresponding to a far-field target, at a bearing of 1200 to the 32-sensor array, with temporal frequency content dominated by the presence of a single 880 Hz tone embedded in Gaussian white noise, was then input to the beamformer. 10 The array of 32 sensors is assumed to contain a block of two malfunctioning sensors close to one of its ends, and a block of three between the array centre and its other end. The two blocks are malfunctioning in the low or the high state. 15 The effects of the application of the malfunctioning sensor detection and compensation unit on the beams formed by the beamformer are shown in Figures 12 and 13. The output of the beamformer is unshaded, and interpolation is performed via frequency-domain techniques 20 using a suite of FFT algorithms. Each figure contains two traces, one for an uncompensated output, the other showing the output after compensation. The output beams are also subject to linear interpolation and temporal averaging in beam space to yield greater statistical reliability. 25 Figure 12 shows the output of the beamformer when the input data is "clean" having a SNR of 100 dB at the array. The red trace is produced with the compensating unit switched off, the blue trace with the compensating unit switched on. A clear improvement in the beams is 30 apparent. Similar improvements are apparent for a set of "noisy" input WO 00/16308 PCT/GB99/03075 90 data shown in Figure 13. The proposed system is thus clearly shown to be capable of processing arrays which contain faulty sensors. In conclusion, it will be understood that the present invention lies, 5 in one aspect, in providing an improved method and apparatus for beamforming which incorporates a plurality of partial beamformers adapted to process incoming sensor data in parallel, before combining the parallel outputs to produce a set of full beams for the array. The invention also resides in the provision of a means for determining the 10 operational status of the sensors and the sampled data in real time, and compensating for any erroneous samples. Furthermore, the invention resides in the production of a fast and simple method of forming beams by tailoring the size of the blocks of information required by the individual beamforming sub-units to process the beam to match the amount of 15 available memory present on a hardware board or unit. The beamformer is self-compensating and provides for wideband or narrowband beams to be formed.

Claims

1. A beamformer technique for use in combination with an array of sensors, each sensor being adapted to produce a channel of output data 5 samples in response to a measurand, characterised by the steps of: (a) storing a block of data samples from m sensors in a first store means; (b) analysing the data samples in the first store means to produce statistics which are interpreted to determine whether the data 10 from a sensor is erroneous; and (c) providing a fault indicating means adapted to store an error code indicative of the function of each sensor in the array in a second store means. 15

2. A beamformer technique according to Claim 1 characterised in that the sampled data is stored in a memory and is in the form of a space-time store comprising m rows of sampled data, each row corresponding to a block of sampled data of temporal length t from a respective sensor. 20

3. A beamformer technique according to Claim 1 or Claim 2 characterised in that step (b) comprises the additional step of dividing the data stored in a row of the store into sub-blocks of data, and analysing each sub-block to produce the statistics which are interpreted to determine whether the individual samples are erroneous. 25

4. A beamformer technique according to any preceding Claim characterised in that the fault indicating means stores the error codes in a sensor activity sequence, the activity sequence comprising a number of error codes which are referenced by an address which corresponds to a 30 unique channel block of data. WO 00/16308 PCT/GB99/03075 92

5. A beamformer technique according to any preceding Claim characterised by further comprising a first or initial detection technique comprising the steps of measuring the median value of the integrated 5 samples in a sub-block, and comparing the individual integrated sample values in the sub-block to the median value to determine if the sample data from the sub-block is erroneous.

6. A beamformer technique according to any preceding claim 10 characterised by the further step of replacing samples of data which are indicated as erroneous by the fault detection means by estimated samples derived from the correct samples neighbouring the faulty samples.

7. A beamformer technique according to any preceding claim 15 characterised by the additional step of replacing a block of sampled data by zero values where there are more than a predetermined maximum number of consecutive sensors producing erroneous data.

8. A beamformer technique according to any preceding claim 20 characterised by the additional step of updating the data samples contained in the store on each sampling interval by adding new samples to one end of-the store and removing old samples from the other temporal end of the store so that the sampled data "ripples" across the store over time. 25

9. A beamformer technique according to any preceding claim characterised by further comprising a second detection process used after an initial set of sensor activity sequence data has been derived using the first or initial detection process. WO 00/16308 PCT/GB99/03075 93

10. A beamformer technique according to Claim 9 characterised in that the second detection process comprises, for each channel of sampled data, the steps of: 5 (a) calculating the signal-to-noise ratio for the whole array of sensors; (b) removing one channel of data and calculating the signal-to noise ratio of the remaining sensors, after the channel has been removed, (c) comparing the two measured signal-to-noise ratios, and 10 (d) determining whether or not a channel of samples is faulty based upon the results of the comparison.

11. A beamformer technique according to Claim 10 characterised by further comprising comparing the value of the signal-to-noise ratio for the 15 removed channel with that of the remaining correctly functioning channels.

12. A beamformer technique according to any one of Claims 1 to 10, characterised in that the individual channel is compared with the signal-to 20 noise ratio of the remaining correctly functioning sensor channels using a block of pre-interpolated data from the store.

13. A beamformer technique according to any preceding claim which operates in real time to detect malfunctioning sensor data samples and 25 calculates estimated data to replace the erroneous samples prior to input to a beamforming routine.

14. A self-compensating beamformer apparatus for use in combination with an array of sensors, each sensor being adapted to produce a channel 30 of output data samples in response to a measurand, comprising: WO 00/16308 PCT/GB99/03075 94 (a) a first store means adapted to store a block of data samples from m sensors; (b) a sample analysing means adapted to analyse the data samples in the store means to produce statistics and interpreting means 5 adapted to interpret the samples to determine whether the data from a sensor is erroneous; and (c) a fault indicating means adapted to store an error code indicative of the function of each sensor in the array in a second store means. 10

15. A self-compensating beamformer apparatus according to Claim 14 characterised by further comprising a memory means adapted to define a space-time store, the store comprising m rows of the sampled data, each rows corresponding to a block of sampled data of temporal length t from a 15 respective sensor.

16. A self-compensating beamformer apparatus according to Claim 14 or Claim 15 which further comprises means for dividing the data stored in a row of the store into sub-blocks of data, and analysing means adapted to 20 analyse each sub-block to produce the statistics which are interpreted to determine whether the individual samples are erroneous.

17. A self-compensating beamformer apparatus according to any one of claims 14 to 16, characterised in that the fault indicating means defines a 25 sensor activity sequence in which the error codes are stored, the sensor activity sequence comprising a number of error codes which are referenced by an address which corresponds to a unique channel block of data. WO 00/16308 PCT/GB99/03075 95

18. A self-compensating beamformer apparatus according to any one of Claims 14 to 17 characterised by further including measuring means adapted to measure the median value of the integrated samples in a sub block, and first comparison means adapted to compare the individual 5 integrated sample values in the sub-block to the median value to determine if the sample data in the sub-block is erroneous.

19. A self-compensating beamformer according to any preceding claim which further includes means for generating estimated samples from the 10 correct samples neighbouring the faulty samples, and means for replacing the faulty samples with the estimated samples.

20. A self-compensating beamformer according to any preceding claim, characterised by further comprising, associated with each channel: 15 a) calculating means adapted to calculate the signal-to-noise ratio for the whole array of sensors; b) removal means adapted to remove one channel of data, the 20 calculating means being adapted to calculate the signal to noise ratio of the remaining sensors; c) second comparison means adapted to compare the two measured signal to noise ratios; and 25 d) processing means adapted to determine whether or not a channel of samples is faulty based on the results of the comparison.

21. A self-compensating beamformer according to any one of Claims 30 14 to 20 characterised by further comprising: WO 00/16308 PCT/GB99/03075 96 a) an interpolating means adapted to interpolate the sampled data prior to beamforming; and 5 b) a decimating means adapted to apply an appropriate decimation factor to the output of the beamformer, in which the decimation factor is dependent upon the number-theoretic properties of the associated processing parameters so that the beamformed data is produced at an 10 output rate which is different to the rate at which the data is sampled from the sensors and different from the rate at which the data is interpolated.

22. A beamformer for use in combination with an array of sensors, each sensor being adapted to produce a channel of output data samples in 15 response to a measurand, in which each beam is produced as a stream of beamformed output samples comprising: (a) an interpolating means adapted to interpolate the sampled data prior to beamforming; (b) a decimating means adapted to apply an appropriate 20 decimation factor to the output of the beamformer, in which the decimation factor is dependent upon the number-theoretic properties of the associated processing parameters so that the beamformed data is produced at an output rate which is different to the rate at which the data is sampled from the sensors and different to the rate at which the data is interpolated. 25

23. A method of beamforming for use in combination with an array of sensors, each sensor being adapted to produce a channel of data samples in response to a measurand, the method comprising the steps of: (a) interpolating the sampled data prior to beamforming; WO 00/16308 PCT/GB99/03075 97 (b) applying an appropriate decimation factor to the output of the beamformer, the decimation factor being dependent upon the number-theoretic properties of the associated processing parameters so that the beamformed data is produced at an output rate which is different 5 to the rate at which the data is sampled from the sensors and different to the rate at which the data is interpolated.