WO1996030780A1

WO1996030780A1 - Method of estimating missing data points at known positions in an input signal

Info

Publication number: WO1996030780A1
Application number: PCT/GB1996/000691
Authority: WO
Inventors: John Bodonyi
Original assignee: The Commonwealth Of Australia
Priority date: 1995-03-24
Filing date: 1996-03-22
Publication date: 1996-10-03
Also published as: GB9506041D0; GB2299214A; AU5151996A; GB2299214B

Abstract

A method of estimating missing data points at known positions in an input signal comprises the steps of: a) Fourier transforming the input signal from one of the frequency or time domains to the other; b) setting a threshold in the other domain so as to separate at least one unambiguous component of the transformed signal from the remaining components including noise; c) inverse Fourier transforming the separated component back to the original domain; d) replacing the missing data points at the known positions in the input signal with data comprising the values of the inverse Fourier transformed separated components at the known positions to generate an improved signal. The steps (a) to (d) are iterated to achieve a desired level of input signal improvement.

Description

METHOD OF ESTIMATING MISSING DATA POINTS AT KNOWN POSITIONS IN AN INPUT SIGNAL

This invention relates to a signal processing method for estimating missing data points in an input signal.

The method is of particular use in receiving systems which comprises a multi channel antenna array and uses digital receiving beam-forming techniques. Within such systems the failure of receiving channels increases side lobe levels. It has in the past been proposed to substitute data from a failed channel with data from an adjacent channel or to simply interpolate from working channels on either side of the failed channel. However such methods have been found beneficial only at low operating frequencies where the phase difference between antenna elements was less than 60 degrees, whilst at the medium and high frequencies the side lobe levels worsened as a result of inaccurate estimates.

This invention provides a method of estimating missing data points in an input signal which comprises the steps of:

a) Fourier transforming the input signal from one of the frequency or time domains to the other,

b) setting a threshold so as to separate at least one unambiguous component of the transformed signal from the remaining components,

c) inverse Fourier transforming the separated component to the one domain, and d) replacing the missing data points in the input signal with data from the separated component to generate an improved signal. The improvement in side lobe levels effected in any single iteration of the method depends on the proportion of missing data points. In an example with 10% of the data point missing a reduction of 20dB in side lobe levels was found as compared to the case when no correction was applied.

As mentioned above, depending on the proportion of missing data, steps (a) to (d) may be iteratively repeated on the improved signal, except that in step (b) a threshold is adeptively set which is lower than that in the preceding iteration so as to separate a greater number of unambiguous components.

In this way an improved estimate of the value of the missing data can be deduced after each iteration.

Depending on the processing power available, the algorithm automatically lowers the threshold in step (b) in each iteration until it is no longer possible to separate any further unambiguous components.

Preferably the threshold in step (b) is set by deducing a mean or median value, of the amplitude of the noise component of the transformed signal and adding a margin to that average value.

If it is desired to reduce the number of iterations the threshold may be set by (i) deducing an average amplitude of the noise in the transform signal within a window (ii) adding a margin to generate a threshold within that window: (iii) moving the window and repeating steps (i) and (ii).

In this way a threshold may be set which more closely follows the spectrum and thus is capable of more easily separating out unambiguous components of the transformed signal.

Where a number of iterations are performed, the missing or estimated data points may simply be replaced by that data estimated by the latest iteration. However in an alternative method the value of the estimated data points may be weighted by a proportion of the value estimated in the preceding iteration. For example the value used may comprise one and a half times the value of an estimated data point generated by a current iteration subtracted from which is a half times the value of that estimated in the preceding iteration. Although the invention has been described with reference to the generation of missing data points caused by channel failure in a multi channel radar receiving array, the same techniques can be used for a wide variety of purposes. For example an array could be designed with a smaller number of channels than would conventionally be required and a method according to the invention employed so as to create a "virtual" antenna array of larger size with a consequential saving in terms of hardware, albeit with a rise in the cost of processing. The advantage of creating new data points outside of the available data set is that the antenna aperture can be virtually enlarged, increasing the received bandwidth or the coherent integration time, all of which lead to enhanced resolution in the respective domain. An alternative use of the method is when a high amplitude noise burst needs to be removed in the time domain, inevitably with the underlying signal, and when the original data must be recreated or restored. For example in HF radar, high amplitude noise bursts during local thunderstorms caused by lightning strikes can corrupt the received radar echoes and, for a finite but possibly vital period, can render the radar blind. In one use of the method the input signal is derived from a radar receiver, and the missing data points correspond to known periods in the time domain for which the receiver radar echo was blinded by a lightning. In order that the invention may be well understood, an embodiment thereof will now be described by way of example with reference to the accompanying drawings, in which

