Deciding whether a received signal is a signal of interest

[0001] The present invention generally pertains to signal communications and is particularly directed to deciding whether an observed received signal, such as a radio communication signal, is a signal of interest (SOI).

[0002] Various different attributes of an observed received signal may be processed to decide whether the observed received signal is a signal of interest.

[0003] The resources that may be used for making such a decision pertain to data processing and computing, receiver gain, the frequency to which the receiver is tuned and the timing of when the receiver is being tuned to a particular frequency. In some circumstances the time constraints for making such a decision are such that there is a need to optimize the use of such resources while maintaining the probability of accurately making the decision.

[0004] When one or more of these resources is limited, it may not be practical to continuously observe the full spectrum of interest; and it may not be practical to continuously process all of the various different attributes of the observed received signal that are typically processed in order to decide whether the observed received signal is a signal of interest. When a network of sensors (such as antenna elements) for observing attributes of the received signal is used, it may not be possible to communicate all of what is observed from each sensor to a computing center, because of the bandwidth limitations of the coupling with the sensors.

[0005] The present invention provides a method of deciding whether an observed received signal is a particular signal of interest (SOI), comprising the steps of:

[0006] (a) providing a statistical model of the particular signal of interest and a statistical model of the environment when the particular signal of interest is not present; and

[0007] (b) with a computer, processing the observed received signal in accordance with a decision rule that uses the provided statistical models to decide whether the observed received signal is the particular signal of interest;

[0008] wherein the decision rule is:

