Modulation Recognition Using Cyclic Cumulants, Part II: Algorithms

An example of how a decision-theoretic approach to a mathematical problem leads to a feature-based algorithmic solution that, viewed in isolation, appears to have sprung from the brow of Zeus fully formed.

At long last, I am picking up the modulation-recognition (MR) topic again. My plan is to describe the core MR algorithm I’ve developed using cyclic cumulants in this post, and then in a third and final post present some performance results using simulated and captured signals.

Recap of Part I Post

The Part-I post introduced the basic MR problem statement, provided four distinct important problem variations, and discussed the problem aspects that render the problem a perennial one.

The four variations are summarized next.

Radio-Frequency Scene Analysis. The most general version of the problem, as described in my introductory post and my wicked-problem post. For an arbitrary RF scene, detect, characterize, and classify all signals in the scene regardless of their durations, bandwidths, signal-to-interference-and-noise ratios (SINRs), types, and parameters. The cocktail-party problem for radios ™.

Conventional Modulation Classification. The least general version of the problem involves determining the type of a single signal in noise, with moderate-to-high SNR, that is at complex baseband or very close to it. This is the problem considered by most modulation-classification papers of the past four decades, including some of mine.

Cognitive-Radio Modulation Classification. In a modern twist on the problem, cognitive-radio classification involves quickly determining the type of a signal in each of a number of assigned frequency channels used by two kinds of users: Primary and Secondary. The Primary users are the licensed users, and they have right-of-way access to the channels. The Secondary users are non-licensed opportunistic users, attempting to make use of the channels while the Primary users are not. Secondaries must be able to determine what kind of signal is occupying a channel, even if the perceived inband SNR of that occupant is low.

Cochannel Modulation Classification. A difficult version of the problem occurs when two or more signals in the scene can be both contemporaneous and spectrally overlapping. This is called the cochannel signal condition. When the two signals have comparable power, then they each experience a signal-to-interference-and-noise ratio (SINR) of less than zero decibels. So typically the cochannel condition implies low SNR. Moreover, for feature-based classifiers, the features must be signal-selective to be effective, and for trained classifiers, such as trained convolutional neural networks, the training set needs to include all expected (or possible, depending on your biases) combinations of the signals, leading to a combinatorial explosion in the size of the training and testing datasets.

In the Part-I post I assert that the MR problem has not been definitively solved because of the large number of problem dimensions, including catalog size, SNR, cochannel complications, data-record length, sampling bandwidth, desire for real-time operation, and the amount and type of available prior information. I elaborated on this in my speculation that the MR problem is a wicked problem. I also discuss commonly used MR performance measures such as the confusion matrix for decision labels and square-root mean-squared error (RMSE) for parameter estimates. I then describe the major classes of solution approaches, including decision-theoretic, feature-based, and machine-learning. Finally, I present simulated examples that show exactly why various features fail in some important real-world situations.

In this post, I’ll describe my mathematical approach to modulation recognition (decision-theoretic leading to a final principled feature-based algorithm). In Part III, I’ll illustrate the approach with numerical examples.

Algorithm Approach

The mathematical approach to the MR problem is sketched in the companion papers My Papers [26,28], and has roots in a less well-motivated approach I took in My Papers [25]. The remainder of this section, and the section on algorithm development, is an explication of [25,26,28] without the heavy constraint of publication page limits.

Key Background Ideas

Hypothesis Testing

A core idea is that of hypothesis testing. When confronted with a data record containing a single signal with unknown type and parameters, we can imagine forming a sequence of hypotheses (proposed explanation) about that signal. On the first hypothesis, the signal is a BPSK signal. On the second, it is a QPSK signal. On the third, it is an OFDM signal. You can imagine going on like this for some time, leading to a set of, say, M hypotheses. Which one, if any, is true? You might add a final hypothesis such as None of the first M hypotheses is true. We would call that the None-of-the-Above (NOTA) or Unknown hypothesis.

The formation of hypotheses can be quite simple, as in the binary hypothesis testing problem (BHTP). In the BHTP, there are only two hypotheses. A typical example is hypothesis zero (called H_0) for which the signal is not present in the data and hypothesis one (H_1) for which the signal is present in the data. The problem is to decide between them in some principled or even optimal manner. I used this most basic problem set up in the cycle-detectors post.

Many MR situations are not simple, and this leads to the consideration of more complex sets of hypotheses. For example, the first hypothesis could be Type-A Signal is Present and the second could be Type-A Signal is Present with Cochannel Type-B Signal. If there are, say, M signal types of interest, and all two-signal cochannel combinations are possible (of interest), then you could formulate \binom{M}{2} distinct two-signal combinations. We could also include all M of the single-signal hypotheses and the NOTA hypothesis. That would lead to a multi-hypothesis problem having

\displaystyle \binom{M}{2} + M + 1 = \frac{M!}{(M-2)!2!} + M + 1 \hfill (1)

hypotheses. If you have M=20 signal types of interest, this leads to 211 hypotheses. (Ouch.)

Complete Probabilistic Characterizations

We’ll have to carefully consider which hypotheses we include in an MR problem statement, but the discussion above avoided concrete examples of the contents of an individual hypothesis. We said things similar to ‘on hypothesis seven, the signal is BPSK.’ But what exactly does ‘signal is BPSK’ mean here, mathematically? This goes to the heart of my musings on BPSKness in several CSP Blog posts, such as All BPSK Signals and Is RFSA a Wicked Problem?

When we say that on the seventh hypothesis, the data is composed of BPSK plus noise, what does that mean? Here are some possibilities:

  1. The inphase component of the data (real part) has a graph (plot of the amplitude versus time) that is consistent with that of a BPSK signal with particular values of the pulse-shaping function, bit rate, carrier frequency offset, and allowable transmitted bit sequences.
  2. The estimated power spectral density is consistent with that of the power spectral density of a BPSK signal with particular parameter values. If you subtract the power-spectrum estimate from the theoretical power spectrum, and summed the absolute value of the resulting difference, the number would be small compared to the energy in the power spectra.
  3. The set of detectable second-order cycle frequencies for the signal in the data conforms to the pattern of second-order cycle frequencies for a BPSK signal with modern modulation-efficient pulses, such as square-root raised-cosine pulses. The non-conjugate pattern is a single non-zero cycle frequency \alpha_1. The conjugate pattern is a triplet of cycle frequencies \{\beta_1, \beta_2, \beta_3\} with \beta_3 > \beta_2 > \beta_1 and \beta_3 - \beta_2 = \alpha_1 = \beta_2 - \beta_1.

But, what is the optimal way to form the hypothesis? What is the complete description of BPSKness? If we can capture the full description of what it means for a data sequence or a continuous voltage record to correspond to BPSK–and only to BPSK–then that description would be a good one to use in hypothesis formulation.

Are any of the three alternatives above complete descriptions, in some relevant sense, of BPSKness? Let’s discuss (remember, no page limits here!).

For the first description alternative, we can conceive of a table containing all possible inphase sampled-data or continuous-time functions that correspond to an instance of the BPSK signal in question. That is, fix the bit rate, carrier frequency offset, carrier phase, symbol clock phase, pulse-shaping function, and power (amplitude), and create a table entry for each possible sequence of transmitted bits. We could then do that for all combinations of bit rate, carrier frequency offset, carrier phase, symbol-clock phase, pulse-shaping function, and power. This leads to a very large table–infinite if we take seriously all possible values of, say, the carrier phase (which is an uncountably infinite set because it is a real number).

In the end, we have a table with infinite rows each containing a finite-length complex-valued sequence or function of time. Here, we can formulate a hypothesis for BPSK such as On Hypothesis Seven, the Signal in the Given Data Appears in the BPSK Table. How to make use of this kind of hypothesis, provided we actually have the table? We could form the difference between the data and each row of the table and retain the index of the row for which the sum of the absolute error is minimum. Then repeat that for similar tables for other modulations of interest (you can see the infinities stacking up here I hope). Finally, choose as the modulation-type decision the one corresponding to the minimum table error. For this to be good, we’d have to ensure that the tables for different signal types don’t have significant numbers of identical entries. (Ponder the question of whether a table for QPSK and a table for BPSK will have any identical entries.)

Alternative one isn’t practical–we can’t make these ideal tables and trying to make approximations will also be difficult due to combinatorial explosion. This idea of BPSKness as a giant list of voltage trajectories does bring to mind modern approaches to MR that attempt to create very large I/Q training datasets, which are nothing more than truncated versions of the tables I’ve described. But we have seen over the past five years where that gets us.

The second description alternative is to describe the essential nature of the signal in terms of its power spectrum. Does the estimated power spectrum for the input data match, in some suitable numerical sense, a stored ideal (theoretical) power spectrum for one of the signals of interest? A problem with this alternative is that the power spectrum is not unique for many sets of signal types of interest. I’ve belabored this in other contexts, and it is easily seen to be true. For example, all unit-power PSK and QAM signals with balanced constellations, identical pulse-shaping functions, and independent and identically distributed symbol sequences have identical power spectra.