Figure 1 is a diagram schematically illustrating a multi channel radar receiving array and showing how in azimuth signals arise in the frequency domain; i.e. the algorithm used for digital beam-forming is equivalent to a Discreate Fourier Transform.

Figure 2a to 2d are schematic time and frequency domain diagrams illustrating use of the method: Figure 3 is a flow chart showing the preiteration algorithm used for generating a threshold for use in a method according to the invention;

Figure 4 is a flow chart showing the method according to the invention; Figure 5 is a diagram illustrating the spectrum of an example input test signal:

Figure 6 is a diagram illustrating the spectrum of Figure 4, after a large proportion of the data points have been removed;

Figure 7 illustrates a spectrum of Figure 5 after a single iteration of the method;

Figure 8 illustrates the spectrum after 17 iterations: Figure 9 illustrates of a running window method for generating a threshold level; and

Figures 10 to 12 are time and frequency domain diagrams illustrating examples where the method may be employed. Referring to Figure 1 a plane wave 1 is incident at an angle θ onto a radar receiver 2 having a multiplicity of separate channels 3, some of which, e.g. 3_3, 3_n-5, have failed. It can be appreciated that the plane wave 1 strikes succeeding channels with a time or phase delay p. This phase delay is treated as a frequency, thus a spectral peak occurs in azimuth at an angle corresponding to the angle θ, the more acute the angle leading to a higher frequency signal. The effect of the finite number of data points, and indeed the missing data points, is to create uncertainty in the position of the main spectral peak s _θ leading to side lobes s, on either side.

The above mentioned case is much simplified. In a practical system each channel is used simultaneously to retrieve information for each of a number of overlapping finger beams in azimuth and also in range and doppler.

For use of the method the input data is a finite set of real or complex numbers from which some data points are missing. The positions of the missing data are known, as in the example of the failed channels in Figure 1. It is assumed that the input signal is composed from the sum of periodic functions, sine or cosine waves, of arbitrary frequencies. However it is a prerequisite that the number of components (sine waves) is significantly less than the number of data points in the set. The randomly distributed missing or zero value points will, in effect, modulate the input signal with noise (multiplicative noise). Figure 2A shows schematically an output signal with time which could be generated by the multi channel array 2 shown in Figure 1, except that targets striking the array from a variety of different azimuth angles give rise to a multiplicity of different frequency components. As an example, noise is generated in time slots corresponding to failed channels 3₃ and 3_n-5, As a first step the time domain signal shown in Figure 2A is Fourier transformed into the frequency domain as shown in Figure 2B. Here the sine waves appear as spectral lines whilst the missing data produce an incoherent noise like spectrum. It is necessary to weight the input data with a window function, in order to produce a lower side lobe level than the permitted noise after data restoration. It is important that the peaks of the noise are less than the largest spectrum line in the input so that at least one component of the signal can be filtered out. Thus referring to Figure 2B a threshold level T is set such that at least one spectrum component f₁ lies above that threshold value, the remaining spectral components being discarded. An inverse Fourier transform is then performed on the f₁ component leading to, for example, the sine wave shown in Figure 2C. From this an estimate for the value of the data signal at time slots corresponding to channel 3₃ and 3_{n -5} can be taken, the estimates now being referred to as 3₃ ^l and 3_n-5 ^l respectively. These values are then reinserted into the waveform shown in Figure 2A so as to generate an improved time domain signal, as shown in Figure 2D.

