Coherent signal power detector using higher-order statistics

A signal detection system divides a data sampling run into blocks and perms a fast Fourier transform on each block, sorting results by frequency. Combinations of results of the fast Fourier transform corresponding to each frequency are processed to derive a test statistic which is unbiased by Gaussian noise while including such combinations of results of the fast Fourier transform which would be redundant over other combinations. Information concerning the frequency behavior of the signal derived in the course of detection, is accomplished with increased sensitivity.

DESCRIPTION 
BACKGROUND OF THE INVENTION 
1. Field of the Invention 
The present invention generally relates to the passive detection of signals 
which may be transmitted acoustically or electromagnetically to a 
detection location and, more particularly, to detection and spectrum 
analysis of signals in the presence of high levels of noise, such as may 
be encountered with hydrophones or in condition-based diagnostics of 
machinery as well as active detectors such as pulse Doppler radar. 
2. Description of the Prior Art 
Many arrangements for detection and enhancement of a signal in the presence 
of noise are well-known in the art; radio and television receivers being 
particularly familiar examples. Coherent radar is another application in 
which enhancement of a signal together with rejection of background noise 
and time-varying clutter is particularly critical. In all of these systems 
and others, however, some parameters of the signal to be detected, such as 
the carrier frequency of the signal of interest, are known. 
When this is not the case, such as in condition-based diagnostics where 
slight changes in the noise output of a complex machine may include one or 
more frequency characteristics which are specific to a potential 
malfunction, and it is desired to detect a completely unknown acoustic or 
electromagnetic signal in the presence of noise, the classical technique 
has employed spectrum analysis. In this technique, discussed in detail in 
"The use of Fast Fourier Transform for the Estimation of Power Spectra: A 
Method Based on Time averaging over Short, Modified Periodograms" by P. D. 
Welsh, IEEE Transactions on Audio Electroacoustics, Vol. AU-15, pp. 70-73, 
June, 1967, which is hereby fully incorporated by reference, a tonal or 
sine-wave signal in a receiver output having low signal to noise ratio 
(SNR) is analyzed at varying resolutions (e.g., bandwidths) to look for 
peaks which increase as bandwidth is reduced (e.g. as resolution is 
increased). If a signal is present having a frequency falling within one 
of a plurality of overlapping bandwidths, decrease of bandwidth will 
reduce the relative signal power attributable to noise, leaving 
substantially only the signal power in the tonal. 
This approach has three principal drawbacks, however. The methodology is 
inherently slow or hardware intensive since a reduction of bandwidth by 
any given factor increases the number of bands which must be processed by 
the same factor for any given resolution. Since the process is carried out 
sequentially at a plurality of resolutions, processing time or hardware 
(whether analog or digital) must be greatly multiplied as resolution 
increases. Further, since the signal of interest may vary in frequency due 
to modulation of frequency or Doppler effects (if either or both of the 
source or receiver are in motion), the signal power of the signal of 
interest may be distributed over several bands as resolution is increased. 
Additionally, to increase spectral resolution, known spectrum analysis 
techniques must increase the amount of data available which, in turn, 
increases processing time. 
Higher-Order Statistics (HOS) methods have been recently employed in 
coherent radar systems for object profiling and velocity measurement. 
Specifically, U.S. Pat. Nos. 5,227,801, and 5,231,403, to Robert D. Pierce 
which are hereby fully incorporated by reference, describe particularly 
effective techniques of achieving these goals. Even more recently, a 
technique of detecting moving and even accelerating targets with coherent 
radar and HOS techniques has been described in U.S. patent application 
Ser. No. 08/127,619, filed Sep. 28, 1993 now U.S. Pat. No. 5,402,131 
issued Mar. 28, 1995, (Navy case No. 75,280) by Robert D. Pierce, which is 
also fully incorporated herein by reference. HOS methods have several 
desirable characteristics of preserving phase information (e.g., 
coherency), are insensitive to linear phase shifts and suppress Gaussian 
noise effects. 
While HOS methods are used to enhance signals in radar systems, some 
aspects of the signals to be detected are necessarily known, as pointed 
out above. To date, no system has significantly increased the noise 
rejection of the classical technique of signal detection, described above, 
for an unknown signal. 
It is also generally recognized that a desired signal can often be detected 
under adverse conditions such as high noise levels if certain adaptations 
are made in the detector, based on the expected nature of the noise and/or 
the signal. However, just as small variations in the frequency of a signal 
can defeat or greatly reduce the effectiveness of known spectrum 
analysis-based detectors, a detector which is specifically adapted to 
certain signal characteristics is generally ineffective to detect signals 
having characteristics which are only slightly at variance therewith. 
Moreover, there is currently no consistent methodology or single apparatus 
capable of detecting a wider variety of signal variation than the spectrum 
analysis-based devices and methods described above which remains the 
methodology of choice where the characteristics of the signal to be 
detected are entirely unknown. Conversely, no technique has heretofore 
been developed to broaden the range of signal variation which will allow 
detection of a signal by a signal analyzer. 
SUMMARY OF THE INVENTION 
It is therefore an object of the present invention to provide a method and 
apparatus for increasing sensitivity and noise rejection in a signal power 
detector, such as a spectrum analyzer, or electromagnetic or acoustic 
receiver, such as a radio receiver or hydrophone system. 
It is another object of the invention to provide a signal spectrum analyzer 
having high speed of operation and capable of tracking a signal which 
varies in frequency over time. 
It is a further object of the invention to provide a signal enhancement 
method and apparatus which will provide rapid detection of an unknown 
signal. 
It is yet another object of the invention to provide a signal analyzer 
capable of detection of an unknown signal having a wide variation in 
frequency. 
It is another further object of the invention to provide a signal detector 
capable of providing information which characterizes the nature of 
frequency variation of the signal. 
In order to accomplish these and other objects of the invention, a signal 
detection method is provided including a computation of a test statistic 
from combinations of values in a block of signal samples including the 
steps of dividing a run of sampled signals into blocks, performing a fast 
Fourier transform on each block and sorting results of the fast Fourier 
transform by frequency, eliminating combinations of the results which are 
redundant, excluding combinations of results in which noise present in 
samples would bias the test statistic from the list of said combinations 
of samples, and determining detection of a signal based on a test 
statistic computed only from remaining combinations of results. 
In accordance with another aspect of the invention, an apparatus for 
detection of a signal in the presence of noise is provided comprising an 
arrangement for dividing a run of samples into blocks, a processor for 
performing a fast Fourier transform on each said block and sorting results 
of said fast Fourier transform by frequency, an index look-up table 
containing combinations of addresses for accessing combinations of 
transform results in which combinations of addresses which are redundant 
over other combinations of addresses in the index look-up table and 
combinations of addresses in which an addresses of a first pair of 
addresses was the same as one of the addresses of a second pair of 
addresses are excluded. 
In accordance with a further aspect of the invention, a spectrum anaylzer 
is provided including apparatus for detection of a signal in the presence 
of noise is provided comprising an arrangement for dividing a run of 
samples into blocks, a processor for performing a fast Fourier transform 
on each said block and sorting results of said fast Fourier transform by 
frequency, an index look-up table containing combinations of addresses for 
accessing combinations of transform results in which combinations of 
addresses which are redundant over other combinations of addresses in the 
index look-up table and combinations of addresses in which an addresses of 
a first pair of addresses was the same as one of the addresses of a second 
pair of addresses are excluded.