The third alternative above involves the pattern of second-order cycle frequencies for a signal. In this description, BPSKness involves two small sets of real numbers together with their differences and the relation between the elements of the two sets. If we consider three different signal types: BPSK, QPSK, and AM-DSB, and denote the symbol rate of the PSK signals as 1/T_0 and the carrier frequency of all three signals as f_c, then the signals are described by these cycle-frequency patterns:

BPSK: \{\{1/T_0\}, \{2f_c - 1/T_0, 2f_c, 2f_c+1/T_0\}\}

QPSK: \{\{1/T_0\}, \emptyset\}

AM-DSB: \{\emptyset, \{2f_c\}\}

The three sets (second-order cycle-frequency patterns) are highly distinct, which is good for MR. But the flaw in this signalness description is similar to the flaw for the second alternative description that used the power spectrum: A large number of distinct signal types of interest also possess the same cycle-frequency pattern as QPSK. Examples include 8QAM, 16QAM, 64QAM, \pi/4-DQPSK, a huge number of continuous-phase-modulated (CPM) signals (My Papers [8]), and others. In fact, there are several distinct signal types that map to the BPSK second-order cycle-frequency pattern: All pulse-amplitude-modulated (PAM) signals with symmetric real-valued constellations. These are sometimes called amplitude-shift-keyed (ASK) signals, such as 4ASK, 8ASK, 16ASK, etc., for which the symbol constellation points lie on the real axis in the complex plane.

So a big obstacle in the way of setting up a multiple-hypothesis MR problem is how to describe the hypotheses such that each hypothesis captures something unique about the associated signal. And at the same time, avoiding descriptions that are infinite, like the tables above. Some plausible descriptions end up forming undesirable equivalence classes–multiple distinct signals will have identical descriptions in their hypotheses.

But we know that, for instance, QPSK and 16QAM are distinct signal types even if their PSD-based descriptions form an equivalence class. Or do we? Why do we think they are distinct? The question becomes: What does a signal label mean? If all the 16QAM signals are crammed into a box labeled 16QAM, how do we construct the contents of that box? Exactly how do we do that?

I draw your attention back to alternative description one, the infinite table of I/Q trajectories. We’ve seen that before! In the context of random processes, it is the ensemble. And we have a way of generating the entries in the table–the elements of the ensemble–using mathematics. We can write an expression that covers all of the uncountable infinity of ensemble elements; that expression is the ensemble generator. For 16QAM, it is

\displaystyle X(t) = A \sum_{k=-\infty}^\infty a_k p(t-kT_0-t_0)e^{-i2\pi f_0 t + i \phi_0}, \hfill (2)

where A, a_k, T_0, t_o, f_0, and \phi_0 can be random variables. These are the amplitude (determines power), transmitted symbols, symbol interval (reciprocal of symbol rate), symbol-clock phase, carrier (offset) frequency, and carrier phase, respectively.

Most generally, a random process is completely specified by the collection of all possible nth-order joint probability density functions for sets \{X(t+\tau_1), X(t+\tau_2), \ldots, X(t+\tau_n)\}. For the random process in (2), these joint density functions can be determined from the probability density functions for each of the involved random variables.

If we wish, we can add more random variables to this process. For example, we can introduce a propagation-channel model such as the discrete multipath channel. For every ensemble element generated by (2), we’d have another uncountable infinity of new ensemble elements, each produced by applying a randomly chosen channel to the chosen element.

So we do have one answer to the question of what constitutes BPSKness: An infinite set of multidimensional (joint) probability density functions. Once those are specified, we can create a random process that will produce all BPSK signals conforming to the underlying random-variable distributions, which is an uncountably infinite set.

The final step in creating hypothesis descriptions is to realize that the set of all possible joint probability density functions for the random process can be expressed in terms of all possible nth-order temporal moment functions–we can sidestep writing down the density functions and just write down the moments. And we already know that we can write down all the joint moments if we know all the joint cumulants. Therefore, our description of signalness can take the following equivalent forms:

  1. The ensemble: A table with an uncountable infinity of rows, each corresponding to one instance of the signal, which means a sample path of the signal’s random-process characterization.
  2. All joint probability density functions for sets like \{X(t+\tau_1), X(t+\tau_2), \ldots, X(t+\tau_n)\}.
  3. All nth-order temporal moment functions for products like \hspace{1.0in} \prod_{j=1}^n X^{(*)_j}(t+\tau_j).
  4. All nth-order temporal cumulant functions for sets like \{X(t+\tau_1), X(t+\tau_2), \ldots, X(t+\tau_n)\}.
Likelihood

Consider a random variable X and an associated hypothesis-testing problem with hypotheses H_k. On hypothesis zero, H_0, the random variable has density function f_X(x|H_0) and on hypothesis one, it has density f_X(x|H_1). The likelihood ratio is given by

\displaystyle \gamma_X (x) = \frac{f_X(x|H_1)}{f_X(x|H_0)}. \hfill (3)

When this ratio is greater than one, the interpretation is that the observed data (x) is more likely to have arisen under hypothesis H_1 than under H_0. If the ratio takes an exponential form, the natural logarithm is often applied to gain access to the exponent, which is reasonable because the logarithm is monotonic, so no decisions will be inappropriately reversed. The log-likelihood ratio (LLR) is then

\displaystyle \mathcal{L}_X(x) = \ln \left[ \gamma_X(x) \right] \hfill (4).

Suppose the random variable X is Gaussian on H_i with mean m_i and variance \sigma_i^2. Then the two densities are easily expressed as

\displaystyle f_X(x|H_i) = \frac{1}{\sqrt{2\pi\sigma_i^2}} e^{-(x-m_i)^2/(2\sigma_i^2)} \hfill (5)

A maximum-likelihood decision would then compute each of the likelihoods and choose the hypothesis corresponding to the maximum of those likelihoods.

Problem Statement Evolution

The goal is now to craft a problem statement that captures the modulation-recognition problem and that is mathematically tractable. We also want the problem statement to be as general as possible, so that it is applicable to a large set of modulated signals. To achieve this goal, we present a sequence of problem statements. At first, the problem statements will be equivalent, but eventually we will introduce variations or approximations that result in problem statements that are tractable, and that are not equivalent to the most general problems from the beginning of the sequence.

We want to process some sampled data and decide if that data contains one of L_B distinct signal types drawn from some catalog of types B. That is, each element of the set B is a label z and there are L_B such labels. Leaving aside the NOTA or Unknown decision, we have L_B hypotheses.

On the ith hypothesis H_i, the data contains a signal with type z_i.

For each signal type z_i in B, we conceive of the collection of all possible z_i-type signals as a random process characterizing signal z_i, which means we know all the joint PDFs for that random process. For convenience, we assume that the autocorrelation function for the z_i random process has a value of one at lag \tau = 0–the elements of B are unit-power random processes. (We do not make this assumption about the power of the received signal under study.)

Now let’s consider some label z \in B (drop the subscript i for convenience). Denote by \mathcal{P}_z the collection of all nth-order probability density functions for signal type z. The (infinite) set \mathcal{P}_z contains the complete probabilistic description of random process z. For any signal (sample path) y(t) we use the notation

\displaystyle y(t) \sim \mathcal{P}_z \hfill (6)

which we read as the signal y(t) is distributed as \mathcal{P}_z.

Problem One

Let’s consider the multiple hypothesis problem given by

H_i: \ \ \ x(t) = a s_i(t-t_0) + w(t), \hfill s_i(t) \sim \mathcal{P}_i \hfill (7)

for i = 1, \ldots, L_B. Choosing hypothesis j is equivalent to declaring that x(t) contains the signal with label j plus white Gaussian noise (w(t)).

Although it is possible to know \mathcal{P}_i for many signal types, it is typically difficult to obtain that information, and it is even more difficult to measure the elements of \mathcal{P}_i from data. We are thus inspired to reformulate the problem.

For an arbitrary probability density function, say

\displaystyle P_s(t, \boldsymbol{\Delta}; n, \boldsymbol{\tau}) = \mbox{\rm PDF} \left\{ s(t+\tau_j) \right\}_{j=1}^n, \hfill (8)

that depends on some vector of random variables \boldsymbol{\Delta}, we can find the characteristic function by Fourier transformation in the \boldsymbol{\Delta} variables,

\displaystyle \Phi_s(t, \boldsymbol{\omega}; n, \boldsymbol{\tau}) = \int_{-\infty}^\infty \cdots \int_{-\infty}^\infty P_s(t, \boldsymbol{\Delta}; n, \boldsymbol{\tau}) e^{i 2 \pi \boldsymbol{\omega}\boldsymbol{\Delta}^\dagger} \, d\boldsymbol{\Delta} .\hfill (9)

Moreover, the characteristic function itself is a Taylor series involving the moment functions for the sets \displaystyle \{s(t+\tau_j)\}_{j=1}^n, or

\displaystyle \left. \frac{\partial} {\partial\boldsymbol{\omega}} \Phi_s(t, \boldsymbol{\omega}; n, \boldsymbol{\tau})\right|_{\boldsymbol{\omega} = \boldsymbol{0}} = (-i)^{-n} M_s(t, \boldsymbol{\tau};n), \hfill (10)