It should be emphasised that signals shown in Figures 2A to 2D are purely schematic. However Figure 2D is intended to show that in the "patched up" or updated set, the large signals are no longer modulated by the randomly distributed zeros and thus the overall noise power is reduced. Should further reduction of the noise be required, then the above steps can be iteratively repeated. Each iteration involves a re-estimation of the noise level and a corresponding adjustment of the threshold level. The noise decreases monitonically with the iterations, as more and more signal components contribute to the estimated values replacing the missing points and asymptotically approach a minimum. This noise floor is determined by parameters such as the ratio of the number of sine wave components to the number of data points, the external noise and the side lobe level of the window function used. Iterations may either be terminated after a fixed number of traversals or adaptively at that point where the noise no longer decreases in subsequent iterations.

Figures 3 and 4 show two flow charts illustrating the sequence of operations required in the method. Figure 3 shows the preiteration computation to estimate the first threshold level and to set up the arrays and pass parameters for the subsequent iterations. Firstly a window is applied to the input data. This weighted data is then kept in memory and updated in subsequent iterations After that a complex Fourier transform is computed using conventional radix 2 or a mixed radix FFT technique. If the number of input data points happens to be a prime number, then a slower Discrete Fourier Transform DFT algorithm must be used instead. The vector of the complex spectrum is kept in the memory for passing to the iteration sub programme. After that the magnitude of the spectrum is computed and normalised to its maximum This normalised vector of magnitudes is again a pass parameter to be updated during the iterations An estimated initial threshold (say between 0 and -10dB) can be used After that those components of the signal which exceed this initial value are removed and the mean of the remaining signal plus any noise is then computed. This mean is then multiplied with a margin (or a margin is added if logorithmic units are used) to obtain the first estimated threshold.

Figure 4 is a flow diagram illustrating the functions performed within a single loop. In the first stage the current threshold and normalised spectrum magnitude is used to generate a new data set. This contains only those complex spectrum lines whose magnitude exceeds the current threshold Below the threshold the signal plus noise is rejected, i.e. set to zero After that an inverse Fourier transform is performed on the chosen complex spectrum points to obtain an estimate of the input signal (when the noise below current threshold is excluded) From this estimate those points which are missing from the original input signal are selected and substituted in to the current input.

Once the input is ameliorated, the reduced noise level is estimated and a new threshold is set up using an identical method as in the preiteration step i .e. the new complex spectrum (CFFT) is computed. After that the old spectrum is updated or replaced. A normalised magnitude is then computed and this is then used update or replace the preceding normalised magnitude. The signal above the threshold value is rejected and the mean noise computed. After that a margin is applied to the mean noise to obtain the new threshold. The new threshold is then compared with the preceding threshold and if the ratio is greater than a given constant then another iteration is performed. If the ratio, or the measure of improvement is less than expected, then computation is stopped and the latest updated input data is output. As an alternative a fixed number of iterations for terminating the computation can be provided.

A more detailed program for implementing the method will now be described with reference to Figures 5 to 8. Figure 5 shows the spectrum of several sine waves of different frequencies and amplitudes. Missing data is simulated by multiplying some of the data points with zeros, which are distributed randomly within one batch of data or waveform period. As will be described below 125 data points are removed from 240 data points - i.e. a reduction of some 52%. As has been described previously the method involves the application of forward and inverse Fourier transform pairs and a thresholding operation between the two. The process for generating the sample test signal is described below.

The effect on the spectrum of removing some 52% of the data points is illustrated in Figure 6 where it can be seen that very few spectral peaks appear unambiguous, the remainder being buried within noise. The preiteration computations are described below.

Once transformed back into the frequency domain the spectrum has the appearance of that shown in Figure 7 where it can be seen that the level of the noise has been much reduced. A further threshold is then set and the process repeated. After seventeen similar iterations, reducing the threshold level each time, the spectrum has the appearance of that shown in Figure 8. Comparing Figure 8 to Figure 5 it can be seen that, in spite of 52% of the data points being missing, the spectrum is faithful to the input test signal. The details of the seventeenth iteration are described below.

In alternative variations of the method, the threshold value may be set predictively using a running window technique. Referring to Figure 9 a window W is moved across the spectrum. For each position of the window the average amplitude of the transformed signal is deduced and a margin added to generate a threshold T. By using this technique the threshold can more closely follow the contours of the spectrum and so separate more unambiguous spectral components during each iteration, thus reducing the number of iterations required.