DETAILED DESCRIPTION OF A PREFERRED EMBODIMENT OF THE INVENTION 
Referring now to the drawings, and more particularly to FIG. 1, there is 
shown an exemplary environment in which the invention may be 
advantageously be employed. For example, as shown, a water-borne vessel 10 
carries a hydrophone transducer 12 for receiving acoustic signals 
transmitted through the water. In this example, the signal of interest is 
a sound made by a whale 14 as the target which is well-known to usually 
include target identifying parameters such as a rising or falling 
frequency tone (a time-frequency representation). The frequency of this 
tone may however be modified by relative motion of the water-borne vessel 
10 and the whale 14. Further in this noisy environment example, other 
background noises or signals derive from the noise of wind, rain 16 and 
breaking waves 18 as well as noise emanating from the hull and wake of the 
water-borne vessel and others. Each of these noise sources may have 
components having particular spectra such as white or pink noise or 
specific characteristic spectral components. Generally, however, signal 
sources of interest will have specific spectral components (e.g. periodic 
signals or tonals). 
It is to be understood that while the major meritorious effect of the 
present invention is in the detection of signal about which nothing or 
very little is known and deriving such information in the course of 
detection, the invention is also applicable to non-passive systems such as 
sonar and radar and other active detection processing. Particularly in the 
latter case, the invention provides a technique of detecting and analyzing 
signals having characteristics which would otherwise defeat pulse Doppler 
radar or other detection systems. 
In order to provide rejection of noise and signal enhancement, the 
invention provides a family of coherent signal detectors which are defined 
over even orders of products of signal samples. Odd order products can 
average to zero and their usefulness has not been fully evaluated. 
Further, the second order product is of very limited practical application 
since it does not provide detection for a varying frequency signal. The 
practical utility of eighth and higher order products has not been fully 
evaluated. 
Therefore, the following discussion will first describe the definition of a 
second order detector in order to convey an understanding of the 
evaluation of samples in accordance with the principles of the invention. 
Then these principles will be extended to fourth order and sixth order 
detectors together with practical implementations thereof; from which 
examples, extension of the principles of the invention and suitable 
implementing structure to higher order products will be evident to those 
skilled in the art as well as implementing structure for the second-order 
detector. 
In the preferred embodiment of the invention, a preprocessor is provided 
which takes a predetermined number (e.g., 32) blocks of samples, each 
containing a predetermined (e.g., 2048) samples. A fast Fourier transform 
(FFT) is performed on each block of samples and the results segregated 
into "bins" by frequency. Thus, each "bin" will contain a number of 
results equal to the number of blocks of samples corresponding to each of 
the blocks but at a single frequency. Since the data in each "bin" 
corresponds to a single frequency and the blocks of samples are separated 
in time, the contents of each bin are, in fact, time separated samples 
which are coherent in the sense of containing amplitude and phase 
information concerning any signal which may be present at that frequency. 
In contrast to classical spectrum analysis techniques which accumulate 
signal power, the present invention maintains the results of the FFT as 
discrete values which are sorted by frequency for further processing by 
HOS methods. If products of these discrete values or results are then 
averaged using HOS methods, any noise will average toward zero and any 
signal present will average toward a constant level, as is shown by the 
following analysis. 
The development of a coherent signal power detector in accordance with the 
invention begins similarly to the classical spectrum analyzer given by 
Welch. As shown in FIG. 2A, a signal, possibly containing a large 
proportion of noise is sampled over time and stored. The sampled time 
history, x((n-1).DELTA.t), with run length, T, is broken into P 
overlapping segments of N samples each. The sample period is .DELTA.T and 
the run length is 
EQU T=N.sub.T .DELTA.t (1) 
such that the total number of samples is 
EQU N.sub.T =N+(P-1)S (2) 
Each segment is shifted by S samples relative to the last segment, such 
that the pth segment is shifted by 
EQU .xi..sub.p =(p-1)S (3) 
samples. Each segment is then windowed and the Fast Fourier Transform (FFT) 
is taken of each segment as shown by the sequence of operations of FIG. 
2B. Where the transform is sampled at each frequency, .omega..sub.j, the 
result of the FFT operation can be expressed as 
##EQU1## 
where 
##EQU2## 
for positive frequencies 
##EQU3## 
such that 
##EQU4## 
At this point, known spectrum analysis techniques average and accumulate 
the magnitude squared of the P transformed segments at each frequency. 
Auto spectrum estimates of the power spectrum are formed by the average 
##EQU5## 
and the cross-spectrum estimate of the power spectrum is given by 
##EQU6## 
The factor, .LAMBDA..sub.K, is an appropriate normalization constant that 
depends on spectral resolution and the shape of the data window applied to 
the overlapping samples. The choice of data window is believed to be of 
substantially the same significance and function as in known spectrum 
analysis techniques in the practice of the invention. The criteria on 
which the choice of data window may be based is essentially unchanged from 
currently known spectrum analysis techniques and therefore need not be 
further discussed. For simplicity in the following discussion, a 
rectangular data window will be assumed. 
The signal of interest is assumed to be periodic but is permitted to change 
in frequency over time in a linear fashion. A linear change in frequency 
with time is often referred to in the art as a "chirp". This signal of 
interest is also assumed to exist in the presence of a significant amount 
of noise. Therefore, following the above analysis, the sampled time 
history is given by 
EQU x(t.sub.n)=A sin(.omega..sub.x t.sub.n +vt.sub.n.sup.2 
+.phi.)+.eta.(t.sub.n) (9) 
where 
EQU t.sub.n =(n-1).DELTA.t (10) 
In complex representation, Equation (9) may be expressed as 
##EQU7## 
The FFT of the pth segment for positive frequencies is 
##EQU8## 
Where .psi..sub.p (.omega..sub.j) in the transform of the noise segment. 
If a rectangular data window is assumed, for simplicity, 
EQU .kappa.(n.DELTA.t)=1,1.ltoreq.n.ltoreq.N (13) 
The blur term given by e.sup.iv(.DELTA.t(n-1)).spsp.2 can be ignored if the 
frequency varies by less than one or two cells or bins during one data 
segment or block (e.g. 2048 samples. Then, 
##EQU9## 
where, to simplify the expression, 
EQU .XI..sub.jp =.omega..sub.x +2v.xi..sub.p .DELTA.t-.omega..sub.j (15) 
The power spectrum is then 
##EQU10## 
and the cross spectrum is 
##EQU11## 
where the last term of the equation corresponds to the noise which is 
correlated between two sensors. Therefore, the ability to estimate the 
phase difference between two channels, as is presented by the direction of 
arrival of a signal at a sensor array, is impaired by the presence of 
noise which is correlated between the two channels 
In accordance with the invention, a second order product of transformed 
samples and conjugates is given by 
##EQU12## 
where the cross terms between the sinusoid and noise (which will average 
to zero) are ignored to simplify the expression. For linear variation of 
frequency over time, this product is non-stationary and it therefore does 
not provide a useful result for detection. However, for constant frequency 
(or velocity of the source or target), the product becomes 
##EQU13## 
This leads to the definition of the second order estimator 
##EQU14## 
subject to the constraints 
EQU b&gt;a (21) 
and 
EQU a.noteq.b (22) 
where the lag, l, is defined as 
EQU l.ident.a-b (23) 
and the number of averages for each lag is 
EQU P.sub.xx =P-l (24) 
The constraint that b must be greater than or equal to a avoids averaging 
of redundant products. 
The expected value of the estimator given by equation (20) is 
##EQU15## 
where biasing by Gaussian noise (including cross terms) is avoided by the 
further constraint of inequality (23) that a may not be equal to b. The 
number of averages, P.sub.xx, varies in dependence on the lag, l. 
In contrast to known spectrum analysis methods where correlated noise the 
estimation of delay or phase shift, this definition is easily extended to 
multiple channels, including cases where time delay or phase shift exists 
between to measurements, at differing locations, of the same signal or 
tonal. For example, 
##EQU16## 
While it can be seen from Equation (26) that phase information is preserved 
by the use of HOS analysis in accordance with the invention, the ability 
to extract phase information is limited because the phase steps in 
multiples of .omega..sub.x .DELTA.tS which is generally greater than 
2.pi.. Therefore, the phase information which can be directly extracted is 
ambiguous. However, the Fourier transform over data which is indexed by 
lag, l, could be used to extract the amplitude squared as long at the 
frequency of the signal of interest is substantially invariant over the 
observation time. Nevertheless, this is a significant restriction on the 
utility of a second order detector in accordance with the invention. The 
example, however, shows the potential for extracting information in a 
simple manner which is not available in known spectrum analysis 
techniques. 
The limitation on frequency variation during the sampling period is avoided 
in the family of fourth-order detectors which will now be described. 
Following the above analysis, the fourth order product of samples and 
conjugates is 
##EQU17## 
where the cross terms, which will average to zero, have again been ignored 
to simplify the expression. The fourth-order estimator is then given by 
##EQU18## 
subject to the constraint 
EQU d=a+b-c (29) 
that makes the fourth-order estimator independent of frequency, and the 
constraints 
EQU 0.ltoreq.a.ltoreq.P,a.ltoreq.b.ltoreq.P 
EQU 0.ltoreq.c.ltoreq.P,c.ltoreq.d.ltoreq.P 
EQU a.sup.2 +b.sup.2 &gt;c.sup.2 +d.sup.2 (30) 
that avoid averaging of redundant products, and the further constraint 
EQU a.noteq.c,b.noteq.c (31) 
that avoids biasing by Gaussian noise. The lag, l, is now defined as 
##EQU19## 
and the number of averages for each lag, P.sub.c, (the subscript 
indicating chirp) can be computed numerically. The mean or expected value 
of the fourth-order estimator is 
##EQU20## 
subject to the same constraints, as will be demonstrated below. 
It will be recognized by those skilled in the art, that the above 
constraints and definition of lag, l, correspond to equations (11), (13), 
(12) and (16) in the aforementioned U.S. Pat. No. 5,402,131 (Navy Case No 
75,280) and the implementation is substantially identical to that 
described therein. Constant frequency of the tonal corresponds to the 
constant velocity case and chirp corresponds to the constant acceleration 
case described therein. The following development of a description of this 
family of detectors generally parallels the development in this copending 
application but is placed in a more generalized notation. However, 
reference is made to the above-incorporated application for further 
details. 
When a constant frequency tonal is present in the jth bin, the expected 
value of the triple average taken over a specified region of values will 
have a real part that is proportional to the amplitude to the fourth power 
of the tonal. The expected value for those bins containing only noise is 
zero. Quantitatively, the expected value for a bin containing a tonal is 
##EQU21## 
and which is independent of lag, l. The non-zero part is real valued. For 
this case, the number of averages which are taken is 
##EQU22## 
Normalizing the fourth-order average by the power spectrum squared gives a 
type of coherency function. This method of normalization has added meaning 
for the constant frequency (or constant velocity) case because the real 
part of the fourth-order average is related to the signal power squared 
and the power spectrum is related to signal power plus noise power. For 
the constant frequency case, the test statistic is the fourth-order 
average divided by the power spectrum squared at each frequency. Therefore 
the ratio of these values is less than unity when noise is present. To 
demonstrate this, let the signal power in a spectral bin containing a 
signal of interest be 
##EQU23## 
and let the noise power in the spectral bin be 
EQU .sigma..sub..eta..sup.2 =G.sub..eta..eta. (.omega..sub.j) (37) 
Therefore, the signal to noise ratio (SNR) is 
##EQU24## 
Then, for moderate to high SNR, 
##EQU25## 
These results are easily extended to multiple channels where signal arrival 
angle introduces a time delay or phase shift between measurements of the 
same tonal at different locations in a sensor array. Three different 
averages are unique: 
##EQU26## 
Thus it is seen that Equation (42) gives only an amplitude and differs 
from Equation (41) by omission of a factor containing only a phase 
difference term. If the amplitude terms A and B are equal, Equation (40) 
differs from Equation (41) by twice this phase difference. Therefore, 
averages over multiple pairs of sensor outputs can readily be used to 
double the effective aperture size when estimating the phase difference. 
Moreover, correlated noise does not bias the result as in known spectrum 
analysis techniques. Therefore, improved estimates of direction of signal 
arrival can be made and achieved more simply than by known spectrum 
analysis techniques. By including a lag term, as described in the 
above-incorporated U.S. Patent Application, similar results can be 
obtained even when the tonal is varying linearly in frequency. 
The standard deviation of the fourth-order estimator normalized by signal 
plus noise power (power spectrum) squared is 
##EQU27## 
where the coefficients A, B, C and D are 
##EQU28## 
It is readily seen from FIG. 3 that the results of a simulation agree very 
well with this calculated value derived from equation (28). 
As further indication of the potential efficacy of the invention, the 
computed power spectrum squared as would be produced by conventional 
spectrum analysis techniques is shown in FIG. 4A and the spectrum detected 
by a fourth-order HOS detector in accordance with the invention is shown 
in FIG. 4B. Both of these Figures show a large peak representing tonal #1 
which is present at a 0.0 dB SNR. FIG. 4B clearly shows tonal #2 which is 
present at -3 dB SNR far above any other peaks which represent noise. This 
result is in sharp contrast with tonal #2 in FIG. 4A which remains buried 
in noise and undifferentiated therefrom. 
To implement a fourth-order signal power detector in accordance with the 
invention, reference is now made to FIGS. 5-7. A preprocessor 500 is shown 
in FIG. 5 which is used to convert the time sampled data into a plurality 
of overlapping blocks of data, perform a fast Fourier transform on each 
block and store the output of the FFT for each block. Specifically, a 
time-history containing a signal of interest and noise is sampled at a 
large number of time instants depending on the desired spectral resolution 
and the amount of overlap as indicated in Equation (2), above. In the 
preferred embodiment of the invention and to limit the sampling period, 
2048 samples are provided in each of 32 blocks. The amount of overlap is 
somewhat arbitrary and can be adjusted to allow for the available data. 
Overlaps greater than 50% are generally regarded in the art to give no 
significant advantage in comparison to processing time required. 
Therefore, for the preferred embodiment of preprocessor 500, register 510 
will require not more than 65,536 (e.g. 64K) digitized samples at the 
desired accuracy or resolution (e.g. number of significant bits). 
It should be understood in regard to FIGS. 5-7 that while a special purpose 
processor could be provided and may provide some performance enhancement, 
it is presently preferred to emulate the functional elements shown in 
FIGS. 5, 6, 6A and 7 with a suitably programmed general purpose processor, 
preferably including a math coprocessor in view of the ease of doing so, 
as is well understood in the art, and the flexibility for modification 
provided thereby. In this regard, it should be noted that certain types of 
noise, such as sea clutter (e.g., Doppler processing) are only 
approximated by Gaussian noise functions and there is much debate among 
those skilled in the art concerning processing enhancements, some of which 
will be discussed below, which will allow rejection of a portion of the 
small amount of noise remaining. Therefore, the flexibility for adding or 
removing some processing steps afforded by the utilization of a general 
purpose processor is preferred. 
The samples in register 510 are broken into blocks by two interconnected 
counters 520 and 530. Block counter 520 counts the number of blocks which 
have been generated at any given time, providing a number to multiplier 
525 which multiplies the number of the previous block by the number of 
samples in the block less the number of samples of overlap and supplies 
the result as an initial count to sample counter 530. From this initial 
count, sample counter 530 is incremented a number of times equal to the 
number (e.g. 2048) samples which will comprise each block. When this 
number of counts is reached by sample counter 530, an output is provided 
to block counter 520 causing the block count to be incremented and a new 
initial count to be provided to sample counter 530. By repeating the 
process in this way, blocks of sample data are built up in registers 540. 
As each of these registers 540 is filled, the contents are provided to a 
fast Fourier transform processor 560. The individual data samples may be 
weighted, if desired, to impose a window on the data as schematically 
indicated at 550, prior to performance of the fast Fourier transform 
operation. The results of this operation are separated by frequency in 
accordance with this operation and the results stored in one of registers 
570.sub.1 -570.sub.32, under control of the block counter 520. It should 
be understood that a section across corresponding bins of registers 
570.sub.1 -570.sub.32 corresponds to register 110 of FIG. 6 and register 
110 need not be separately provided. It should also be noted that if the 
signal in register 510 is derived from an active detection system or 
otherwise should contain complex values (e.g. I and Q samples from a pulse 
Doppler radar system) each of registers 570.sub.1 -570.sub.32 would be 
twice as long since the Fourier transform of complex data signals does not 
necessarily have conjugate symmetry around zero frequency, as will be 
understood by those skilled in the art. 
Referring now to FIG. 6, the major storage elements of implementation 100 
as shown in FIGS. 6 and 7 include a sample register 110, a data sample 
look-up table 118 and an index look-up table (LUT) 120. Data storage space 
is set aside in high-speed Random Access Memory (RAM) accessible by the 
processor during initialization. Then, index look-up table 120 is created. 
It is considered important to the preferred implementation of the invention 
to create the index LUT 120 in such a way as to eliminate all data which 
does not meet the constraints of Equations (29)-(32), above. Specifically, 
for P=32 complex data in a block, this reduces the amount of data stored 
and the consequent number of computations of averages from 32,768 to 2600, 
in accordance with equation (35). Additionally, it is both hardware 
efficient and computationally efficient to compute and store lag indices, 
l, during creation of the index LUT. Therefore, even though for a constant 
P, the addresses provided into the data sample LUT and lag indices for 
various combinations of samples would always be the same, the possibility 
of change of P and the size of the index LUT's which would result make it 
preferable to create the index LUT anew each time the processing in 
accordance with the invention is performed in view of the simplicity of 
doing so in combination with computation of l for each address. 
Specifically, the index LUT is preferably created by four nested loops with 
exit conditions for a, b or c being equal to P and branching conditions 
which discard combinations of a, b and c which do not meet the conditions 
specified in equations (29), (30) and (31), above. The preferred procedure 
for creating the index LUT requires only a relatively few FORTRAN.TM. 
commands and may be expressed as follows, wherein npts=P: 
EQU npas=0 
Main loop - loop over all combinations (states) of ja, jb, jc and jd that 
will be tested. 
EQU do 20 ja=1,npts 
Restrict jb to be equal to or greater than ja to remove redundancy in 
accordance with equation (30). 
EQU do 20jb=ja, npts 
EQU do 20jc=1, npts 
Restrict jd to be equal to or greater than jc to remove redundancy in 
accordance with equation (30). 
EQU do 20jd=jc,npts 
Then apply equation (29) to remove combinations of ja, jb, jc and jd which 
do not satisfy it. Of course, jd could be computed directly from ja, jb 
and jc. Then discard combinations of samples which would cause sensitivity 
of the average by source frequency. 
EQU jta=ja+jb-jc-jd 
EQU if(jta.eq.0) go to 30 
EQU go to 40 
EQU 30 continue 
Then, combinations are excluded which would cause biasing of the average by 
Gaussian noise as given by Equation (31). 
EQU if (ja.eq.jc) go to 40 
EQU if (jb.eq.jc) go to 40 
Then, the lag is calculated, as follows, in accordance with equation (32) 
and is restricted to greater than or equal to 0 in accordance with the 
last condition of equations (30) to remove redundant states. 
EQU jtb=ja**2+jb**2-jc**s-jd**2 
EQU if(jtb.lt.0) go to 40 
Then l is computed (as jtb2) using jtb, above, and a maximum is iteratively 
found. 
EQU jtb2=jtb/2 
EQU if(jtb2.gt.imax) imax=jtb2 
The index LUT 120 is completed by storage as an address to the index LUT 
(npas) varies from 1 to 2600. 
EQU npas=npas+1 
EQU jja(npas)=ja 
EQU jjb(npas)=jb 
EQU jjc(npas)=jc 
EQU jjd(npas)=jd 
EQU jab(npas)=jtb2 
The process is exited when looping is completed to ensure that the index 
LUT will not exceed required limits and the number of possible states 
which satisfy all conditions of equations (29), (30) and (31) are counted 
and stored to the index LUT. The number of states which produce the same l 
are also counted and stored for future use in constant change of frequency 
with time (e.g. constant chirp or constant acceleration) processing. The 
return loop commands allowing discarding of combinations is: 
##EQU29## 
and the main loop return is provided by: 
##EQU30## 
At this point, the index LUT 120 will be complete to provide access to the 
data sample LUT 118. A block of P (=32) data samples of the complex values 
of the returned signals can now be input to the data sample register 110 
(FIG. 6). Then, all allowed combinations of the data samples are selected 
by multiplexer 112 and multiplied by multiplier 114 and the products 
placed in the data sample look-up table 118 in accordance with pairs of 
addresses ranging from 1 to P. 
To compute the averages, a sequencer 130, which effectively functions as a 
clock which is synchronized to the processing, provides a logical 0 or 1 
to dual selector 132, 132' such that, at any instant, one pair of values 
of a and b or c and d will be provided to the data sample look-up table. 
Of course, two such tables could be provided to allow a single access of 
each. However, it is considered more efficient to provide only a single 
table from which a single product is retrieved at a time since the data 
sample LUT is large and stored data would be redundant (unless the 
conjugates of the products are pre-computed for storage in a second 
table). However, since some products will not be used, it is more 
efficient to compute the conjugate later and only for those products 
actually retrieved. 
The accessed product is provided to a demultiplexer 140 which is also 
driven by sequencer 130 and thus synchronized with dual selector 132, 
132'. If selector 132 provides addresses to the data sample LUT, the 
retrieved product is directly stored to register 142. If selector 132' 
provides addresses to the data sample LUT, the conjugate of the product is 
computed at 141 and the result stored to register 142'. The values in 
registers 142 and 142' are then multiplied by multiplier 144 and the 
result accumulated. 
In order to accommodate processing for both constant frequency and constant 
change of frequency with time (e.g. constant chirp), it is preferable to 
accumulate the result of multiplication in two different ways 
simultaneously. For the constant frequency case where lags are, by 
definition, the same, accumulation at accumulator 146 is without regard to 
lag, l. At the conclusion of the process, the real part of the accumulated 
result is divided by 2600 and normalized by the power spectrum squared as 
the HOS average for the constant frequency case. For the constant chirp 
(or constant source or target acceleration) case, accumulation is done by 
lags with the calculated value of l, stored in the index LUT and divided 
by the number of values accumulated for each value of l for each set of 
data to be retrieved, in sequence, with l used as an address in 
accumulator array 150 in which the accumulation is to be done. 
Referring now to FIG. 6A in which accumulator array 150 is repeated for 
reference, the accumulated values corresponding to values of l from 1 to 
64 are mapped into a register 160. Each value is directly mapped into a 
stage of register 160 from 65-128 and the conjugate is computed and stored 
in a mirror-symmetrical location in stages 1 to 64. As pointed out above, 
there will be no accumulated value for a lag of zero or for certain other 
values of l, such as 62 and 63. These values are supplied by interpolation 
by interpolators 152, 154 and 156. Any interpolation technique is suitable 
for this purpose, and the details of the interpolation function performed 
are unimportant to the practice of the invention. 
Then, if desired (e.g. to de-emphasize the values toward the ends of the 
register), each value can be multiplied by a value to impose a window on 
the data in register 160 by multiplier window 170. It should be noted that 
the application of a window of a particular shape is common practice in 
the art to remove side lobes and improve resolution. In the invention, at 
the present time, a rectangular data window is used to achieve a slight 
improvement in performance. However, a cosine window can be used on the 
theory that the density of interpolated values (for which there is no l 
produced by a combination of a, b, and c) will increase near the ends of 
register 160. The cosine window results in a slight loss of performance 
over a constant weighting (e.g. rectangular window) but has value when 
resolving the chirp between detected targets. That is, a rectangular 
window should be preferable for detection and a tapered window such as a 
cosine window should be preferable for resolving a plurality of sources or 
targets within the same frequency bin. 
The symmetrical mapping of values and conjugates into register 160, as 
shown, causes an artifact to be developed by performance of a fast Fourier 
transform 180 which can be eliminated after the sign of each value is 
inverted at 190 after the fast Fourier transform is performed at 180. The 
results then correspond to chirps or accelerations, sorted into "bins" 
from which a maximum real part can be selected or the real parts of the 
values in the "bins", normalized by the power spectrum squared, displayed 
as a chirp or acceleration spectrum at 195. The display of a spectrum is 
preferred, as discussed above, particularly where plural sources/targets 
having different chirps/accelerations may be present in the same frequency 
bin. 
Then to complete the process, the result of either the constant frequency 
case and/or the constant chirp case which has been divided by the power 
spectrum squared and then compared to a threshold, as shown in FIG. 7, to 
determine whether detection has occurred. In practice, the threshold is 
empirically adjusted (e.g., manually), or in accordance with predicted 
performance as will be discussed with reference to FIG. 8, until the false 
alarm rate falls to an acceptable level such as 10.sup.-3. However, noise 
which may be encountered in practical use of the invention may be variable 
or modulated by some natural mechanism as a function of frequency. It has 
been found by the inventor that some dynamic adjustment of the test 
statistic prior to thresholding is particularly advantageous in 
combination with the present invention to provide additional sensitivity 
and reliability of detection. Utilization of such dynamic adjustment is 
assumed in the foregoing discussion. This dynamic adjustment is referred 
to as constant false alarm rate (CFAR). However, CFAR is considered to be 
a perfecting feature of the invention and is not necessary to the practice 
thereof. CFAR will now be discussed with reference to FIG. 7. 
As explained in the above-incorporated U.S. Patent Application, if each 
value in the sample input register 110 is multiplied by its conjugate and 
the results summed and then divided by the number of samples, P (=32), an 
estimate of the power spectrum of the "signal-plus-noise" (or noise 
variance over a bandwidth associated with a bin if a target is not 
present) of the data samples (or more properly, an estimate of the signal 
power (if any) plus noise variance) can be derived at each frequency bin. 
In summary, the power spectrum squared at a frequency bin has been divided 
by the HOS estimator to derive a ratio that is used as the test statistic 
value. This will adjust the test statistic to a fixed threshold such that 
the dynamic variations in the noise statistic maintain the false alarm 
rate to quite closely track a constant value, on average. 
Since operation of the detector in accordance with the invention at a 
substantially constant false alarm rate is thus made possible, it is 
preferable to establish detection thresholds graphically from the 
predicted false alarm rate as shown in FIG. 8. FIG. 8 shows the predicted 
false alarm rates of both the constant frequency (810) and the linearly 
varying (constant chirp) frequency (820) cases, discussed above, derived 
from simulations. If the desired sensitivity is to be a false alarm rate 
of one false alarm per thousand tries, then the threshold is determined by 
constructing a line horizontally at a false alarm rate of 10.sup.-3 which 
yields the SNR thresholds of -4.5 dB for the constant frequency (or 
velocity) case and -1.5 dB for the linearly varying frequency (constant 
acceleration or chirp) case. If the threshold comparator of FIG. 7 is set 
to provide a positive indication when the actual SNR exceeds one of these 
values (e.g. -4.5 dB, a high degree of certainty of detection is implied 
by a positive output of the threshold comparator, as shown in FIG. 9. 
In comparison to the power spectrum or quadratic detector of FIG. 10, at a 
10.sup.-3 false alarm rate, it is seen that the invention has similar 
performance to a quadratic detector for the linearly varying frequency 
case. However, it must be recalled that this level of performance of the 
quadratic detector requires knowledge of the characteristics of the noise 
whereas the invention delivers a comparable level of signal detection 
performance without requiring such knowledge. Therefore, the invention 
provides a capability heretofore unavailable in signal power detectors at 
very slight cost in performance (about 7% probability of detection at the 
SNR of the set threshold, as shown in FIGS. 9 and 11). Further, the 
invention is able to characterize the nature of the chirp and quantify the 
rate of change of frequency of the detected signal by virtue of processing 
of data indexed by lags in the same manner that chirp or acceleration 
could be quantified in the above-incorporated U.S. Patent Application. 
The invention is also capable of applying a consistent methodology to a 
wide variety of signal characteristics and variations thereof by the 
simple expedient of altering the spectral resolution (e.g. by adjustment 
of the sampling period and/or the FFT size corresponding to the number of 
samples included in a block of samples) which will alter the frequencies 
allocated to each bin when the FFT is performed as shown in FIG. 5. For 
example, the sampling rate and block size can be adjusted to keep the 
frequency variation to less than two bins, as indicated by Equations (1), 
(5) and (6). Overlap is not critical and has been found to be optimally 
about 50% and can be adjusted somewhat to fit the available data without 
significant loss of performance. That is, while reduction of overlap 
corresponds to discarding of data performance losses are negligible for 
most overlap adjustments which may be required. Overlap adjustment is 
readily accomplished by changing the number by which the block number is 
multiplied at 525 of FIG. 5. By the same token, and approximately 50% 
overlap will allow a 64K sample memory to accommodate 32 blocks of 4096 
samples, each, corresponding to a practical upper limit on the spectral 
resolution which may be required since increased spectral resolution 
limits the amount of frequency variation over time which will permit 
detection. 
To demonstrate the efficacy of the present invention, FIGS. 8 and 9 shows 
the receiver operating characteristics of a fourth order detector in 
accordance with the invention for both constant and linearly varying 
frequency. In FIG. 9, these characteristics are approximately the same, as 
indicated by only slight variation at 910 and 920, for a substantial 
variation of SNR. In other words, for a SNR threshold (the horizontal axis 
of FIG. 9 which also is the same scale as the threshold SNR in dB 
indicated at the right side of FIG. 4B) set to the actual SNR (represented 
by the curves 910, 920 and 930 of FIG. 9), the probability that detection 
being made by the invention is only slightly less than 40%. If a constant 
frequency were to be detected and the threshold set in accordance with a 
-4.5 dB SNR as shown in FIG. 8 to correspond to a 10.sup.-3 false alarm 
rate, a signal having a -2 dB SNR would be detected about 85% of the time. 
Similarly, by interpolation, a signal having a 0 dB SNR would be detected 
in virtually every instance. 
The performance of a sixth-order detector in accordance with the invention 
is similar to that indicated above for the fourth-order detector. However, 
the definition and preferred implementation will be illustrative of 
variations of the invention which are deemed to fall within the scope 
thereof. Essentially, some loss of potential performance is unavoidable as 
the order of the detector is increased since the number of constraints to 
avoid redundancy and biasing of the processing result by Gaussian noise 
reduces the number of products which can be averaged. On the other hand, 
this reduces storage requirements and provides some simplification of the 
processing required, particularly if a special purpose processor is 
provided. Higher order processing can also provide additional information 
concerning frequency variation over time and can also allow extension of 
the methodology to numbers of channels greater than two. 
Specifically, consider the sixth-order product 
##EQU31## 
where the cross terms between the sinusoid and noise are omitted for the 
same reasons as in the fourth-order embodiment. This leads to the 
definition of the sixth-order estimator as 
##EQU32## 
subject to the constraint that makes the sixth-order estimator independent 
of frequency 
EQU f=a+b+c-d-e (50) 
and the constraint that makes the sixth-order estimator independent of 
linear frequency sweep 
EQU a.sup.2 +b.sup.2 +c.sup.2 =d.sup.2 +e.sup.2 +f.sup.2 (51) 
and the constraint that avoids the averaging of redundant products 
EQU 0.ltoreq.a.ltoreq.P,a.ltoreq.b.ltoreq.P,b.ltoreq.c.ltoreq.P 
EQU 0.ltoreq.d.ltoreq.P,d.ltoreq.e.ltoreq.P,e.ltoreq.f.ltoreq.P 
EQU a.sup.3 +b.sup.3 +c.sup.3 .gtoreq.d.sup.3 +e.sup.3 +f.sup.3 (52) 
and the constraint that avoids biasing by Gaussian noise 
EQU a.noteq.d,b.noteq.d,c.noteq.d (53) 
If an accelerating frequency is now allowed, then a lag, l, would be 
defined as 
##EQU33## 
The number of averages P.sub.s, is computed numerically and for P=32, 
P.sub.s =1166. By comparison, it will be recalled that P.sub.c was 2600 in 
the constant frequency case of the fourth-order detector. The expected 
value of this estimator is 
##EQU34## 
subject to the above constraints. It can be appreciated that this 
estimator, allows simultaneous consideration of a larger number of 
combinations of sensor outputs and sensor array geometries from which 
phase information could be extracted and, for example, can be analyzed to 
provide quantification of numerous types of motion of the source of the 
signal of interest. Analysis of up to six channels are possible with a 
sixth-order detector. 
A practical implementation of a sixth-order detector is shown in FIG. 
6.sub.1. This implementation differs from the fourth-order embodiment 
shown in FIG. 6 by the omission of storage of precomputed products of 
pairs of data. The reduced number of products which must be computed after 
constraints are applied is reduced to levels at which the preferred 
computation arrangement of FIG. 6 does not engender an equal degree of 
comparative convenience. 
More specifically, look up table 120' now provides six addresses in 
register 110 (which is effectively the same structure as a bin from each 
one of registers 570 of FIG. 5) to selector 112' to provide six values to 
register 143. The conjugates of the latter three of these values is 
computed at conjugate calculators 141'. These values are then multiplied 
at multiplier 144' and the results simultaneously accumulated for constant 
frequency detection at 146' and the individual results accumulated in 
accordance with particular values of lag, l, in lag table 150'. These 
accumulated values may then be indexed by lag and further processed by the 
arrangement of FIG. 6A precisely as for the fourth-order embodiment of the 
invention. 
These methods and practical implementations thereof may be readily extended 
to higher order averages for sinusoids that have phase terms containing 
higher order polynomials. For example, in an eighth-order embodiment, 
constant frequency in the fourth-order embodiment would correspond to a 
first order polynomial in phase; a linear chirp to a second order 
polynomial in phase, and so on for higher order embodiments. These 
techniques can thus be used to decompose polynomial phase signals to 
determine the coefficients of these polynomials. However, as the order of 
the detector increases, a higher SNR is expected to be required to 
correctly estimate the coefficients. 
It should be appreciated from the foregoing that during each data segment, 
the frequency of the signal of interest can vary up to two bin widths 
without significant distortion of its FFT. Such distortion will take the 
form of a frequency shifts for frequency variations of less than two bin 
widths but will have little effect on the shape of the spectrum. However, 
as the transform at each bin is averaged in accordance with the invention, 
the signal of interest will have spent only part of the time in or near 
that bin. Therefore, it is likely that a higher SNR will be required for 
isolating a signal of interest where the frequency change is wide or rapid 
since a wider processing bandwidth is required to contain the variation. 
The invention also provides for this possibility by allowing control of 
the sample data into blocks or segments at multiplier 525 and counter 530. 
For example, the sample run could be divided into two separate runs and 
each processed at half the spectral resolution by reducing the size of the 
FFT (e.g., to 512 samples). The results from each run can then be 
coherently averaged or displayed as a function of time. If this is not 
sufficient to capture the frequency variation, the sample run could be 
divided into four parts and the spectral resolution divided by four and so 
on. However, it must be recognized that reducing the spectral resolution 
increases the bandwidth of each frequency bin and increases the noise 
power. Therefore, it is anticipated that an increased SNR may be necessary 
for frequency variations of greater size or rate. Nevertheless, the 
reduction of Gaussian noise provided by the selection of values for 
calculating the averages provides detection capability of unknown signals 
by a consistent methodology which was not heretofore available. 
In view of the foregoing, it can be seen that the invention provides a 
signal detector of improved performance and sensitivity at low signal to 
noise ratios and capable of rapid signal analysis for detection, even when 
the signal varies in frequency over time. The invention is capable of 
detecting signals about which nothing is known in advance of the analysis 
and results in the production of substantial information about the signal 
which can be used to advantage in other signal analysis techniques, 
including classical spectrum analysis by allowing the portion of the 
spectrum containing the signal of interest to be isolated for detailed 
analysis. Therefore, in view of the simple and consistent methodology of 
the invention and implementations thereof, it is anticipated that the 
invention will provide features, when included in signal analyzers of 
currently known types, which will enhance user interaction therewith as 
well as adding the capability of detecting and tracking an unknown signal. 
For example, user specified limitation of the portion of the spectrum to 
be analyzed is often provided in commercially available spectrum 
analyzers. If the present invention is combined therewith in a single 
instrument or test arrangement, this feature could be used to limit 
processing in accordance with the invention to frequency bins of interest 
to increase speed of detection. The detector in accordance with the 
invention is completely applicable to any periodic signal, including radar 
and sonar signals which may be transmitted and reflected from objects as 
well as passive detection of signals, regardless of source or origin, 
including condition-based diagnostics of machinery and passive tracking of 
the motion of signal sources with hydrophones and the like. 
While the invention has been described in terms of a single preferred 
embodiment, those skilled in the art will recognize that the invention can 
be practiced with modification within the spirit and scope of the appended 
claims.