where the nth-order time-varying moment is given by

\displaystyle M_s(t, \boldsymbol{\tau};n) = E\left[ \prod_{j=1}^n s(t+\tau_j) \right]. \hfill (11)

All the possible moments are found in the multidimensional Taylor-series expansion of the characteristic function. For example, the pth-order moment corresponding to the product of possibly repeated terms

\displaystyle \prod_{j=1}^r s(t+\tau_j)^{q_j}

where

\displaystyle \sum_{j=1}^r q_j = p, \ \ \ q_j \ge 1, \hfill (12)

is given by the coefficient in the series that corresponds to the term

\displaystyle \frac{\displaystyle i^p \prod_{j=1}^r \omega_j^{q_j}}{\displaystyle \prod_{j=1}^r q_j!}

The point here is that the characteristic functions are completely specified by the set of all possible nth-order moments for s(t). In turn, due to the uniqueness of the Fourier transform (under mild conditions on the involved functions), the characteristic functions completely specify the set of all possible joint probability density functions.

Let’s use some compact notation to try to capture these ideas, which will be critical for our succession of hypotheses-laden problem formulations. Let the symbol \mathcal{M}_i denote all possible moment functions for signals with the ith label in B. Then the following two statements are equivalent:

\displaystyle s_i(t) \sim \mathcal{P}_i \hfill (13)

\displaystyle s_i(t) \sim \mathcal{M}_i. \hfill (14)

Saying that a signal is characterized by an infinite set of probability density functions is equivalent to saying that the signal is characterized by an infinite set of nth-order moments, provided the moments correspond to the density functions.

Finally, we make the connection to cumulants. As we saw in the SPTK post on the characteristic function, the natural logarithm of the characteristic function has a series representation in which the coefficients of the various orders of the frequency variables are the cumulants. That is, in the notation we’re using here,

\displaystyle \left. \frac{\partial}{\partial \boldsymbol{\omega}} \ln \Phi_s (t, \boldsymbol{\omega}; n, \boldsymbol{\tau}) \right|_{\boldsymbol{\omega}=\boldsymbol{0}} = (-i)^{-n} C_s (t, \boldsymbol{\tau}; n). \hfill (15)

So if you knew all the cumulants, you could construct the logarithm of the characteristic function, and if you know the logarithm of the characteristic function, you can use exponentiation to obtain the characteristic function itself, and if you know the characteristic function, you know the probability density functions. Let’s denote the infinite set of cumulants by \displaystyle \mathcal{C}. Then we arrive at the equivalencies

\displaystyle s_i(t) \sim \mathcal{P}_i \Longleftrightarrow s_i(t) \sim \mathcal{M}_i \Longleftrightarrow s_i(t) \sim \mathcal{C}_i. \hfill (16)

I emphasize “equivalencies” here–this is not approximately true, it is exactly true. We can now formulate a second problem.

Problem Two

The original problem (7), posed in terms of the infinite set of joint probability density functions \mathcal{P}, is equivalent to a problem posed in terms of an infinite set of moment functions \mathcal{M},

H_i: \ \ \ x(t) = as_i(t-t_0) + w(t), \hfill s_i(t) \sim \mathcal{M}_i \hfill i = 1, \ldots, L_B. \hfill  (17)

On the ith hypothesis, the data x(t) contains a noisy signal that has moments \mathcal{M}_i. We have good reasons for preferring cumulants over moments, which leads to Problem Three.

Problem Three

The problems (7) and (17) are equivalent to a problem posed in terms of an infinite set of cumulant functions \mathcal{C},

H_i: \ \ \ x(t) = as_i(t-t_0) + w(t), \hfill s_i(t) \sim \mathcal{M}_i \hfill i = 1, \ldots, L_B. \hfill  (18)

We can also pose problems in terms of the distribution, moments, or cumulants of the data itself, rather than just for the signal component of the data. Let’s use the notation \mathcal{P}_{i+w} to refer to all the probability density functions for the sum s_i(t) + w(t). To fully capture the probabilistic nature of the data, let the symbol \mathcal{P}_{i+w}(a, t_0) denote the collection of all probability density functions for a s_i(t-t_0) + w(t). The symbols \mathcal{M}_{i+w}(a, t_0) and \mathcal{C}_{i+w}(a, t_0) similarly denote the sets of moments and cumulants, respectively.

Problem Four

H_i: \ \ \ x(t)  \sim \mathcal{P}_{i+w}(a, t_0) \hfill i = 1, \ldots, L_B. \hfill  (19)

Problem Five

H_i: \ \ \ x(t)  \sim \mathcal{M}_{i+w}(a, t_0) \hfill i = 1, \ldots, L_B. \hfill  (20)

Problem Six

H_i: \ \ \ x(t)  \sim \mathcal{C}_{i+w}(a, t_0) \hfill i = 1, \ldots, L_B. \hfill  (21)

Problems Four, Five, and Six are actually more general than the previous problems because they also include the full probability structure of the noise w(t).

We would now like to develop a problem statement that deals with measurable quantities and one with a solution that can be implemented (in software) for a catalog B of interest. It turns out that both probability density functions and higher-order moments are difficult to estimate or use for various reasons in many scenarios, but that cumulants are relatively easy to estimate and lead to useful algorithm generalizations in scenarios involving multiple cochannel signals and noise with unknown parameters.

Elements of \mathcal{C}_i take the form

\displaystyle C_{s_i}(t, \boldsymbol{\tau};n) = \sum_{P_n} k(p) \prod_{j=1}^n R_{s_i}(t, \boldsymbol{\tau}_{\nu_j}; n_j) \hfill (22)

where P_n is the set of all distinct partitions of the index set \{1, 2, \ldots, n\}. Each element of P_n is a collection of p subsets of the index set such as \{ \nu_1, \nu_2, \ldots, \nu_p\} where the union of the subsets \nu_k is the index set and the p subsets are all distinct (their intersections are the emptyset). The size of \nu_k is n_k, so that the sum of the n_k is n.

In (22), the symbol R_s(t, \boldsymbol{\tau}_{\nu_j}; n_j) denotes the n_jth-order temporal moment function for s(t). It will be helpful in the following progression of problem statements to recall that cumulants for Gaussian random variables, processes, and signals are zero for n \ge 3.

Problem Seven

H_i: C_x(t, \boldsymbol{\tau};n) = a^n C_{s_i}(t-t_0, \boldsymbol{\tau}; n) + C_w(t, \boldsymbol{\tau};n), \  i = 1, \ldots, L_B,\   n=1, 2, \ldots  \hfill  (23)

Here we are saying that on the ith hypothesis, the observed data has a set of cumulants that are equal to the sum of the cumulants for the ith signal, s_i(t), plus the cumulants for the noise component of the data.

As an aside, we can imagine at this point extending Problem Seven to include cochannel signals. We could have a matrix of hypotheses H_{i,j}, for example, where the (i,j)th hypothesis means that the data has a set of cumulants equal to the sum of those for signal type i, signal type j, and the noise component,

H_{i,j}: C_x(t, \boldsymbol{\tau};n) = a_i^n C_{s_i}(t-t_1, \boldsymbol{\tau}; n) + a_j^n C_{s_j}(t-t_2, \boldsymbol{\tau}; n) + C_w(t, \boldsymbol{\tau};n),   n=1, 2, \ldots  \hfill  (24)

Problem Seven gets us closer to a set of hypotheses that involve measured quantities, but isn’t quite there yet, since the specified cumulants are the theoretical cumulants, and so require an infinite ensemble or an infinite-duration signal to average over.

The cumulants can be estimated directly from the data given knowledge of their variation with t. Moreover, the standard estimates are asymptotically independent and Gaussian.

There is a serious dimensionality problem remaining regarding the lag vector \boldsymbol{\tau} and the order variable n. To move toward mitigating the high dimensionality of the quantities involved in the hypotheses, we now make the move from time-varying cumulants to their Fourier-series components, the cyclic cumulants.

For the RF communication signals of typical interest to RF scene analysis and modulation recognition, the cumulants are periodically or polyperiodically time-variant, and so we have the Fourier-series representation encountered in the cyclic cumulant post,

\displaystyle C_{s_j}(t, \boldsymbol{\tau}; n) = \sum_{\beta_j (n)} C_{s_j}^{\beta_j(n)} (\boldsymbol{\tau};n) e^{i 2\pi \beta_j(n) t} \hfill (25)

where \displaystyle \beta_j(n) are the cycle frequencies for the jth signal type for order n.

We note the equivalence of the time-varying cumulant for signal type j and the set of cyclic cumulants for signal type j,

\displaystyle C_{s_j}(t, \boldsymbol{\tau}; n) \Longleftrightarrow \left\{ \left\{ C_{s_j}^{\beta_j(n)} (\boldsymbol{\tau}; n) \right\}, \left\{ \beta_j(n) \right\} \right\}. \hfill (26)

Problem Eight