As described, the missing or estimated data points are simply replaced by that data estimated by the latest iteration. However in an alternative method the value of the estimated data points may be weighted by a proportion of the value estimated in a preceding iteration. For example the value used may comprise one and a half times the value of an estimated data point generated by a current iteration subtracted from which is a one half times the value of that estimated in the preceding iteration. Although the invention has been described with reference to the generation of missing data points caused by channel failure in a multi channel radar receiving array, the same techniques can be used for a wide variety of purposes. For example an array could be designed with a smaller number of channels than would be conventionally required and a method according to the invention employed so as to create a "virtual" antenna array of larger size with a consequential saving in terms of hardware. The advantage of creating new data points outside of the available data set is that the antenna aperture can be virtually enlarged, increasing the received bandwidth or the coherent integration time, all of which lead to enhanced resolution in the respective domain. Another use is explained with reference to Figure 10 showing a high frequency spectrum extending between 3 and 30MHZ. Certain parts of the spectrum C are used for communication purposes by other users. Conventionally only a part of the spectrum located between a single pair of communications bands C can be used for radar purposes. However using a technique according to the invention all parts of the spectrum identified by R may be used for radar purposes, and that pan of the signal which would have been occupied by the communications bands C is generated using a method according to the invention. A high bandwidth signal allows high resolution, and that if returning signals occupying the bands C were simply ignored then this would lead to large range side lobes. An alternative use of the method is when a high amplitude noise burst needs to be removed in the time domain, inevitably with the underlying signal and when the original data must be recreated or restored. Figure 1 1 shows a received signal varying with time in an HF radar, a high amplitude noise burst B caused by a lightning strike can corrupt the received radar echo and for a finite but possibly vital period render the radar blind. The method according to the invention can be used to estimate the data which would have been received during this period of time. A similar technique may be used to restore scratched gramophone records, for example.

A mono-static radar operates by turning on a radar receiver, when the transmitter is turned off and vice versa. Figure 11 indicates the periods of time for which the receiver RX is switched on. A method according to the invention can be used to generate data which would have been received during the time that the transmitter TX is on, and the receiver RX is off, allowing the advantages previously described.

Claims

1. A method of estimating missing data points in an input signal which comprises the steps of.

b) setting a threshold so as to separate at least one unambiguous component of the transformed signal from the remaining components. c) inverse Fourier transforming the separated component to the one domain, and

d) replacing the missing data points in the input signal with data from the separated component to generate an improved signal.

2. A method according to claim 1 comprising iteratively repeating steps (a) to (d) on the improved signal except in step (b) setting a threshold lower than in the preceding iteration so as to separate a greater number of unambiguous components.

3. A method according to claim 2 comprising lowering the threshold in step (b) in each iteration until it is no longer possible to separate any further unambiguous components.

4. A method according to any preceding claim comprising setting the threshold in step (b) by deducing an average value of the amplitude of the transformed signal and adding a margin to that average value.

5. A method according to claim 4 in which the average value comprises the mean value.

6. A method according to claim 4 in which the average value comprises the median value.

7. A method according to any preceding claims in which the threshold in step (b) set by (i) deducing an average amplitude of the transformed signal within a window, (ii) adding a margin to generate a threshold within that window; and (iii) moving the window and repeating steps (i) and (ii).

8. A method according to any of claims 2 to 7 comprising weighting the value of the estimated data points by a proportion of the value estimated in the preceding iteration.

9. A method according to any preceding claim in which the input signal is derived from a multi channel radar receiving array, and in which the missing data points correspond to missing channels.

10. A method according to claim 9 in which the missing channels are deliberately omitted from the array to create a virtual array.

1 1. A method according to any of claims 1 to 8 in which the input signal is derived from a radar receiver, and in which the missing data points correspond to known periods in the time domain for which the received radar echo was blinded.

12. A method according to any of claims 1 to 8 in which the input signal is derived from a radar or communication receiver and in which the missing data points correspond to gaps in the frequency domain occupied by other users.

13. A method according to any of claims 1 to 8 in which the input signal is derived from a mono-static radar, in which reception is interrupted during transmission, and in which the missing data points correspond to the periods of time when the receiver is switched off.

14. A method substantially as described with reference to the drawings.