[0009] $\begin{array}{cc}\mathrm{decision}=\{\begin{array}{cc}\mathrm{SOI}& \mathrm{if}\mathrm{LLR}\left(y\right)>+\lambda \\ ~\mathrm{SOI}& \mathrm{if}\mathrm{LLR}\left(y\right)<-\lambda \\ \mathrm{unsure}& \mathrm{otherwise}\end{array}\mathrm{where}\mathrm{LLR}\left(y\right)=\log \frac{P\left(a_1\left(y\right)·a_2\left(y\right),\dots |\mathrm{SOI}\right)}{P\left(a_1\left(y\right)·a_2\left(y\right)\dots |~\mathrm{SOI}\right)},& \left[\mathrm{Eq}.1\right]\end{array}$
where LLR is log likelihood ratio, P denotes a probability, a₁(y) and a₂(y) are functions of first and second attributes of the observed signal y in a sequential order that said attribute functions are presented for computation of LLR(y), and +λ and −λ are real numbers,

[0010] wherein the statistical model of the particular signal of interest includes functions of attributes that individually have a probability distribution function f₁ when the observed received signal is the particular signal of interest such that
KL(f₁,f₀)>0,  [Eq. 2]
where f₀ is a probability distribution function (PDF) of the respective individual attribute when the observed received signal is not a signal of interest.

[0011] In one aspect of the method of the present invention, at least one of the attribute functions is expressed as:
a(y_C,y_R)f_c=f_c
a(y_C,y_R)Δf=Δf
a(y_C,y_R)t_c=t_c
a(y_C,y_R)Δt=Δt,
where y_C is a sample of the observed signal received by a commutated element of a commutated antenna array, y_R is a sample of the observed signal received by a commutated element of the antenna array, f is frequency, and t is time.

[0012] In another aspect of the method of the present invention, at least one of the attribute functions is expressed as one of the following:

[0013] $\begin{array}{cc}{a\left(y_C,y_R\right)}_{\frac{S}{N}}=\frac{{\langle y_C,y_R\rangle }^2}{{y_C}^2{y_R}^2-{\langle y_C,y_R\rangle }^2},& \left(i\right)\end{array}$
where y_c is a sample of the observed signal received by a commutated element of a commutated antenna array and y_R is a sample of the observed signal received by a commutated element of the antenna array;
a(y_C,y_R)AOA=arg max_ΦColor*U(U*diag(Energy)U)⁻¹U*Color
where
U_n=u_kn(Φ)
Color_n=<y_C,y_R>_Rn
Energy_n=<y_R,y_R>_Rn  (ii)
y_C is a sample of the observed signal received by a commutated element of a commutated antenna array and y_R is a sample of the observed signal received by a commutated element of the antenna array;

[0014] $\begin{array}{cc}{a\left(y,k\right)}_{\mathrm{Sinusoids}}=\frac{{Y\left(k\right)w}^2}{{y\left(k\right)-Y\left(k\right)w}^2}\mathrm{where}w={\left(Y\left(k\right)*Y\left(k\right)\right)}^{-1}Y\left(k\right)*y\left(k\right)& \left(\mathrm{iii}\right)\end{array}$
and k is the k in k-means;
a(y,k)_{PSK k-ary 1}=a(y^k,k)Sinusoids,  (iv)
and k is the k in k-means;

[0015] $\begin{array}{cc}{a\left(y,k\right)}_{\mathrm{PSK}k-\mathrm{ary}2}=\frac{\sum _i{\mathrm{Mean}\left(C_i\right)}^2}{\sum _i\mathrm{Variance}\left(C_i\right)}\mathrm{where}\left\{C_i\right\}=k-\mathrm{means}\left(\Delta y\right)& \left(v\right)\end{array}$
and k is the k in k-means;

[0016] $\begin{array}{cc}{a\left(y,k\right)}_{\mathrm{PSK}k\text{-}\mathrm{ary}3}={\max }_{\tilde{T}\in \left[\frac{1}{2\Delta f},\frac{2}{\Delta f}\right]}\frac{\sum _i{\mathrm{Mean}\left(C_i\right)}^2}{\sum _i\mathrm{Variance}\left(C_i\right)}\mathrm{where}\left\{C_i\right\}=k\text{-}\mathrm{means}\left({\bigcup }_i{s}_i^{\prime }\right),& \left(\mathrm{vi}\right)\end{array}$
and k is the k in k-means;

[0017] $\begin{array}{cc}{a\left(y,k\right)}_{\mathrm{QAM}k\text{-}\mathrm{ary}1}=\frac{\sum _i{\mathrm{Mean}\left(C_i\right)}^2}{\sum _i\mathrm{Variance}\left(C_i\right)}\mathrm{where}\left\{C_i\right\}=k\text{-}\mathrm{means}\left(y^2\right),& \left(\mathrm{vii}\right)\end{array}$
and k is the k in k-means;

[0018] $\begin{array}{cc}{a\left(y,k\right)}_{\mathrm{QAM}k\text{-}\mathrm{ary}2}={\max }_{\tilde{T}\in \left[\frac{1}{2\Delta f},\frac{2}{\Delta f}\right]}\frac{\sum _i{\mathrm{Mean}\left(C_i\right)}^2}{\sum _i\mathrm{Variance}\left(C_i\right)}\mathrm{where}\left\{C_i\right\}=k\text{-}\mathrm{means}\left({\bigcup }_i{s}_i^{\prime }\right),& \left(\mathrm{viii}\right)\end{array}$

[0019] k is the k in k-means, f is the frequency of the signal, and U denotes a union operator;

[0020] $\begin{array}{cc}{a\left(y\right)}_{\mathrm{AM}}=\frac{{<V\left({\tilde{\omega }}_c,y,\omega \right),V\left({\tilde{\omega }}_c,y,-\omega \right)>}^2}{\begin{array}{c}{V\left({\tilde{\omega }}_c,y,\omega \right)}^2{V\left({\tilde{\omega }}_c,y,-\omega \right)}^2-\\ {<V\left({\tilde{\omega }}_c,y,\omega \right),V\left({\tilde{\omega }}_c,y,-\omega \right)>}^2\end{array}}{\tilde{\omega }}_c=\arg {\max }_{{\omega }_c}{<V\left({\omega }_c,y,\omega \right),V\left({\omega }_c,y,-\omega \right)>}^2& \left(\mathrm{ix}\right)\end{array}$
and ω is frequequncy.

[0021] The following is an explanation of the meaning of the term “function of an attribute of the received signal”.

[0022] Let x(M) denote a transmitted signal where M is the transmitted message. Let y=C(x(M)) denote an associated received signal where C denotes the effect of the channel. Let x(M)∈H⊂^N for some set H where H denotes some “class” of signals. For example, H may be the set of all 4-ary PSK or AM modulated signals, in which case H is a proper subset of ^N. Let y∈^N. x(M) and y are random variables. Let a:^N→ denote an attribute function and a(y) an attribute of y. a(y) is used to decide whether or not x(M)∈H. Let f₁ denote the PDF of a(y) when x(M)∈H, and f₀ denote the PDF of a(y) when x(M)∉H. For an attribute to be useful it is necessary that f₁ and f₀ are dissimilar. A very important measure of similarity is the Kullback-Leibler divergence defined as,

[0023] $\begin{array}{cc}\mathrm{KL}\left(f_1,f_0\right)=\int f1\left(u\right)\log \frac{f_1\left(u\right)}{f_0\left(u\right)}du& \left[\mathrm{Eq}.3\right]\end{array}$
Notionally, a(y) is “useful” if
KL(f₁,f₀)>0.  [Eq. 4]
(KL(f₁,f₀) is always non-negative.)

[0024] The present invention additionally provides non-transitory computer readable storage media that include computer executable program instructions for causing one or more computers to perform and/or enable one or more of the steps of the above-described methods. Each storage medium may include a plurality of non-transitory media.

[0025] Additional features of the present invention are described with reference to the detailed description.

[0037] Referring to FIG. 1, an exemplary system in which the method of the present invention is performed includes an array of commutated antenna elements 10, a reference antenna element 11, a receiver 12, a computer 14 and a computer display monitor 16. In some embodiments the receiver 12 and the computer 14 are combined in a single assembly. The computer 14 contains at least one digital signal processor and non-transitory computer readable storage medium/media that include computer executable program instructions for causing the computer to perform and/or enable the various processing steps that are described herein. These instructions are stored in the computer readable storage medium/media of the computer when the computer is manufactured and/or upon being downloaded via the Internet or from a portable non-transitory computer readable storage medium/media containing such instructions.

[0038] Samples 20 of an observed signal received by the array of commutated antenna elements 10 are obtained by the receiver 12 and samples 22 of the observed received signal are obtained by the receiver 12 from the reference antenna element 11. Preferably, the samples 20, 22 of the observed received signal are observed in zero-mean additive white Gaussian noise (AWGN).

[0039] Individual samples 20 are obtained from different elements of the array of commutated antenna elements 10. The samples 20, 22 are obtained at different times, at different frequencies and/or by using different CDMA signal access codes. The usefulness of attributes of the received signal increases with the number of samples because the probability distributions functions (PDFs} become more localized with one caveat. From an information-theoretic point-of-view the Sampling Theorem must be respected. More samples that are not statistically independent are not constructive.

[0040] The samples 20, 22 of the observed received signal are processed by the computer 14 to decide whether the observed received signal is a particular signal of interest.

[0041] Also referring to FIG. 2, a statistical model 24 of the particular signal of interest and a statistical model 25 of the environment when the particular signal of interest is not present are provided, as shown in FIG. 2. The statistical model of the particular signal of interest includes functions of attributes that individually have a probability distribution function f₁ when the observed received signal is the particular signal of interest such that
KL(f₁,f₀)>0,  [Eq. 2]
where f₀ is a probability distribution function of the respective individual attribute when the observed received signal is not a signal of interest.

[0042] The computer 14 processes the observed received signal samples 20, 22 in accordance with a decision rule that uses the provided statistical models 24, 25 to decide whether the observed received signal is the particular signal of interest, as shown at 26. The decision rule is:

[0043] $\begin{array}{cc}\mathrm{decision}=\{\begin{array}{cc}\mathrm{SOI}& \mathrm{if}\mathrm{LLR}\left(y\right)>+\lambda \\ ~\mathrm{SOI}& \mathrm{if}\mathrm{LLR}\left(y\right)<-\lambda \\ \mathrm{unsure}& \mathrm{otherwise}\end{array}\mathrm{where}\mathrm{LLR}\left(y\right)=\log \frac{P\left(a_1\left(y\right)·a_2\left(y\right),\dots |\mathrm{SOI}\right)}{P\left(a_1\left(y\right)·a_2\left(y\right)\dots |~\mathrm{SOI}\right)},& \left[\mathrm{Eq}.1\right]\end{array}$
where LLR is log likelihood ratio, P denotes a probability, a₁(y) and a₂(y) are functions of first and second attributes of the observed signal y in a sequential order that said attribute functions are presented for computation of LLR(y), and +λ and −λ are real numbers.

[0044] P(a₁(y),a₂(y), . . . |SOI) is a model for SOI; and P(a₁(y),a₂(y), . . . |˜SOI) is a model for ˜SOI. For example, when y is a SOI, one would expect LLR(y) to be very large and positive, and specifically, LLR(y)>>λ. Opposite expectations apply when y is ˜SOI.

[0045] The decision rule defined by Equation 1 is ternary. The possible decisions are SOI, ˜SOI, and “unsure”.

[0046] Referring to FIG. 3, when it is determined, as shown at 27, that processing by the computer 14 of a particular observed received signal in accordance with a decision rule that includes the statistical model for the particular SOI has resulted in a false decision either that the observed received signal is a signal of interest, as shown at 28, or that the observed signal is not a signal of interest, as shown at 30, the operator causes the computer 14 to construct a revised statistical model for the particular observed signal, as shown at 32, and/or to revise at least one of the parameters +λ and −λ of the decision rule that is used in the decision processing step for the particular observed signal, as shown at 34. A determination that such processing has resulted in a false decision can be made by a system operator or by some agent external to the system.

[0047] Preferably, the decision processing step 26 is incremental in that LLR(y) is computed incrementally for cumulatively ordered sampled attributes a₁(y),a₂(y), . . . until it is decided whether or not the observed received signal is a signal of interest, as shown in FIG. 4.

[0048] The computer 14 first computes

[0049] ${\mathrm{LLR}\left(y\right)}_1=\log \frac{P\left(a_1\left(y\right)|\mathrm{SOI}\right)}{P\left(a_1\left(y\right)|~\mathrm{SOI}\right)},$
as shown at 36. When LLR(y)₁>+λ or >−λ, it is decided that the received signal is either SOI or ˜SOI and the decision processing step is complete, as shown at 38. However, when |LLR(y)₁|<λ and >−λ, an “unsure” decision has been made, and the decision processing step continues by next computing

[0050] ${\mathrm{LLR}\left(y\right)}_2=\log \frac{P\left(a_1\left(y\right),a_2\left(y\right)|\mathrm{SOI}\right)}{P\left(a_1\left(y\right),a_2\left(y\right)|~\mathrm{SOI}\right)},$
as shown at 40.

[0051] When |LLR(y)_n| is either >λ or <−λ, a decision as to whether or not the received signal is a SOI has been made and the decision processing step is complete.

[0052] It is possible that after some finite number of computations of LLR(y)_n, the decision may remain “unsure”. Nonetheless, the decision processing step will be completed with a decision of SOI or ˜SOI provided that the attributes functions are defined reasonably and are presented for computation of LLR(y) in an order that minimizes the average time for making a decision as to whether or not the observed received signal is a signal of interest. When the received signal y is a SOI, LLR(y)_n usually increases “almost” monotonically.

[0053] Preferably, the attributes functions a₁(y),a₂(y), . . . are presented for computation of LLR(y) in an order for incremental computation that minimizes the average time for making a decision as to whether or not the observed received signal is a signal of interest. This criterion derives from a classical construction in statistics, the Bayesian Hypothesis Test, wherein it is desired to make a decision among a candidate set of hypotheses. Here the hypotheses are SOI and ˜SOI. (The ternary and incremental decision is an extension.) For the Bayesian Hypothesis Test a defined cost function is specified. There is a cost for deciding SOI when the correct hypothesis is ˜SOI, and there is a cost for deciding ˜SOI when the correct hypothesis is SOI. The Bayesian Hypothesis Test construction minimizes the expected cost. It is also necessary to specify the a priori probabilities of each hypothesis. The Bayesian Hypothesis Test is described at http://en.wikipedia.org/wiki/Bayesian_inference.

[0054] For example, one may select the attribute functions a(y)_Energy and a(y_C,y_R)_fc and detect a target signal by observing energy over a threshold within the frequency band of a SOI.

Initial Identification of a Target Signal as a SOI

[0055] Some approaches for identifying a target signal as a signal of interest SOI for use establishing criteria that are used for constructing the statistical models include the following:

[0056] First, after making several passes of the environment while assuming that no SOI is present, a statistical model P(a₁(y),a₂(y), . . . |˜SOI) of the environment is constructed and stored in a SOI database. When an unusual signal y is detected the unusual signal is processed in relation to the statistical model of the environment to determine whether the unusual signal y is an outlier.

[0057] Note that if one were to use only a(y)_Energy as an attribute function, ordinary energy detection may be used with a constant false-alarm rate. Here λ has the slightly different interpretation of testing whether the PDF of the attribute is sufficiently small to classify the unusual signal y is an outlier. Note also that it is not sufficient to use a cumulative distribution function CDF of the attribute, as such a technique would fail for a PDF of an attribute that is concentrated at multiple locations.

[0058] Secondly, a system operator may specify characteristics of a SOI. For example, this may be something like, a(y)_Energy>τ and a_AOA(y_C,y_R)∈(−15°,+15°). A reasonable approach is to regard this as a statement about a PDF. For example, assume that a(y)_Energy is uniformly distributed in some bounded range greater than τ and a_AOA(y_C,y_R) is uniformly distributed in the interval (−15°,+15°). The detection process allows these PDF's to be refined once the operator affirms a labeled target as a SOI.

[0059] Once the operator affirms a target signal as a SOI, a statistical model P(a_i(y)|SOI) of the target signal is entered into the SOI database as a particular signal of interest. A SOI database entry for a particular signal of interest may be revised upon subsequent identifications of the particular target signal. This calculation is more expensive and preferably is performed in the background. After that, the earlier discussion applies.

Estimation of +λ

[0060] An initial value of +λ is estimated by initially labeling target signals in a graphical user interface (GUI) as a SOI in accordance with the estimated value of +λ. The initial value of +λ may be estimated prior to deployment of the system in the field. A target signal that is labeled as a SOI is entered into a SOI database. When it is subsequently determined by a system operator that a particular target signal which has been labeled as a SOI for the estimated +λ is ˜SOI, pursuant to feedback from the operator, the particular target signal is removed from the SOI database, and the value of +λ is increased; whereby in the future, it may be more difficult to decide if a received signal is a SOI.

[0061] The operator also may provide feedback for the other error type by indicating an instance when a particular target signal is an SOI after the particular target signal had been implicitly labeled as ˜SOI by not being labeled as SOI in accordance with the estimated value of +λ.

[0062] The estimation of +λ is made in relation to a lower value region of +λ in which it is desired that the decision be “unsure”. If one wants to make it more difficult for the decision to be SOI, one should change the value of +λ to increase the region in which the decision is “unsure”. One may not want to make a corresponding change to the region in which the decision is ˜SOI; whereupon +λ and −λ would not be positive and negative values of the same real number.

[0063] For the Bayesian Hypothesis Test the optimal +λ depends on the cost function and the a priori probabilities of the hypotheses.

Determining that Attributes are Useful Attributes

[0064] Whenever a system operator affirms that a target signal is an SOI, there is an opportunity to determine which attributes of the target signal are useful attributes by estimating the probability distribution function (PDF) of the each of the individual attributes. Such estimates are made by using kernel density estimation, which is described at http://en.wikipedia.org/wiki/Kernel_density_estimation. For kernel density estimation, the PDF is approximated as a linear combination of translates of a kernel function. The kernel function integrates to 1. For a linear combination of N kernel functions dividing by N ensures that the integral of the linear combination integrates to 1, as required for a PDF. The bandwidth parameter of the kernel function stretches the kernel function along its domain. A smaller bandwidth concentrates the kernel function in a smaller region of its domain. Hence smaller bandwidth produces a linear combination that varies more rapidly. Conversely, larger bandwidth produces a linear combination that varies more smoothly.

[0065] The preferred kernel function is the PDF of a Normal distribution. This is probably the most common function used for kernel density estimation. There is no particular significance to the Normal distribution, other than as a PDF, this kernel integrates to 1, as required.

[0066] The bandwidth parameter is determined by a cross-validation technique. For a modest number of samples S, the samples are partitioned as S=S₁∪S₂. For a given bandwidth, one estimates the distribution by using S₁ and then assesses the efficacy of the asserted bandwidth by computing the likelihood of S₂.

[0067] The statistical model P(a₁(y),a₂(y), . . . |SOI) is a joint distribution. To construct a reasonable estimate of a distribution, there must be enough samples. For a single attribute distributed along a line, this is simple enough. But for even a joint distribution with two attributes, it may be difficult to get enough samples in a plane. This is the so-called curse of dimensionality. It only gets worse with many attributes since the decision rule uses a product of PDF's of the individual attributes. This may be regarded as a joint PDF where it is assumed that the attributes a₁(y),a₂(y), . . . are statistically independent. This has been found to be a generally satisfactory assumption, albeit not exactly correct.

Criteria for Screening Attributes

[0068] When screening attributes that may be considered for representation in the statistical model, an attribute m is preferred over an attribute n for construction of the statistical model of a particular signal of interest when C(n){P(m)−1}>C(m){P(n)−1}. C(i) is the time cost of computing attribute i and P(i) is the probability of satisfying the associated criteria.

Presentation of Attribute Functions in an Optimal Order for Computation of LLR(y)

[0069] As discussed above, the decision processing step is incremental in that LLR(y) is computed incrementally for cumulatively ordered sampled attributes functions a₁(y), a₂(y), . . . until it is decided whether or not the observed received signal is a signal of interest; and it is preferable that the attribute functions are presented for computation of LLR(y) in an order that minimizes the average time for making a decision as to whether or not the observed received signal is a signal of interest.

[0070] In deciding whether to present a particular attribute function a₁(y) for computation before another attribute a₂(y), one should determine whether attribute a₁(y) is likely to increase LLR(y)₁ more rapidly than a₂(y) and/or the relative time costs of observing and computing a₁(y) and a₂(y). The expected increase in LLR(y)₁ associated with presenting a₁(y) is the quantity KL(P(a₁|SOI), P(a₁|˜SOI)). The time cost of observing and computing a_i(y) is determined in accordance with:

[0071] $\begin{array}{cc}i=\arg {\max }_i\frac{\mathrm{KL}\left(P\left(a_i|\mathrm{SOI}\right),P\left(a_i|~\mathrm{SOI}\right)\right)}{C_i}& \left[\mathrm{Eq}.5\right]\end{array}$

[0072] An attribute function a(y) that is based upon observations of a particular attribute over different time intervals may be treated in the statistical model as two different attribute functions. For example, an AOA estimate computed over 10 ms takes less observational computation time than the same AOA estimate computed over 20 ms, but may be less useful because it may have a higher variance.

[0073] Another consideration is an assumption of statistical independence. When independence of attributes is assumed the PDF of an attribute is an estimate of the marginal PDF with respect to other attributes. Essentially, it is an average over the other attributes. In particular, it is marginal with respect to S/N. The instantaneous S/N may differ considerably from the average and it is appropriate to revise the estimate of the rate of change of LLR(y) for a particular attribute selection once that attribute has been selected. This is very similar to the exploitation/exploration paradigm in reinforcement learning.

Reinforcement Learning

[0074] In some exemplary embodiments, reinforcement learning is used to allocate computer resources for performing the decision processing step 26.

[0075] Reinforcement learning is an area of machine learning in computer science, concerned with how an agent ought to take actions in an environment so as to maximize some notion of cumulative reward. Reinforcement learning involves trial-and-error and feedback to arrive at the cumulative reward. Examples of cumulative rewards that are applicable to the present method of determining whether a particular observed signal is a signal of interest are: (1) minimizing the time cost of computation, such as by minimizing the average time for making a decision as to whether or not the observed received signal is a signal of interest; (2) minimizing the time cost of observing the attributes; and (3) minimizing the time cost of communicating the observed attributes.

[0076] Referring to FIG. 5, during the construction of a statistical model for a particular signal of interest, the attribute functions a₁(y),a₂(y), . . . are so ordered, as shown at 42, that during the decision processing step 26 the attribute functions a₁(y),a₂(y), . . . are presented for computation of LLR(y) in an order that minimizes the average time for making a decision as to whether or not the observed received signal is a signal of interest.

[0077] The average time to make a decision that completes the decision processing step 26 is estimated by using the formula

[0078] $\frac{\mathrm{KL}\left(f_1,f_0\right)}{\mathrm{Time}\mathrm{Cost}},$
as shown at 44. The estimate of the average time to make a decision is revised in accordance with cumulative estimates, as shown at 46.

[0079] In some exemplary embodiments, reinforcement learning is used to allocate spectral tuning resources and/or temporal tuning resources and/or receiver gain resources for the statistical model construction step.

[0080] The time cost of communicating observed attributes is relevant to a situation in which both the time of arrival and the angle of arrival of the observed received signal are communicated from a sensor network. Computing the time of arrival requires samples on the order of the bandwidth (i.e. the Nyquist Theorem). This may involve a lot of data; and the time cost of sending this data may be very large.

[0081] A useful reinforcement learning method is the so-called N-armed bandit paradigm, which is described at http://en.wikipedia.org/wiki/Multi-armed_bandit.

Distinguishing Different Observed Received Signals of Interest From One Another

[0082] Referring to FIG. 6, in an exemplary embodiment, different observed received signals of interest are distinguished from one another by using the method of the present invention. Statistical models that are different for a plurality of different particular signals of interest received signals are separately provided, as shown at 50; and the decision processing step 26 (FIG. 2) is performed, as shown at 52, by separately using the provided statistical models of the different particular signals of interest. The different observed received signals are then distinguished from one another in accordance with the differences in the decisions made when separately using the provided statistical models for the different particular signals of interest, as shown at 54.

[0083] Referring to FIG. 7, in an exemplary embodiment, a statistical model is constructed for the union of a particular set of observed signals and other observed signals by using attributes associated with selected criteria for distinguishing the particular set of observed signals from the other observed signals, as shown at 56.

Attribute Functions

[0084] The following is a discussion of how some attribute functions may be used in the statistical models in exemplary embodiments of the present invention for different types of attributes. An attribute is a function of an instance of an observed signal.

Energy or S/N

[0085] For any rectangle one may compute,

[0086] $\begin{array}{cc}{a\left(y_C,y_R\right)}_{\underset{N}{S}}=\frac{{\langle y_C·y_R\rangle }^2}{{y_C}^2{y_R}^2-{\langle y_C·y_R\rangle }^2}& \left[\mathrm{Eq}.6\right]\end{array}$
where <. . . > denotes inner-product and y_c and y_R denote samples from the commutated and reference elements respectively. From earlier material and here stated without proof, this is an estimate of

[0087] $\frac{S}{N}$
for 2 antenna elements in AWGN and a single coordinate descent step. The denominator is always non-negative. As the energy of the noise approaches zero, y_C→sy_R for some complex scalar s and a(y_C,Y_R)→∞.
Ordinary Energy
a(y)_Energy=∥y∥²  [Eq. 7]

[0088] Notice that instead of a(y_C,y_R) this attribute function has been written as a(y). This highlights the fact that this attribute function calculation does not require both the commutated and reference element samples, and either y=y_R or y=y_C may be used. Or, one could use both for a slight improvement. Under the assumption of AWGN for the channel, the distribution of a(y) is easy to understand. a(y) has a chi-square distribution for x(M)∉H. If we may also assume x(M)∈H has constant energy, then a(y) is non-central chi-square for x(M)∈H. Typical f₁ and f₀ are shown in FIGS. 9A and 9B, where one can see that a(y) is useful as KL(f₁,f₀)>0.

Spectral and Temporal Support.

[0089] The support of a function refers to the subset of its domain where the function is non-zero. For the rectangular decomposition of a received signal, as shown in FIG. 8, the spectral support is characterized by f_c and Δf, and the temporal support is characterized by t_c and Δt.

[0090] In some embodiments, energy and f_c could be the only attributes used to identify an SOI. This is clearly unnecessarily restrictive.

For the present invention one may use any of the following attribute functions to identify a SOI:
a(y_C,y_R)f_c=f_c
a(y_C,y_R)Δf=Δf
a(y_C,y_R)t_c=t_c
a(y_C,y_R)Δt=Δt,  [Eq. 8]
As before, one may also use a(y) rather than a(y_C,y_R).
Here x(M)∈H for some SOI, and otherwise x(M)∉H for other signals that are not of interest. Generally, one would expect a(y_C,y_R) to have a broader distribution for x(M)∉H (, as characterized by f₀). For example, in some concepts of operation, the notion is that entirety of the probability mass of f₁ is located near f_c.
Angle of Arrival (AOA)

[0091] In the rectangular decomposition of a received signal, as shown in FIG. 8, for any rectangle R the color of the rectangle Color(R) is defined as the complex scalar,
Color(R)=<y_C,y_R>_R.  [Eq. 9]

[0092] The subscript R on the inner-product indicates that the computation of the inner-product is restricted by the spectral and temporal support of the rectangle R. What is rendered in FIG. 8 is Arg{Color(R)}. It has been previously demonstrated that the color sequence in the rectangular decomposition is a sufficient statistic for computing an MLE (most likely estimation) of the AOA. Accordingly, one may define an attribute function as a point estimate of the AOA in accordance with,
a(y_C,y_R)_AOA=arg max_ΦColor*U(U*diag(Energy)U)⁻¹U*Color
where
U_n=u_kn(Φ)
Color_n=<y_C,y_R>_Rn
Energy_n=<y_R,y_R>_Rn  [Eq. 10]
u_kn(Φ) is taken from the calibration tables for the antenna, the response of antenna element k_n o a signal from angle Φ. Note that this is a multiple angle calculation and to compute a single angle Φ is restricted to a single tuple. As in the case of the above discussed spectral and temporal support attributes, x(M)∈H for some SOI, and otherwise x(M)∉H for other signals that are not of interest. Generally, one would expect a(y_c,y_R) to have a broader distribution for x(M)∉H (as characterized by f₀). A related attribute is the a posteriori PDF of the AOA. Here the range of a(y_C,y_R) is the set of all non-negative functions that integrate to 1 defined on the interval [−π,+π].

[0093] For the next group of attributes functions, specific structural assumptions are made about x(M). These attribute functions include: (a) multiple sinusoids to include CW and FSK, (b) PSK, (c) FM and QAM, (d) AM and DSB, and (e) SSB. For these attribute functions, the channel is more problematic for several reasons. First, frequency error is unconstrained. Secondly, digital modulation samples obtained away from the center of a symbol may not exhibit the characteristics of the constellation.

CW and FSK.

[0094] Spectrum estimation is a classic signal processing problem and there are many different approaches, each with advantages and disadvantages.

Suppose x(M)_n=s e^iωn where s is a complex scalar that models amplitude and phase and ω is the frequency of the sinusoid. Of course, s and ω are unknown. Let y=x(M)+v where v is AWGN.

One approach is to compute the MLE for ω and s and then compute the associated

[0095] $\frac{S}{N}.$
Here,
ω_MLE=arg max_ω|ℑ(y,ω)|²  [Eq. 11]
s_MLE=ℑ(y,ω_MLE).  [Eq. 12]
ℑ denotes the Fourier transform.
With

[0096] $\begin{array}{cc}{\tilde{x}}_n=s_{\mathrm{MLE}}{e}^{i{\omega }_{\mathrm{MLE}}n},\frac{S}{N}=\frac{{\tilde{x}}^2}{{y-\tilde{x}}^2}.& \left[\mathrm{Eq}.13\right]\end{array}$

[0097] While this method enjoys certain statistical optimality and this is a simple calculation for a single sinusoid, the product space of the parameters quickly becomes problematic for multiple sinusoids.

[0098] One alternative method is AR spectral modeling. Here, y(n) is approximated as a linear combination of immediately preceding samples. One earlier sample is required for each sinusoid. For example, for one sinusoid y(n)=ω₁y(n−1) where ω₁=e¹⁰⁷ . In general, y=Yw for an appropriately defined Y. Generally, w is estimated by ordinary LS as, w=(Y*Y)⁻¹Y*y and the computational complexity of this algorithm is determined by the computation of w. For a small number of sinusoids Cramer's rule may be used. For larger problems notionally Gaussian elimination may be used, but recognizing that Y*Y is Toeplitz, the Levinson-Durbin algorithm is preferred for O(N²) complexity.

[0099] $\frac{S}{N}=\frac{{\mathrm{Yw}}^2}{{y-\mathrm{Yw}}^2}.$
In AR spectral modeling the estimated frequencies of the sinusoids are determined by the locations of the zeros of the FIR filter determined by w.

[0100] Using an AR spectral model,

[0101] $\begin{array}{cc}{a\left(y,k\right)}_{\mathrm{Sinusoids}}=\frac{{Y\left(k\right)w}^2}{{y\left(k\right)-Y\left(k\right)w}^2}\mathrm{where}w={\left(Y\left(k\right)*Y\left(k\right)\right)}^{-1}Y\left(k\right)*y\left(k\right)& \left[\mathrm{Eq}.14\right]\end{array}$
Here k denotes the order of the AR model and the number of sinusoids. k is the number of tuples of w. For example, for CW k=1 and for binary FSK k=2.
Notice that because of modulation, binary FSK is not a sum of two pure sinusoids. An AR model is essentially an FIR filter, and symbol transitions induce a transient of length k−1. This degrades

[0102] $\frac{S}{N}$
and a more elegant solution would be to detect these transients (by anomalous energy) and delete these samples from the calculation of w. This produces much better behavior for a(y,k), but also destroys the Toeplitz structure of Y*Y and thereby significantly increases computational difficulty.
PSK

[0103] Consider first the idealized case where the receiver samples only at the center of a transmit symbol. For k-ary PSK,

[0104] ${x\left(M\right)}_n=s{e}^{iwn}{e}^{i\frac{2\pi }{k}{m}_n}$
where s is an unknown complex gain and ω is an unknown frequency error. m_n∈ is the unknown transmit data sequence. Let y=x(M)+v where v is AWGN.

[0105] A first approach is to observe that {x(M)_n}^k=s′e^ikωn and absent noise (y_n)^k={x(M)_n}^k. Let y^k denote this vector. y^k is a sinusoid at frequency kω and it is reasonable to define,
a(y,k)_{PSK k-ary 1}=a(y^k,k)_Sinusoids  [Eq. 15]

[0106] A second approach is to observe that

[0107] ${x\left(M\right)}_n\stackrel{\_}{{x\left(M\right)}_{n-1}}=s^{\prime }{e}^{iw}{e}^{i\frac{2\pi }{k}{m}^{\prime }}$
and absent noise

[0108] $y_n\stackrel{\_}{y_{n-1}}={x\left(M\right)}_n\stackrel{\_}{{x\left(M\right)}_{n-1}}.$
Let

[0109] ${\left(\Delta y\right)}_n=y_n\stackrel{\_}{y_{n-1}}.$
Δy is the same k-ary constellation rotated by ω.

[0110] The efficacy of a 4-ary PSK model may be determined by applying k-means (with k=4) and computing

[0111] $\frac{S}{N}.$
Let {C_i}=k-means(Δy) where {C_i} is a set of clusters which represent a partition of the points of Δy. For any single cluster C we may define

[0112] $\frac{S}{N}=\frac{{\mathrm{Mean}\left(C\right)}^2}{\mathrm{Variance}\left(C\right)}.$
Similarly, for the set of all clusters {C_i} we may define

[0113] $\frac{S}{N}=\frac{\sum _i{\mathrm{Mean}\left(C_i\right)}^2}{\sum _i\mathrm{Variance}\left(C_i\right)}.$
Accordingly,

[0114] $\begin{array}{cc}{a\left(y,k\right)}_{\mathrm{PSK}k-\mathrm{ary}2}=\frac{\sum _i{\mathrm{Mean}\left(C_i\right)}^2}{\sum _i\mathrm{Variance}\left(C_i\right)}\mathrm{where}\left\{C_i\right\}=k-\mathrm{means}\left(\Delta y\right)& \left[\mathrm{Eq}.16\right]\end{array}$

[0115] For the purpose of implementation, it is desirable to seed each k-means trial with k points that are uniformly spaced on a circle with random rotation. The result is independent of the radius of the circle used for this seeding.

[0116] Now consider what happens when samples are not necessarily taken at the center of transmit symbols. Suppose that the center of the first transmit symbol is at time t₀ and the center of the n^th symbol is at time t₀+(n−1)T. There are two issues here. The receiver does not know t₀ or T. Generally, the difficulty here is to determine an appropriate set of sampling instants, by a reasonable computation. The following method has been found to satisfy this requirement.

[0117] Begin by fixing an estimate of T. Let {tilde over (T)} be an estimate of T. Consider

[0118] $\tilde{T}\in \left[\frac{1}{2\tilde{\Delta }f}·\frac{2}{\tilde{\Delta }f}\right]$
where {tilde over (Δ)}f=a(y_C,y_R)_Δf. Let

[0119] $S_i={\bigcup }_{t\in \left[i\tilde{T}·\left(i+1\right)\tilde{T}\right]}y\left(t+\tilde{T}\right)\stackrel{\_}{y\left(t\right)}.$
The density with which samples are taken in a {tilde over (T)} interval may vary, but generally four samples are satisfactory. For T={tilde over (T)} there is at least one sample of the form

[0120] $s^{\prime }{e}^{\mathrm{i\omega }}{e}^{i\frac{2\pi }{k}{m}^{\prime }}.$
Other samples will be in a neighborhood of this value and generally there will be no bias in the error because of source compression (i.e. adjacent symbols are independent).

[0121] Let {C_i}=k-means(∪_iS_i) and μ_i=Mean(C_i). For the cluster {C_i} determine a partition wherein ∪_iS_i and μ_i are the cluster centers. FIG. 10 is an exemplary histogram of clusters.

[0122] Let s′_i=arg min_u∈Si|μ_i−u|². That is, s′_i is that element of S_i that is closest to μ_i. Let μ′_i=k-means(∪_is′_i). This process may be iterated. It appears that typically five iterations are sufficient to determine the correct subset of the samples. That is, to determine the centers of the transmit symbols. Thus,

[0123] $\begin{array}{cc}{a\left(y,k\right)}_{\mathrm{PSK}k-\mathrm{ary}3}={\max }_{\tilde{T}\in \left[\frac{1}{2\tilde{\Delta }f}·\frac{2}{\tilde{\Delta }f}\right]}\frac{\sum _i{\mathrm{Mean}\left(C_i\right)}^2}{\sum _i\mathrm{Variance}\left(C_i\right)}\mathrm{where}\left\{C_i\right\}=k-\mathrm{means}\left({\bigcup }_i{s}_i^{\prime }\right)& \left[\mathrm{Eq}.17\right]\end{array}$

[0124] It is still appropriate to retain a(y,k)_{PSK k-ary 2} as defined previously, for instances where the knowledge of symbol timing is present.

FM and QAM.

[0125] The treatment of FM and QAM is similar to the above discussed treatment of PSK. Consider first the idealized case where the receiver samples only at the center of a transmit symbol. For QAM, x(M)_n=s e^iωnu(m_n) where s is an unknown complex gain and ω is an unknown frequency error. u(m_n) is an unknown transmit data symbol, and element of the QAM constellation. Let y=x(M)+v where v is AWGN.

[0126] A first approach is to observe that |x(M)_n|²∈A where #A=k and absent noise |y_n|²=|x(M)_n|². (# denotes cardinality.) Let y² denote this vector. Thus (y²)_n∈A where #A=k.

k is the number of distinct amplitude rings. For QAM k>1 and typically small (i.e. 2 or 3). For FM k=1 and FM is distinguished from k-ary PSK because FM does not exhibit the k-ary PSK constellation.

[0127] It follows that one may use,

[0128] $\begin{array}{cc}{a\left(y,k\right)}_{\mathrm{QAM}k-\mathrm{ary}1}=\frac{\sum _i{\mathrm{Mean}\left(C_i\right)}^2}{\sum _i\mathrm{Variance}\left(C_i\right)}\mathrm{where}\left\{C_i\right\}=k-\mathrm{means}\left(y^2\right)& \left[\mathrm{Eq}.18\right]\end{array}$

[0129] As above, now consider what happens when samples are not necessarily taken at the center of transmit symbols. Begin by fixing an estimate of T and consider

[0130] $\tilde{T}\in \left[\frac{1}{2\tilde{\Delta }f}·\frac{2}{\tilde{\Delta }f}\right].$
Let S_i=∪_{t∈[i{tilde over (T)},(i+1){tilde over (T)}]}|y(t)|². For T={tilde over (T)} there exists s∈S_i such that s∈A.

[0131] Let {C_i}=k-means(∪_iS_i) and μ_i=Mean(C_i). The cluster {C_i} determine a partition of ∪_iS_i and μ_i are the cluster centers.

As before, let s′_i=arg min_u∈Si|μ_i−u|² and μ′_i=k-means(∪_is′_i), and iterate.

Thus,

[0132] $\begin{array}{cc}{a\left(y,k\right)}_{\mathrm{QAM}k-\mathrm{ary}2}={\max }_{\tilde{T}\in \left[\frac{1}{2\tilde{\Delta }f}·\frac{2}{\tilde{\Delta }f}\right]}\frac{\sum _i{\mathrm{Mean}\left(C_i\right)}^2}{\sum _i\mathrm{Variance}\left(C_i\right)}\mathrm{where}\left\{C_i\right\}=k-\mathrm{means}\left({\bigcup }_i{s}_i^{\prime }\right)& \left[\mathrm{Eq}.19\right]\end{array}$
Note that this algorithm is identical to the earlier algorithm for PSK. The difference is in the definition of S_i.

[0133] For PSK

[0134] $S_i={\bigcup }_{t\in \left[i\tilde{T}·\left(i+1\right)\tilde{T}\right]}y\left(t+\tilde{T}\right)\stackrel{\_}{y\left(t\right)}$
and for QAM S_i=∪_{t∈[i{tilde over (T)},(i+1){tilde over (T)}]}|y(t)|². There are also some minor differences in an appropriate implementation of k-means. It is still appropriate to retain a(y,k)_{QAM k-ary 1} as defined previously, for instances where the knowledge of symbol timing is present.
AM or DSB.

[0135] For x_n=s e^iωcnu_n where ℑu_n=0. That is, u_n is purely real. Note the use of ω_c rather than the earlier ω.

[0136] Recall, from the inverse Fourier Transform u_n=∫_ωs₋(ω)e^−ωn+s₊(ω)e^+ωn. ℑu_n=0 implies s₋(ω)e^−ωn and s₃₀(ω)e^+ωn are conjugates, and this implies that s₋(ω) and s₊(ω) are conjugates. Thus, |ℑ(u,ω)| is an even function of ω and Arg ℑ(u,ω) is an odd function of ω. The important observation is ℑ(u,ω)=ℑ(u,−ω). This is the characteristic that will be used to identify AM (or DSB).

[0137] Suppose one knows ω_c. Letting V(w_c,y,ω)=ℑ(e^−ωcny_n,ω) one may test for whether or not V(ω_c,y,ω)≈V(ω_c,y,−ω). This is the same relationship that exists between the commutated and reference samples in the presence of a target signal (in AWGN), in the construction of the rectangular decomposition. Hence, the following estimate is appropriate:

[0138] $\begin{array}{cc}{a\left(y\right)}_{\mathrm{AM}}=\frac{{<V\left({\tilde{\omega }}_c,y,\omega \right),V\left({\tilde{\omega }}_c,y,-\omega \right)>}^2}{\begin{array}{c}{V\left({\tilde{\omega }}_c,y,\omega \right)}^2{V\left({\tilde{\omega }}_c,y,-\omega \right)}^2-\\ {<V\left({\tilde{\omega }}_c,y,\omega \right),V\left({\tilde{\omega }}_c,y,-\omega \right)>}^2\end{array}}{\tilde{\omega }}_c=\arg {\max }_{{\omega }_c}{<V\left({\omega }_c,y,\omega \right),V\left({\omega }_c,y,-\omega \right)>}^2& \left[\mathrm{Eq}.20\right]\end{array}$
The first term in the denominator has no dependence on ω_c (by Parseval's theorem). Hence, a(y)_AM is maximized over ω_c when the numerator is maximal. This justifies the estimate {tilde over (ω)}_c.

[0139] Practically, the computation of {tilde over (ω)}_c is much simpler than it appears. One need only compute three FFT's. Compute one FFT transform, flip this over the ω-axis (a reversal), and correlate this with the unmodified transform to determine {tilde over (ω)}_c. This correlation may again be computed with two FFT's.

Implementation of Operation

[0140] Referring to FIG. 11, operation of the method of the present invention is implemented by a scan policy controller 60 in combination with a machine interface (MI) 62, a computer 64, signal observation resources 66 and a SOI database 68.

[0141] The MI 62 is embodied in the computer 64 (14 in FIG. 1) and the computer monitor (16 in FIG. 1) or may be embodied separately therefrom.

[0142] Also referring to FIGS. 1 and 2, the computer 64 (14 in FIG. 1) is programmed for constructing the statistical models of particular signals of interest 24 and the statistical model of the environment 25, and for processing the observed received signal during step 26 in accordance with the decision rule to decide whether the observed received signal is a SOI.

[0143] The signal observation resources 56 include antennas and receivers 10, 11, 12 that can be controlled to access specified attributes of the observed received signal for use in constructing the statistical models 24, 25.

[0144] The scan policy controller 60 is a collection of apparatus that control the signal observation resources 66 and the operation of the computer 64. The scan policy controller 60 is embodied in the computer 64 (14 in FIG. 1) and the computer monitor (16 in FIG. 1) or may be embodied separately therefrom.

[0145] The SOI database 66, which may be included in the computer 64, stores a statistical model for the environment 25 and statistical model for particular SOI's. Statistical models of target signals that are identified as signals of interest are entered into the SOI database upon affirmation by the operator that a particular target signal that has been identified as a SOI is indeed a SOI. The statistical models for particular SOI's are specified by the system operator.

[0146] When the operator desires a decision by the computer 64 as to whether an observed received signal is a particular SOI, the operator causes the statistical model for the particular SOI to be accessed by the scan policy controller 60 from the SOI database 66 and used for construction of the statistical model for the observed received signal; and the scan policy controller 60 decides whether the observed received signal is the particular SOI.

[0147] The MI 62 facilitates operation of the scan policy controller by a system operator by enabling the operator to specify and/or revise the attributes of a target signal that are used to construct the statistical model and/or to revise parameters of the decision rule in response to decisions made by the computer 64 (14 in FIG. 1) as to whether the observed received signal is a SOI.

[0148] A significant benefit of the methods described herein is reducing operator workload. This is realized by automating operator tasks that have historically been manual, wherever possible.

[0149] The benefits specifically stated herein do not necessarily apply to every conceivable embodiment of the present invention. Further, such stated benefits of the present invention are only examples and should not be construed as the only benefits of the present invention.

[0150] While the above description contains many specificities, these specificities are not to be construed as limitations on the scope of the present invention, but rather as examples of the preferred embodiments described herein. Other variations are possible and the scope of the present invention should be determined not by the embodiments described herein but rather by the claims and their legal equivalents.

[0151] Regarding the method claims, except for those steps that can only occur in the sequence in which they are recited, and except for those steps for which the occurrence of a given sequence is specifically recited or must be inferred, the steps of the method claims do not have to occur in the sequence in which they are recited.

• [0152] AM Amplitude modulation
• [0153] AOA Angle-of-arrival
• [0154] AR Autoregressive
• [0155] AWGN Additive white Gaussian noise
• [0156] CDF Cumulative distribution function
• [0157] CW Continuous wave
• [0158] DSB Double sideband
• [0159] FIR Finite impulse response
• [0160] FM Frequency modulation
• [0161] FSK Frequency shift keying
• [0162] GUI Graphical user interface
• [0163] KL Kullback-Leibler divergence
• [0164] LLR Log likelihood ratio
• [0165] LS Least-squares
• [0166] MLE Maximum-likelihood estimate (or estimator)
• [0167] MI Machine interface
• [0168] PDF Probability density function
• [0169] QAM Quadrature amplitude modulation
• [0170] SOI Signal-of-interest

[0171] $\frac{S}{N}$
Signal-to-noise ratio