In this next problem formulation, the hypotheses relate to sets of cyclic cumulants rather than to sets of time-varying cumulants, which eliminates the time variable t from our setup,

\displaystyle H_i: C_x^{\beta_i(n)} (\boldsymbol{\tau}; n) = a^n C_{s_i}^{\beta_i(n)} (\boldsymbol{\tau}; n) e^{-i2 \pi \beta_i (n)t_0} + C_w^{\beta_i(n)}(\boldsymbol{\tau};n) \hfill (27)

for i= 1, \ldots, L_B, n=1, 2, \ldots, all cycle frequencies \beta_i(n) for which the signal type has a non-zero cyclic cumulant, and \boldsymbol{\tau} \in \mathbb{R}^n.

So far we have not sacrificed any optimality since the information exploited in Problem Eight is identical to that exploited in Problem One, and in both cases the information consists of the complete (exhaustive) probabilistic description of the involved signal types.

However, the dimensionality of the involved quantities is very large for n > 2. The high-dimensional variables include the n-dimensional delay vector \boldsymbol{\tau}, all orders n, and the potentially infinite sets of cycle frequencies \beta_i(n). While most communication signals possess a finite set of cycle frequencies for any finite order n, some, such as the rectangular-pulse PSK signals. possess infinite sets of cycle frequencies.

We have two basic choices for dimensionality reduction: restrict the order n or restrict the lag vector \boldsymbol{\tau}. Before delving into the dimensionality-reduction issue, however, let’s make one final modification to our problem statement.

Problem Nine

The modification consists of deleting (neglecting) the cyclic cumulant that corresponds to signal power, which means the cyclic cumulant for n=2 and \beta_i(2) =0. This is also known as the autocorrelation function.

\displaystyle H_i: C_x^{\beta_i(n)}(\boldsymbol{\tau}; n) = a^n C_{s_i}^{\beta_i(n)}(\boldsymbol{\tau}; n)e^{-i 2 \pi \beta_i(n)t_0}, \hfill (28)

for i= 1, \ldots, L_B, n=1, 2, \ldots, all cycle frequencies \beta_i(n) for which the signal type has a non-zero cyclic cumulant except \beta_i(2)=0, and \boldsymbol{\tau} \in \mathbb{R}^n.

We have left off any mention of the cyclic cumulants for the white Gaussian noise process w(t) in Problem Nine because of the peculiarities of the higher-order cumulants for Gaussian variables and signals,

\displaystyle C_w^{\beta_i(n)} (\boldsymbol{\tau};n) = \left\{ \begin{array}{ll} 0, & n \ge 3 \\ 0, & \beta_i(n) \neq 0 \\ 0, & \boldsymbol{\tau} \neq 0 \\ \sigma_n^2, & \mbox{\rm otherwise}. \end{array} \right. \hfill (29)

We thereby avoid the “second-order stationary” probabilistic parameter but also lose optimality, and we gain sweet simplicity. Moreover, a case of great practical interest corresponds to C_{s_i}^0(\boldsymbol{\tau};2) = f(\boldsymbol{\tau}) for multiple values of i. That is, the signals have identical power spectra and therefore the power spectrum and, equivalently, the autocorrelation, do not help to distinguish these signals from each other.

Another good reason for leaving the autocorrelation out of the hypotheses is that when we generalize the hypotheses to the case of cochannel signals, that parameter reflects all of the cochannel signals, not just the ith. This is also true of cyclic cumulants for (n,m,k) = (n, n/2, 0), which are the traditional stationary-signal higher-order statistics.

Problem Ten: The Order-Reduced Classifier

A major simplification of Problem Nine is to severely restrict the order n. In particular, we allow a single value of the order: n=2. This leads to the suboptimal (because we are throwing away almost all of the probabilistic information about the involved signals) hypotheses given by

\displaystyle H_i: C_x^{\beta_i(2)}(\boldsymbol{\tau}; 2) = a^2 C_{s_i}^{\beta_i(2)}(\boldsymbol{\tau};2) e^{-i2\pi \beta_i(2)t_0}, \hfill (30)

for i=1, 2, \ldots, L_B, \boldsymbol{\tau} \in \mathbb{R}^n, and \beta_i(2) \neq 0.

Problem Eleven

The set of hypotheses (30) suggests the following practical problem, where the limit probabilistic parameter for the received data is replaced by a noisy estimate,

\displaystyle H_i: \hat{C}_x^{\beta_i(2)}(\boldsymbol{\tau};2) = a^2 C_{s_i}^{\beta_i(2)} (\boldsymbol{\tau};2) e^{-i2\pi\beta_i(2)t_0} + v(\beta_i(2), \boldsymbol{\tau}),

for i=1, 2, \ldots, L_B, \boldsymbol{\tau} \in \mathbb{R}^n, \beta_i(2) \neq 0, and where v(\cdot) is white Gaussian noise.

Problem Eleven is recognized as a known-signal-in-noise problem (a matched-filtering problem). The problem is easy when the cycle frequencies for the different hypotheses are disjoint–that is, when the signals all have different cycle frequencies. We focus on the harder problem, which is of substantial practical interest, where all the signals have the same cycle frequency sets,

\displaystyle \left\{ \{\beta_i(2)\}\right\}_{i=1}^{L_B} = \{ \beta(2)\} for all indices i = 1, 2, \ldots, L_B.

Problem Eleven is then simplified to

\displaystyle H_i: \hat{C}_x^{\beta(2)}(\boldsymbol{\tau};2) = a^2 C_{s_i}^{\beta(2)} (\boldsymbol{\tau};2) e^{-i2\pi\beta(2)t_0} + v(\beta(2), \boldsymbol{\tau}), \hfill (31)

for i=1, 2, \ldots, L_B, \boldsymbol{\tau} \in \mathbb{R}^n, \beta(2) \neq 0.

Further simplifications are possible. Note that the information in a single cyclic cumulant function requires only n-1 lag variables \tau_j to capture. That is, we usually work with the reduced-dimension cyclic temporal moment and cumulant functions

\displaystyle C_s^\beta(\boldsymbol{\tau}-\boldsymbol{1}\tau_n; 2) = \bar{C}_s^\beta (\boldsymbol{u};n) e^{-i2 \pi \beta \tau_n} \hfill (32)

where the n-1-dimensional lag vector \boldsymbol{u} is given by

\displaystyle \boldsymbol{u} = \left[ \tau_1-\tau_n, \tau_2-\tau_n, \ldots, \tau_{n-1}-\tau_n\right]. \hfill (33)

More generally, we simply use the n-1-dimensional function \bar{C}_s^\beta(\boldsymbol{u};n) with arbitary \boldsymbol{u}.

Let’s name the second-order cycle frequencies \gamma_j as in \displaystyle \left\{ \gamma_j \right\}_{j=1}^{N_2} = \left\{ \beta(2) \right\}. Then Problem Eleven can be expressed as

\displaystyle H_i: \hat{\bar{C}}_x^{\gamma_j} (\tau; 2) = a^2 \bar{C}_{s_j}^{\gamma_j}(\tau; 2) e^{-i2\pi \gamma_j t_0} + v(\gamma_j, \tau), \hfill (34)

for i=1, 2, \ldots, L_B, \tau \in \mathbb{R}, \gamma_j \neq 0, and j = 1, 2, \ldots, N_2.

Finally, we make the simplifying assumption that the noise v(\gamma_j, \tau) is Gaussian, independent of j, and white so that we can simplify the notation to v(\gamma_j, \tau) = v(\tau).

Solution to Problem Eleven: The ORC

If the amplitude a and the delay (relative to the assumed signal model’s temporal origin) t_0 are known, then (34) is truly a known-signal problem. However, these parameters need to be estimated in practice because we can’t assume we know the amplitude and delay of a signal we are trying to recognize. We find these later; for now, assume we have available estimates \hat{a}_i and \hat{t}_{0,i}.

For each \gamma_j and \tau, \hat{\bar{C}}_x^{\gamma_j} (\tau; 2) is an independent Gaussian random variable with mean value

\displaystyle \hat{a}_i^2 \bar{C}_{s_i}^{\gamma_j} (\tau, 2) e^{-i2\pi \gamma_j \hat{t}_{0,i}}, \hfill (35)

and variance \sigma_v^2. Assume each cyclic cumulant has approximately finite support contained within the interval [-T, T] for T large enough. Form the vector \hat{\bar{\boldsymbol{C}}} by separating the real and imaginary parts of the cyclic cumulant,

\displaystyle \hat{\bar{\boldsymbol{C}}} =  \left[ \hat{\bar{C}}_{re} \ \ \hat{\bar{C}}_{im} \right]

where

\displaystyle \hat{\bar{C}}_{re} = \Re \left[ \underbrace{\hat{\bar{C}}_x^{\gamma_1} (-T, 2) \cdots \hat{\bar{C}}_x^{\gamma_1} (T, 2)}_{\gamma_1} \cdots    \underbrace{\hat{\bar{C}}_x^{\gamma_{N_2}} (-T, 2) \cdots \hat{\bar{C}}_x^{\gamma_{N_2}} (T, 2)}_{\gamma_{N_2}} \right] \hfill (36)

\displaystyle \hat{\bar{C}}_{im} = \Im \left[ \underbrace{\hat{\bar{C}}_x^{\gamma_1} (-T, 2) \cdots \hat{\bar{C}}_x^{\gamma_1} (T, 2)}_{\gamma_1} \cdots    \underbrace{\hat{\bar{C}}_x^{\gamma_{N_2}} (-T, 2) \cdots \hat{\bar{C}}_x^{\gamma_{N_2}} (T, 2)}_{\gamma_{N_2}}\right] \hfill (37)

Through this sequence of definitions, approximations, and assumptions (some quite probably questionable), we arrive at a vector \hat{\bar{\boldsymbol{C}}} where each component is a real random variable with known mean and variance \sigma_\nu^2/2. Let’s use the simplified notation

\displaystyle \hat{\bar{C}}_x^{\gamma_j} (k, 2)_r = \Re \left[ \hat{\bar{C}}_x^{\gamma_j} (k, 2) \right] \hfill (38)

\displaystyle \hat{\bar{C}}_x^{\gamma_j} (k, 2)_i = \Im \left[ \hat{\bar{C}}_x^{\gamma_j} (k, 2) \right] \hfill (39)

to enable us to write the joint probability density function for the elements of \hat{\bar{\boldsymbol{C}}} as the product of 2\times (2T+1) \times N_2 marginal Gaussian distribution functions,

\displaystyle P(\hat{\bar{\boldsymbol{C}}}|H_i) = \prod_{j=1}^{N_2} \left[ \prod_{k=-T}^T \left(\frac{1}{\sqrt{2}\pi (\sigma_\nu^2/2)} \right) e^{\left(-\frac{[\hat{\bar{C}}_x^{\gamma_j}(k,2)_r - E[\hat{\bar{C}}_x^{\gamma_j}(k,2)_r] ]^2}{\sigma_\nu^2}\right)}\right.

\displaystyle \times \left. \prod_{l=-T}^T \left(\frac{1}{\sqrt{2}\pi (\sigma_\nu^2/2)} \right) e^{\left(-\frac{[\hat{\bar{C}}_x^{\gamma_j}(l,2)_i - E[\hat{\bar{C}}_x^{\gamma_j}(l,2)_i] ]^2}{\sigma_\nu^2}\right)}\right] \hfill (40)

Let’s simplify the notation for the mean values,

\displaystyle D_{s_i}^{\gamma_j} (\tau, 2)_r = \Re \left[ \hat{a}_i^2 \bar{C}_{s_i}^{\gamma_j} (\tau, 2) e^{-i2\pi \gamma_j \hat{t}_{0,i}} \right] \hfill (41)

\displaystyle D_{s_i}^{\gamma_j} (\tau, 2)_i = \Im \left[ \hat{a}_i^2 \bar{C}_{s_i}^{\gamma_j} (\tau, 2) e^{-i2\pi \gamma_j \hat{t}_{0,i}} \right] \hfill (42)

\displaystyle P(\hat{\bar{\boldsymbol{C}}}|H_i) = \prod_{j=1}^{N_2} \left[ \prod_{k=-T}^T \left(\frac{1}{\sqrt{2}\pi (\sigma_\nu^2/2)} \right) e^{\left(-\frac{[\hat{\bar{C}}_x^{\gamma_j}(k,2)_r - D_{s_i}^{\gamma_j}(k,2)_r ]^2}{\sigma_\nu^2}\right)}\right.

\displaystyle \times \left. \prod_{l=-T}^T \left(\frac{1}{\sqrt{2}\pi (\sigma_\nu^2/2)} \right) e^{\left(-\frac{[\hat{\bar{C}}_x^{\gamma_j}(l,2)_i - D_{s_i}^{\gamma_j}(l,2)_i ]^2}{\sigma_\nu^2}\right)}\right] \hfill (43)

Our density function can be slightly simplified to

\displaystyle P(\hat{\bar{\boldsymbol{C}}}|H_i) =\left(\frac{\sqrt{2}}{\pi (\sigma_\nu^2/2)} \right)^{2N_2(2T_1)} \prod_{j=1}^{N_2} \left[ \prod_{k=-T}^T  e^{\left(-\frac{[\hat{\bar{C}}_x^{\gamma_j}(k,2)_r - D_{s_i}^{\gamma_j}(k,2)_r ]^2}{\sigma_\nu^2}\right)}\right.

\displaystyle \times \left.   e^{\left(-\frac{[\hat{\bar{C}}_x^{\gamma_j}(k,2)_i - D_{s_i}^{\gamma_j}(k,2)_i ]^2}{\sigma_\nu^2}\right)}\right] \hfill (44)

We can now compute the logarithm of the density function (which is also a likelihood function), which will simplify finding its extremal values, leading to a method of signal classification. That is, we want to find the index i for which the density function is maximum.

The logarithm converts the products to sums,

\displaystyle \ln P(\hat{\bar{\boldsymbol{C}}}| H_i) = K_1 + \frac{1}{\sigma_\nu^2} \sum_{j=1}^{N_2} \sum_{k=-T}^T \left[ -\hat{\bar{C}}_x^{\gamma_j}(k,2)_r^2 - \hat{\bar{C}}_x^{\gamma_j}(k, 2)_i^2 -D_{s_i}^{\gamma_j} (k,2)_r^2 - D_{s_i}^{\gamma_j}(k,2)_i^2 \right.

\displaystyle \left. +2\hat{\bar{C}}_x^{\gamma_j}(k,2)_rD_{s_i}^{\gamma_j}(k,2)_r + 2\hat{\bar{C}}_x^{\gamma_j}(k,2)_i D_{s_i}^{\gamma_j}(k,2)_i \right] \hfill (45)

\displaystyle = K_1 - \frac{1}{\sigma_\nu^2} \sum_{j=1}^{N_2} \sum_{k=-T}^T \left| \hat{\bar{C}}_x^{\gamma_j} (k,2) - D_{s_i}^{\gamma_j} (k,2) \right|^2 \hfill (46)

Substituting the definition of D back into this compact expression leads to

\displaystyle \ln P(\hat{\bar{\boldsymbol{C}}}| H_i) = K_1 - \frac{1}{\sigma_\nu^2} \sum_{j=1}^{N_2} \sum_{k=-T}^T \left| \hat{\bar{C}}_x^{\gamma_j} (k, 2) - \hat{a}_i^2 \bar{C}_{s_i}^{\gamma_j} (k, 2)e^{-i2\pi \gamma_j \hat{t}_{0,i}} \right|^2 \hfill (47)

If we want to maximize the conditional density function over the various signal-type hypotheses indexed by i, we can maximize the log-likelihood function, which leads to minimizing the right term of (47), which cannot be positive. Therefore, our signal-class index estimate is

\displaystyle \hat{i} = \mbox{\rm arg}\, \min_{i \in B} \sum_{i=1}^{N_2} \sum_{k=-T}^T  \left| \hat{\bar{C}}_x^{\gamma_j} (k, 2) - \hat{a}_i^2 \bar{C}_{s_i}^{\gamma_j} (k, 2) e^{-i2\pi \gamma_j \hat{t}_{0,i}} \right|^2 \hfill (48)

Discussion

The order-reduced version of our general problem results in a relatively simple distance-based classifier embodied by Equation (48). This classifier compares all cyclic autocorrelation functions in the catalog of known types with the measured cyclic autocorrelation functions obtained from the data at hand. It can work provided no entries of the catalog set B are scaled versions of each other:

\displaystyle \mathcal{C}_i(2) \neq \mathcal{C}_j(2) \ \ \forall i \neq j, \ \ i, j \in B \hfill (49)

where \displaystyle \mathcal{C}(2) is the set of cumulants restricted to order n=2.

Unfortunately, there are many signals for which the condition (49) does not hold. In particular, all PSK signals with four or more points in their constellation, and where the constellation is doubly symmetric in the complex plane, have identical second-order statistics. This can be referred to as an equivalence class. The signals for which (49) does not hold form a subset of signal types that are equivalent with respect to the lens of second-order limit parameters (time-varying mean and time-varying autocorrelation).

To accurately distinguish between signals that fall into such an equivalence class, we can consider a delay-reduced version of the general problem instead of the order-reduced version (Problem Nine). For this new problem variant, we restrict the range of the n-dimensional delay vector \mathbf{\tau} to the origin and restrict the order n to some integer N_{max} > 2. This thinking leads us to Problem Twelve.

Problem Twelve

Our hypotheses are now given by

\displaystyle H_i: C_x^{\beta_i(n)} (\boldsymbol{0}; n) = a^n C_{s_i}^{\beta_i(n)}(\boldsymbol{0};n) e^{-i 2 \pi \beta_i(n) t_0} \hfill (49)

for i=1, \ldots, L_B, n=1, 2, \ldots, N_{max}, and \beta_i(2) \neq 0. That last constraint just means we don’t want to include the conventional (non-cyclic) autocorrelation (equivalently, we want to exclude the power spectrum). For each value of the order n, let the set of cycle frequencies \beta_i(n) be represented by the \gamma as before,

\displaystyle \beta_i(n) = \left\{ \gamma_{i,n}(j) \right\}_{j=1}^{N_n} \hfill (50)

where N_n is the number of nth-order cycle ferquencies.

The formulation in (49) can be expressed in terms of the specific cycle frequencies \gamma,

\displaystyle H_i = \bar{C}_x^{\gamma_{i,n}(j)} (\boldsymbol{0}; n) = a^n \bar{C}_{s_i}^{\gamma_{i,n}(j)} (\boldsymbol{0}; n) e^{-i 2\pi \gamma_{i,n}(j) t_0} \hfill (51)

where i = 1, \ldots, L_B, n=1, \ldots, N_{max}, and j=1, \ldots, N_n. The variant of Problem Twelve in (51) is again recognized as a known-signal-in-noise problem. Let’s assume that we can find hypothesis-dependent estimates of the amplitude and delay given by \hat{a}_i and \hat{t}_{0,i}.

We next introduce the concept of cumulant warping in which each nth-order reduced-dimension cyclic cumulant is raised to the power 2/n. This attenuates the typical rapid growth of the cyclic-cumulant magnitude with order n, which prevents the cumulants with the largest order N_{max} from dominating the eventual classification statistics (which you haven’t yet seen). The warping means that each considered cyclic cumulant is scaled by a^2 and so all the various involved cumulants are commensurate with signal power.

The next Problem Twelve variant is

\displaystyle H_i: \left[ \bar{C}_x^{\gamma_{i,n}(j)}(\boldsymbol{0};n) \right]^{2/n} = \hat{a}_i^2 \left[ \bar{C}_{s_i}^{\gamma_{i,n}(j)}(\boldsymbol{0}; n) \right]^{2/n} e^{-i4\pi \gamma_{i,n}(j)\hat{t}_{0,i}/n} \hfill (52)

where i = 1, \ldots, L_B, n=1, \ldots, N_{max}, and j=1, \ldots, N_n.

Let’s now (naively) assume that these nonlinearly transformed cyclic cumulants can be estimated with variance independent of both n and \gamma_{i,n}(j)

\displaystyle H_i: \left[ \hat{\bar{C}}_x^{\gamma_{i,n}(j)} (\boldsymbol{0};n) \right]^{2/n} = \hat{a}_i^2 \left[ \bar{C}_{s_i}^{\gamma_{i,n}(j)} (\boldsymbol{0};n) \right]^{2/n} e^{-i4 \gamma_{i,n}(j) \hat{t}_{0,i}/n} + \nu(\gamma_{i,n}(j), n) \hfill (53)

where i = 1, \ldots, L_B, n=1, \ldots, N_{max}, and j=1, \ldots, N_n. Here the random variable \nu models the estimation error, and it has zero mean and a variance independent of n and the cycle frequency \gamma_{i,n}(j).

For convenience, let’s compact the notation

\displaystyle \hat{D}_x^{\gamma_{i,n}(j)}(n) = \left[ \hat{C}_x^{\gamma_{i,n}(j)}(\boldsymbol{0};n) \right]^{2/n} \hfill (54)

\displaystyle D_{s_i}^{\gamma_{i,n}(j)} (n) = \hat{a}_i^2 \left[ \bar{C}_{s_i}^{\gamma_{i,n}(j)} (\boldsymbol{0};n) \right]^{2/n} e^{-i 4 \pi \gamma_{i,n}(j) \hat{t}_{0,i}/n} \hfill (55)

Then our final variant of Problem Twelve is simply

\displaystyle H_i: \hat{D}_x^{\gamma_{i,n}(j)}(n) = D_{s_i}^{\gamma_{i,n}(j)}(n) + \nu(n) \hfill (56)

where i = 1, \ldots, L_B, n=1, \ldots, N_{max}, and j=1, \ldots, N_n.

The interpretation of Problem Twelve is that on the ith hypothesis, a particular (large) set of measured cyclic cumulants is equal to a noisy version of the ideal cyclic cumulants for the signal corresponding to i. The trick is that the cyclic cumulants are severely restricted in terms of the delay vector \boldsymbol{\tau} and the maximum considered order n.

Solution to Problem Twelve

Since (56) is a form of the known-signal-in-noise problem, we can immediately write the solution,

\displaystyle \hat{i} = \mbox{\rm arg}\min_{i \in B} \sum_{n=1}^{N_{max}} \sum_{j=1}^{N_n} \left| \hat{D}_x^{\gamma_{i,n}(j)}(n) - D_{s_i}^{\gamma_{i,n}(j)}(n) \right|^2 \hfill (57)

Or, back substituting for the D variables to express it directly in terms of cyclic cumulants,

\displaystyle \hat{i} = \mbox{\rm arg} \max_{i \in B} \sum_{n=1}^{N_{max}} \sum_{j=1}^{N_n} \left| \left[ \hat{\bar{C}}_x^{\gamma_{i,n}(j)} (\boldsymbol{0}; n) \right]^{2/n} - \hat{a}_i^2 \left[ C_{s_i}^{\gamma_{i,n}(j)} (\boldsymbol{0}; n) \right]^{2/n} e^{-i 4 \pi \gamma_{i,n}(j) \hat{t}_{0,i} /n} \right|^2 \hfill (58)

The interpretation of the solution to Problem Twelve is that the best decision index corresponds to the minimum of the sum of squared distances between the signal’s ideal cyclic cumulants and the corresponding measured cyclic cumulants, or what I call cyclic-cumulant matching.

Reverie

That was a lot. And, I am sorry to say, not one figure has yet been placed into this exposition. We ended up with two likelihood-based (decision theoretic) signal classifiers that are each a major simplification of a much more complete and imposing signal classifier. Along the way, we thought hard about how to capture the probabilistic essence of a manmade communication signal. As I say elsewhere on the CSP Blog, what makes a BPSK signal a BPSK signal? What is BPSKness? I opted to pursue a series of problem statements that focused on the notion that a signal’s unique signalness can be completely captured by the infinite set of all of its cumulants.

In this post we are able to take a lot of time, and a lot of words and equations, to get us to a potentially useful algorithm. Now compare what we did here with what I wrote in My Papers [28], which also presents the order-reduced and delay-reduced classifiers. I simply can’t go through all the intermediate steps in a published paper–there is always a severe page limit relative to the writing and equations required by a complicated R&D algorithm-development project.

Extensions

Power-Constrained Classification

The restriction on distinct catalog elements for both the order-reduced and delay-reduced classifiers can be eased if prior information on the signal power is available.

Order-Reduced Classifier with Power Constraint

Since all catalog elements arise from signal models with unit power, the variable a^2 represents the power of the received signal. The variable we’ve been calling \hat{a}_i^2 is an estimate of the signal power conditioned on hypothesis i–this power estimate is dependent on information contained in the ith hypothesis.

Suppose, though, that we obtain a hypothesis-independent estimate of a denoted by \hat{a}. Then we can use this estimate directly in the classifier, substituting each dependent estimate with the independent one, or we can constrain the decision to be one of the signal-type hypotheses for which the dependent power estimate is close to the independent one,

\displaystyle \hat{i} = \mbox{\rm arg} \max_{i \in B} \sum_{j=1}^{N_2} \sum_{k=-T}^T \left| \hat{\bar{C}}_x^{\gamma_j} (k;2) - \hat{a}_i^2 \bar{C}_{s_i}^{\gamma_j} (k;2) e^{-i 2\pi \gamma_j \hat{t}_{0,i}} \right|^2 \hfill (59)

subject to

\displaystyle \left| \hat{a}^2 - \hat{a}_i^2 \right|^2 < \epsilon \hfill (60)

We choose the maximum likelihood catalog entry subject to a constraint on the estimated power level; estimated power levels that are too far removed from the independent estimate \hat{a}^2 are grounds for exclusion of the corresponding signal hypotheses.

We have increased our tolerance to similar catalog entries for the order-reduced classifier, but it still fails when two catalog entries are actually identical. Such entries would yield identical dependent power levels, rendering the constraint useless to discriminate between them.

Delay-Reduced Classifier with Power Constraint

We can introduce the same constraint to the delay-reduced classifier to accommodate warped cyclic-cumulant catalog entries that are scaled versions of each other,

\displaystyle \hat{i} = \mbox{\rm arg} \min_{i \in B} \sum_{n=1}^{N_{max}} \sum_{j=1}^{N_n} \left| \left[ \hat{C}_x^{\gamma_{i,n(j)}} (\boldsymbol{0};n) \right]^{2/n} - \hat{a}_i^2 \left[ \hat{C}_{s_i}^{\gamma_{i,n}(j)} (\boldsymbol{0};n) \right]^{2/n} e^{-i4\pi \gamma_{i,n}(j) \hat{t}_{0,i}/n} \right|^2 \hfill (61)

subject to

\displaystyle \left| \hat{a}^2 - \hat{a}_i^2 \right|^2 < \epsilon \hfill (62)

Discussion

The use of the power constraint in the two derived classifiers might seem strange. Why not simply replace the hypothesis-dependent power estimates \hat{a}_i^2 with the hypothesis-independent power estimate \hat{a}^2 and compute the likelihoods? This will work just fine, I’ll wager, and is essentially equivalent to the constraint. However, that method will not be available when we generalize the problem to multiple cochannel signals because we can only obtain a hypothesis-independent estimate of the sum of the involved signal powers, not estimates of their individual power levels. That may be a bit strong of a statement. When the signals have different spectral shapes and/or their relative center frequencies are significantly different (so that they only overlap partially spectrally), you might be able to get at power estimates for each signal. But when the signals have similar spectral shapes and their centers are closely spaced, you will just observe a smooth lump of energy in the power spectrum. How much does each signal contribute?

If I give you the sum of two students’ test scores (say, 152), can you tell me the two scores?

Multiple Cochannel Signals

Our cochannel signal model is a simple generalization of our original model (7),

\displaystyle x(t) = \sum_{m=1}^M a_m s_m(t-t_m) + w(t) \hfill (63)

The cochannel model for which the s_m(t) are disjoint in either time or frequency or both reduces to a set of single-signal-in-noise problems of the sort previously considered. This can be accomplished by segmenting the input data in time and frequency using algorithms such as that described in the post on automatic spectral segmentation. So we assume here that the M signals overlap in time and at least partially in frequency so that they cannot be separated by time gating and/or linear time-invariant filtering.

We now develop a sequence of hypothesis-testing problems that parallels the previous single-signal development. In this development, we assume the number of signals M is known or estimated prior to signal classification. (In practice, we don’t invoke this assumption. The number of signals is estimated jointly with blind estimation of all the cycle frequencies exhibited by the data.)

In the single-signal case, a hypothesis consisted of the identification of a probability structure for the signal,

\displaystyle H_i: \ \ \ s(t) \sim \mathcal{P}_i, \ \ \ i \in B. \hfill (64)

We now identify a hypothesis with the selection of a probability structure that incorporates the individual probability structures for each of the M constituent signals. Let’s introduce some vector notation here,

\displaystyle x(t) = \boldsymbol{s}(t) \boldsymbol{a}^\top + w(t) \hfill (65)

where \displaystyle \boldsymbol{s}(t) a vector of signals

\displaystyle \boldsymbol{s}(t) = \left[ s_1(t-t_1) \ s_2(t-t_2) \ldots s_M(t-t_M)\right] \hfill (66)

and \boldsymbol{a} is a vector of amplitudes

\displaystyle \boldsymbol{a} = [a_1 \  a_2 \ldots a_M] \hfill (67)

and \top denotes a matrix transpose.

We say that \boldsymbol{s}(t) has probability structure \boldsymbol{\mathcal{P}}_{\boldsymbol{i}},

\displaystyle \boldsymbol{s}(t) \sim \boldsymbol{\mathcal{P}}_{\boldsymbol{i}} \hfill (68)

\displaystyle \boldsymbol{i} = [i_1 \ i_2 \ldots i_M] \in B^M = \underbrace{B \times B \times \ldots \times B}_{M} \hfill (69)

\displaystyle \boldsymbol{\mathcal{P}}_{\boldsymbol{i}} = \left[\mathcal{P}_{i_1} \ \mathcal{P}_{i_2} \ldots \mathcal{P}_{i_M} \right] \hfill (70)

\displaystyle s_j(t) \sim \mathcal{P}_{i_j}, \ \ \ j=1, 2, \ldots, M. \hfill (71)

Problem Thirteen

Our first cochannel-signal hypothesis is simply

\displaystyle H_{\boldsymbol{i}}: \ \ \ \boldsymbol{s}(t) \sim \boldsymbol{\mathcal{P}}_{\boldsymbol{i}}, \ \ \ \boldsymbol{i} \in B^M \hfill (72)

We can form equivalent problem statements using either moments or cumulants,

\displaystyle H_{\boldsymbol{i}}: \ \ \ \boldsymbol{s}(t) \sim \boldsymbol{\mathcal{M}}_{\boldsymbol{i}}, \ \ \ \boldsymbol{i} \in B^M \hfill (73)

\displaystyle H_{\boldsymbol{i}}:\ \ \ \boldsymbol{s}(t) \sim \boldsymbol{\mathcal{C}}_{\boldsymbol{i}}, \ \ \ \boldsymbol{i} \in B^M \hfill (74)

The cumulant formulation is preferred because the cumulant structure of x(t), which contains multiple statistically independent signals, is vastly simpler than the moment structure. This is a consequence of the additive (“cumulants accumulate”) nature of cumulants, which nature the moments do not possess,

\displaystyle \boldsymbol{s}(t) \sim \boldsymbol{\mathcal{C}}_{\boldsymbol{i}} \Longleftrightarrow s_m(t) \sim \mathcal{C}_{i_m}, \ \mbox{\rm with}\  m = 1, 2, \ldots, M. \hfill (75)

Since the signals are presumed to be independent, the cumulants for the sum signal x(t) are just the sums of the individual-signal cumulants,

\displaystyle C_x(t,\boldsymbol{\tau};n) = \sum_{m=1}^M D_{s_m}(t, \boldsymbol{\tau}; n) + C_w(t, \boldsymbol{\tau}; n) \hfill (76)

Let’s use the notation \displaystyle x(t) \sim \boldsymbol{\mathcal{C}}_{\boldsymbol{i}+w} to indicate that \boldsymbol{s}(t) \sim \boldsymbol{\mathcal{C}}_{\boldsymbol{i}} for our data signal x(t). “The data is distributed as a particular (infinite) set of (infinite) sets of cumulants indexed by \boldsymbol{i} = [i_1\ i_2 \ldots i_M].

Problem Fourteen

In this problem, our hypotheses take the form of “the data on hypothesis i is distributed as a collection of cumulant sets, one for each signal that is present,”

\displaystyle H_{\boldsymbol i}: \ \ \ x(t) \sim \boldsymbol{\mathcal{C}}_{\boldsymbol{i}+w}, \ \ \ \boldsymbol{i} \in B^M \hfill (77)

This is equivalent to the following formulation in terms of temporal cumulant functions,

\displaystyle H_{\boldsymbol{i}}: \ \ \ C_x(t, \boldsymbol{\tau};n) = \sum_{m=1}^M D_{s_{m,i_m}}(t, \boldsymbol{\tau};n) + C_w(t, \boldsymbol{\tau}; n) \hfill (78)

for \displaystyle \boldsymbol{i} = [i_1\ i_2 \ldots i_M] \in B^M and n=1, 2, \ldots, \boldsymbol{\tau} \in \mathbb{R}^n. Recall that in (78), the D_s(\cdot) function is a more compact way of writing the temporal cumulant,

\displaystyle D_{s_m} (t, \boldsymbol{\tau};n) = a_m^n C_{s_m}(t-t_m, \boldsymbol{\tau};n) \hfill (79)

Now we need to represent each time-varying cumulant in terms of the cyclic cumulants. We take C_{s_{m,i_m}}(\cdot) to be the cumulant for the mth signal when that signal takes on catalog label i_m. Therefore the cycle frequencies for s_{m,i_m} depend on i_m. However, the most interesting and difficult case corresponds to cycle frequencies being dependent on m only, which is similar to our philosophy in the single-signal case. Recall there we deemphasized the situation where each catalog entry depends on the numerical values of the cycle frequencies–so you could detect the cycle frequencies, then simply look up which catalog entry corresponds to those particular numerical values.

Using the Fourier-series decomposition of the time-varying cumulant,

\displaystyle C_{s_{m,i_m}} (t, \boldsymbol{\tau};n) = \sum_{j=1}^{N_{n,m}} D_{m,i_m}^{\gamma_j(n,m)} (\boldsymbol{\tau};n) e^{i 2 \pi \gamma_j(n,m) t} \hfill (80)

We reformulate the problem using the D notation as follows

\displaystyle H_{\boldsymbol{i}}: C_x^{\gamma_j(n,m)}(\boldsymbol{\tau};n) \sum_m D_{s_{m,i_m}}^{\gamma_j(n,m)}(\boldsymbol{\tau};n) + C_w^{\gamma_j(n,m)} (\boldsymbol{\tau};n) \hfill (81)

where \displaystyle \boldsymbol{i}\in \prod_{m=1}^M B, \displaystyle n=1, 2, \ldots, \boldsymbol{\tau} \in \mathbb{R}_n, j=1, 2, \ldots, N_{n,m}, m=1, 2, \ldots, M.

If, as before, we agree to avoid the “second-order stationary” statistic, that is, we don’t include \gamma_j(2,m) = 0 in our sum, we obtain the next variation on the cochannel problem formulation

\displaystyle H_{\boldsymbol{i}}: C_x^{\gamma_j(n,m)}(\boldsymbol{\tau};n) \sum_m D_{s_{m,i_m}}^{\gamma_j(n,m)}(\boldsymbol{\tau};n) \hfill (82)

where \displaystyle \boldsymbol{i}\in \prod_{m=1}^ M B, n=1,2, \ldots, \boldsymbol{\tau} \in \mathbb{R}^n, j=1, 2, \ldots, N_{n,m}, and \displaystyle \gamma_j(2,m) = 0.

Because we’ve assumed the noise w(\cdot) is stationary and Gaussian, the exclusion of the “second-order stationary” cycle frequency of zero means that no cyclic cumulants for w(\cdot) appear in the problem formulation.

Now we assume that all nonzero cycle frequencies are unique to one signal, and we continue to avoid the second-order cycle frequency of zero. This leads to the problem formulation

\displaystyle H_{\boldsymbol{i}}:  C_x^{\gamma_j(n,m)} (\boldsymbol{\tau};n) = D_{s_{m,i_m}}^{\gamma_j(n,m)} (\boldsymbol{\tau};n) \hfill (83)

where \displaystyle \boldsymbol{i}\in \prod_{m=1}^M B, \displaystyle n=1, 2, \ldots, \boldsymbol{\tau} \in \mathbb{R}^n, j=1, 2, \ldots, N{n,m}, m=1, 2, \ldots, M, and \gamma_j(2,m) \neq 0.

The case in which some signals share some or all of their cycle frequencies is treated another time.

Using the spelled-out expression for the D(\cdot) function,

\displaystyle D_{s_{m,i_m}}^{\gamma_j(n,m)} (\boldsymbol{\tau};n) = a_m^n C_{s_{m,i_m}}^{\gamma_j(n,m)}(\boldsymbol{\tau};n) e^{-i 2\pi \gamma_j(n,m) t_m} \hfill (84)

and naively assuming that we can estimate the C_x^\gamma(\cdot) functions with variance independent of n, \boldsymbol{\tau}, and \gamma, we formulate the next problem variant as

\displaystyle H_{\boldsymbol{i}}: \hat{C}_x^{\gamma_j(n,m)} (\boldsymbol{\tau};n) = a_m^n C_{s_{m,i_m}}^{\gamma_j(n,m)} (\boldsymbol{\tau};n) e^{-i2\pi \gamma_j(n,m) t_m}  + v(\boldsymbol{\tau};n, \gamma_j) \hfill (85)

Because the cumulant can grow rapidly with increasing n, we propose the nonlinearly transformed version of (84) where each nth-order cyclic cumulant is raised to the power 2/n,

\displaystyle H_{\boldsymbol{i}}:  \left[ \hat{C}_x^{\gamma_j(n,m)}(\boldsymbol{\tau};n) \right]^{2/n} = a_m^2 \left[ C_{s_{m,i_m}}^{\gamma_j(n,m)} (\boldsymbol{\tau};n) \right]^{2/n} e^{-i 4 \pi \gamma_j(n,m) t_m/n} + v(\boldsymbol{\tau}, n, \gamma_j) \hfill (86)

A final helpful modification is to replace the cyclic cumulants with the reduced-dimension cyclic-cumulants,

\displaystyle H_{\boldsymbol{i}}:  \left[ \hat{\bar{C}}_x^{\gamma_j(n,m)} (\boldsymbol{\tau};n) \right]^{2/n} = a_m^2 \left[ \bar{C}_{s_{m,i_m}}^{\gamma_j(n,m)} (\boldsymbol{\tau};n) \right]^{2/n} e^{-i 4 \pi \gamma_j(n,m) t_m/n} + v(\boldsymbol{\tau}, n, \gamma_j). \hfill (87)

As before, we will derive the order-reduced and the delay-reduced classifiers. Unlike before, however, we do not assume we have a prior estimate of a_m^2, only of \sum_{m=1}^M a_m^2 = a^2, the total combined power of the several involved cochannel signals. This we can do with automatic spectral segmentation and simple power-signal measurements. On the other hand, due to the disjoint nature of the cyclic cumulants (we have imposed by assumption), we can obtain hypothesis-dependent estimates of the individual-signal power levels \hat{a}_{m,i_m}. This in turn means we can compare the prior estimate of a^2 with the sum of the hypothesis-dependent power levels, producing a goodness-of-fit relative to power for any particular multiple-signal hypothesis.

Multiple-Signal Order-Reduced Classifier

Restricting the maximum cumulant order to two and computing the likelihoods for the problem variant (87) leads to the multiple-signal ORC,

\displaystyle \hat{\boldsymbol{i}} = \mbox{\rm arg} \min_{\boldsymbol{i}=[i_1, i_2, \ldots, i_M]} \sum_{m=1}^M \sum_{\stackrel{j=1}{\gamma\neq 0}} \left| \hat{\bar{C}}_x^{\gamma_j(2,m)} (\tau) - \hat{a}_{m,i_m}^2 \bar{C}_{s_{m,i_m}}^{\gamma_j(2.m)} (\tau) e^{-i 2\pi \gamma_j(2,m)\hat{t}_{m,i_m}} \right|^2 \hfill (88) 

subject to

\displaystyle \left| \sum_{m=1}^M \hat{a}_{m,i_m}^2 - \hat{a}^2 \right| < \epsilon \hfill (89)

which is called the sum-of-powers (SOP) constraint. The SOP constraint says to pick the catalog indices for the elements of B such that the sum of the power levels implied by the catalog selections and measured features is close to the measurable power of the non-noise input-signal component.

The multiple-signal order-reduced classifier will fail when some of the catalog entries (the \displaystyle \bar{C}_s^\gamma(\tau) functions, which are the cyclic autocorrelation functions) are identical. This can happen, for instance, for square-constellation digital QAM signals such as QPSK, 16QAM, 64QAM, etc. Therefore we will also need the multiple-signal version of the delay-reduced classifier.

Multiple-Signal Delay-Reduced Classifier

Here we restrict the maximum cumulant order to N_{max} and the cumulant delay vector to the origin, \boldsymbol{\tau} = \boldsymbol{0}. The hypotheses are

\displaystyle H_{\boldsymbol{i}}: \left[ \hat{\bar{C}}_x^{\gamma_j(n,m)} (\boldsymbol{0};n) \right]^{2/n} = \hat{a}_{m,i_m}^2 \left[ \bar{C}_{s_{m,i_m}}^{\gamma_j(n,m)} (\boldsymbol{0};n) \right]^{2/n} e^{-i 4\pi \gamma_j(n,m) \hat{t}_{m,i_m}/n} + v(n, \gamma_j) \hfill (90)

subject to the SOP constraint,

\displaystyle \left| \sum_{m=1}^M \hat{a}_{m,i_m}^2 - \hat{a}^2 \right| < \epsilon \hfill (91)

If I’ve done my job well here, you should be able to write down the expression for the multiple-signal DRC decision-index-vector estimate. I’ll leave you to it!

Other Extensions

The extension of these ideas and results to complex-valued signals is straightforward. The notation gets even more complicated because we have to consider all relevant conjugation configurations–how many factors are conjugated in the second-, fourth-, and sixth-order cumulants used by the delay-reduced classifier? Symmetry results are useful there.

The cases where two or more cochannel signals share (1) some of their cycle frequencies or (2) all of their cycle frequencies appear in practical settings and are good candidates for generalizing the results. Get to it! Because I’m tired of this post and you are probably tired of it too.

Blind Versus Non-Blind Modulation Recognition

In practice, we often do not know any of the involved signals’ cycle frequencies or modulation types: we are blind to these values. To apply the order-reduced and delay-reduced classifiers, you need to know the cycle frequencies so you can compute the relevant cyclic cumulants from the available data record. Therefore you must find a way to blindly estimate the cycle frequencies. I assert, with the benefit of substantial practical experience implementing these classifiers and applying them to millions of simulated and captured signals, that this is possible.

Closing Remarks

So that’s it. This long post is meant to provide a detailed behind-the-scenes look at the admittedy coarse sketch of results shown in My Papers [28].

Looking back at [28], published in 2000, I find that I did not number the equations. I am not happy with my younger self. But since I can’t easily refer to some of the equations in [28] here, it gives me an opportunity to include a figure, finally, in this post.

Figure 1 shows the order-reduced and delay-reduced classifier expressions as I original published them in [28]. There I’ve reexpressed this post’s (48) in the frequency domain so that the ORC is clearly a spectral-correlation-matching algorithm. I figured the spectral correlation function was likely more familiar to more readers than the cyclic autocorrelation or its cousin the second-order cyclic temporal cumulant.

Figure 1. The order-reduced and delay-reduced classifiers as presented in their original form in My Papers [28].

Have I achieved my goal of showing that complex CSP-based algorithms may appear to be conjured from nothing by some sorcery aided by a dark lore, but are actually often founded on rather dry mathematical facts and logical principled decision-making approaches such as decision theory?

Only you can be the judge of that.

Author: Chad Spooner

I'm a signal processing researcher specializing in cyclostationary signal processing (CSP) for communication signals. I hope to use this blog to help others with their cyclo-projects and to learn more about how CSP is being used and extended worldwide.

Leave a Comment, Ask a Question, or Point out an Error

Discover more from Cyclostationary Signal Processing

Subscribe now to keep reading and get access to the full archive.

Continue reading