In this post, I discuss a signal-processing algorithm that has almost nothing to do with cyclostationary signal processing. Almost. The topic is automated spectral segmentation, which I also call band-of-interest (BOI) detection. When attempting to perform automatic radio-frequency scene analysis (RFSA), we may be confronted with a data block that contains multiple signals in a large number of distinct frequency subbands. Moreover, these signals may be turning on an off within the data block. To apply our cyclostationary signal processing tools effectively, we would like to isolate these signals in time and frequency to the greatest extent possible using linear time-invariant filtering (for separating in the frequency dimension) and time-gating (for separating in the time dimension). Then the isolated signal components can be processed serially.

It is very important to remember that even perfect spectral and temporal segmentation will not solve the cochannel-signal problem. It is perfectly possible that an isolated subband will contain more that one cochannel signal.

The basics of my BOI-detection approach are published in a 2007 conference paper (My Papers [32]). I’ll describe this basic approach, illustrate it with examples relevant to RFSA, and also provide a few extensions of interest, including one that relates to cyclostationary signal processing.

### The Problem

The efficacy of our cyclostationary signal processing (CSP) algorithms for any particular data record depends on several factors, including the noise power, interference power, processing block length, and signal type. In general, it is better to apply CSP algorithms to a signal that has been filtered so that adjacent-channel signals and out-of-band noise are minimized. This just means that we want to minimize the “total SINR” of the data prior to CSP processing. By “total SINR” I simply mean the ratio of the signal power to the power of everything else in the data record, whether it has frequency content in the signal’s band or not.

This aim isn’t so hard to achieve when the input data consists of a single signal in white Gaussian noise (WGN). Then we need only apply a linear time-invariant filter to the data, where the center frequency and bandwidth of the filter match the center frequency and bandwidth of the signal. For example, consider a QPSK signal in noise. The PSD looks like this:

If we apply a simple Fourier-based linear time-invariant filter to that QPSK data, where the center frequency is zero and the width is Hz, we remove almost all of the out-of-band noise and retain the signal:

Removing the out-of-band noise won’t help much with estimating the spectral correlation function for the QPSK signal itself, because that will depend on the spectral components of the data that do in fact lie in the signal band. On the other hand, the SCF for frequencies that lie well outside the signal’s band is ideally zero, and will be easier to estimate when the out-of-band noise is removed.

But what about the temporal parameters like the cyclic autocorrelation and the cyclic cumulants? Well, they will be more sensitive to the presence of the out-of-band noise, because each point estimate takes into account all values of the input data block, so it is important to remove as much of it as possible prior to applying the estimator for the time-domain parameter of interest.

So we have reasons for wanting to apply simple signal processing methods to simple RF scenes like the QPSK signal above: increasing the total SNR improves estimator performance for fixed estimator parameters such as block length, spectral resolution etc. And increasing the total SNR means effectively locating the signal in the frequency domain and then applying a bandpass filter to that location.

The problems for our estimators get worse, however, when the scenarios aren’t simple.

#### Densely Populated RF Bands

In RF scene analysis, the goal of the signal processing algorithm is to blindly detect, classify, and characterize all signals present in the scene. No matter what the scene contains. Here is a sequence of PSDs for scenes of interest, where the scene complexity–measured by things like the number of signals present, the SINRs, and the deviation from ideal of the noise floor–is increasing:

By the time we get to the final few examples, we see collection-system filter roll-offs on the edges of the band, non-flat noise floors, high dynamic ranges in signal power and signal bandwidth, and closely spaced signals in frequency.

#### Distortion of Estimated Cumulants Due to Time-Gating

There is another reason to isolate a signal in the time and frequency domains prior to processing with CSP: distortion of sets of cyclic cumulants caused by partial signal temporal occupancy within a processed data block. (Awkwardly put, I know.)

Suppose we have a signal that exists only during a subset of the time instants corresponding to some sampled-data record. Let the subset occupy a fixed fraction of the data record, which itself spans the continuous-time interval ,

Here and is stationary white Gaussian noise. We use because it allows us to fix the *fraction* of time that the signal is on, relative to the interval , so we can talk sensibly about the values of estimators as . For example, if and , then and we have our usual case where the signal persists over the entire observation (integration) interval. On the other hand, if and , then and the signal persists only over half of the observation interval.

Now we know that the temporal cumulant function is a particular sum of products of lower-order temporal moment functions,

(See the post on cyclic cumulants for more information on this notation.)

Now the cumulant for is the sum of the cumulant for and , and the latter is zero for , so we focus on the cumulant for . We’ll need the various lower-order cyclic moments (Fourier components of the time-varying moment ),

Or, using and ,

In the limit as , we have

So the signal that persists only over a fraction of the observation interval is a scaled version of the signal that persists over the entire interval, and the scaling factor is quite reasonably .

This implies that the temporal moment function is also scaled by ,

Keeping in mind that for higher orders the cumulant for is that for , we have

Unlike the temporal moment, the temporal cumulant is **not** simply a scaled version of that fully persistent signal–it is warped or distorted by the inclusion of , which is dependent on the partition size unless (the case in which the signal is persistent over the whole block).

So, processing involving linear combinations of multiple orders of moments is insensitive to , but processing involving cumulants is not, because that involves *different products of moments*. Since my modulation recognition algorithm is cumulant-based (see here, but more coming in a future post as well), this is an actual problem (for me).

All of this implies that we need to include an *interval-of-interest (IOI) estimator* as well as a band-of-interest estimator in our spectral-segmentation algorithm. That way our calculations will always involve once we estimate the IOI and time-gate the signal component accordingly.

### Potential Solutions

A few others have created automatic spectral segmentation algorithms, or the equivalent, such as automatic noise-floor estimators that work for dense (crowded) RF scenes. One is The Literature [R96], which applies morphological image-processing operators (like erosion and dilation) to spectral data. Others, like The Literature [R97] advocate sophisticated high-resolution/high-accuracy spectrum estimators such as the multi-tapering method of The Literature [R98]. I think we do not need high spectrum-estimate accuracy here; the goal is to ** estimate the frequency intervals **for which non-noise energy exists, not to obtain a super accurate spectrum estimate for all frequencies in the analysis band. So I’m of the opinion that high-performance spectrum estimators are out of place in the general band-of-interest estimation problem.

### My Solution: Multiresolution Adaptive Spectrum Estimation

I outline the solution in My Papers [32]. The idea is to recursively apply a low-cost simple non-parametric power spectrum estimator to the input data record, which is assumed to contain complex samples and correspond to the complex envelope of some RF data (like the various spectra above).

Why “recursively?” In some cases, a wideband RF scene may contain several signals that are closely spaced in frequency. To see that there is a small gap between them requires that the spectral resolution of the spectrum estimator lies in a particular range. On the other hand, the scene may also (or only) contain a wideband signal with significant spectral variation across its bandwidth. In this latter case, there is only one band of interest, but the variability of a high-resolution spectrum estimate can make it appear that it is composed of multiple closely spaced bands of interest.

In other words, the required spectral resolution depends on the number and spacing of the signals in the scene, and we don’t want to make many assumptions about that.

What we do instead is to start out by applying a coarse spectral resolution to the data record, find the minima to estimate the noise floor, and declare everything else as bands of interest. Then, we apply the same algorithm to each of the detected BOIs, after filtering, frequency shifting, and decimating each one. This process will reveal which of the coarse lumps from the first spectrum estimate are really composed of multiple closely spaced BOIs and which are not.

I’ve found that this is a tricky idea to get across, so before we outline the mathematics of the algorithm, I want to further describe the multi-resolution idea with some pictures.

Consider the following multi-signal RF scene:

This scene contains individual signals, none of which overlap in frequency. This first PSD estimate for our -signal scene used the frequency-smoothing method and a smoothing-window width (spectral resolution) of Hz ( Hz in normalized frequency). With this resolution, our eyes can count bands of interest, and we suspect the ones near kHz may contain several more. But the resolution is too coarse to reveal them in this picture.

More generally, let’s look at a sequence of PSD estimates for this scenario:

In this sequence of PSD estimates, the spectral resolution is decreasing from top to bottom (the smoothing-window width in the FSM is increasing). The bands of interest are detected using only a single PSD estimate–no recursive processing is allowed. This shows that basic spectral segmentation algorithms fail when the PSD spectral resolution is not well-matched to the spectral variation exhibited by the underlying signals in the data.

The idea behind the automatic multi-resolution spectral segmentation algorithm is recursively applying a standard non-parametric PSD estimator (like the FSM), which I also call *telescoping* spectral estimation. For this -signal example, telescoping spectral segmentation is illustrated by the following diagram:

When used on our -signal example, we obtain the following sequence of results:

The indicated spectral resolution applies to the initial application of the non-parametric PSD estimator. Subsequent applications, as the telescoping action proceeds, use an effective resolution that has the same relation to the local sampling rate that the initial resolution has to the original sampling rate of kHz.

Here we see that the ability to correctly detect the bands of interest and to estimate the noise floor is only weakly dependent on the choice of spectral resolution.

### Algorithm

First we provide a fixed-resolution algorithm that finds bands of interest by defining them as arbitrary-shaped bumps that rise up from the noise floor. The multi-resolution version will follow.

In this algorithm statement, the th BOI has endpoints given by .

- Obtain data samples , .
- Specify resolution (normalized frequency units).
- Let . This is the width of the smoothing window in frequency bins.
- Compute spectrum estimate with resolution product , denote by , .
- is the minimum of for .
- Denote by the noise-floor threshold .
- For each :
- Set .
- For each , if and , set . If and , set and increment by .
- If , then . (Handles BOI at left edge of PSD.)
- If , then . (Handles BOI at right edge of PSD.)
- Estimate noise floor by computing average value of PSD for all frequency points not assigned to a BOI interval.
- Estimate power for BOI by summing PSD values for all frequency points assigned to BOI .
- Remove weak BOIs, resulting in BOIs.
- BOI intervals are , .
- Stop.

We build on this fixed-resolution algorithm to create the recursive multi-resolution algorithm in the following way:

- Obtain data samples , .
- Specify nominal resolution .
- Set , , and .
- Set .
- Set .
- Set .
- If , set and Goto Step 13.
- Extract the BOI data from that corresponds to . Bandpass filter, downconvert, and decimate this data so that its fractional bandwidth is maximized and it is centered at zero frequency. The resulting data is .
- Estimate the PSD for .
- Use the fixed-resolution BOI detection algorithm with to find BOIs in . Let the number of detected BOIs be .
- Assign detected-BOI values to , .
- If and , set , indicating the BOI is completed, else set for .
- .
- .
- If , then Goto Step 7.
- .
- Set .
- If for any , Goto Step 6 else Stop.

When we apply this algorithm to the -BOI data file in a Monte Carlo simulation, we obtain the following performance graphs:

Here the total signal power is about dB, and the total power of the noise varies from dB to dB.

### Examples

Let’s apply the algorithm to the RF scenes I showed in the first several figures of the post. In all cases, the number of processed samples is (regardless of the sampling rate) and the spectral resolution at each stage of recursion is times the effective sampling frequency at that stage.

### Extensions

#### Colored Noise

Several of the RF scenes above possess noise floors that are not flat throughout the analysis bandwidth. In particular, they have noise floors that decrease rather rapidly at the edges of the band. This is typically due to the action of a receiving system, which applies filtering prior to sampling. The applied filter bandwidth is smaller than the complex sampling rate, so that the band edges end up containing little energy. The noise is therefore not white (relative to the sampling rate), and so it is typically called colored. There are other deviations from white noise as well, such as tilted noise floors or those that contain ripples of various sorts. The most common is the filter-roll-off band-edge coloration, and the band-of-interest detector can handle that at the expense of increased computational cost.

The idea is to apply the basic algorithm to the data but restrict it to operate on a subband such as . Keep a count of the number of bands that are detected for each of several values of . Retain only those bands that correspond to the maximum. When , the stated algorithm above will return a single band of interest, which is nearly the entire analysis band, because it will identify the noise floor as the value of the spectrum estimate very near . In other words, the entire spectrum looks like one big band of interest having a very small noise floor. On the other hand, when is such that only the flat part of the noise floor is considered, it will operate effectively and return all the true bands of interest that actually exist in the flat part, which is greater than or equal to one.

#### Interval-of-Interest Estimation

To perform the localization of signal components in the time domain, which I have called *interval-of-interest estimation*, I simply apply the algorithm developed for the band-of-interest estimation to the time-domain data itself. This is done for each detected band of interest separately. Interval-of-interest estimation is just the time-frequency dual algorithm to band-of-interest estimation.

#### Multi-Resolution Spectral Correlation Estimation for Targeted-Signal BOI Detection

The algorithm as discussed so far employs only power spectrum estimation–no cyclostationarity! A CSP extension is to replace the spectrum estimates with spectral correlation estimates. In fact, the entire algorithm above can be viewed as exploiting non-conjugate spectral correlation at a cycle frequency of . Why not choose a different value for ? Such as one that corresponds to a cycle frequency for a signal that is being searched for? Then all bands that do not possess the desired cyclostationarity are ignored and are not processed, potentially saving significant computation.

There are a few complications, but the algorithm works for non-zero . Let’s apply it to each of the scenes above, tailoring the choice for according to the signals that we know exist in the scenes.

First, the simple single-signal scene involving QPSK. It has a symbol rate of , so we use that as the cycle frequency:

Next up is the hybrid FH/DSSS signal, which possesses many cycle frequencies, but a prominent one is . Using this cycle frequency results in the following detection:

Notice here there is a single detected band of interest rather than the two distinct bands produced by the main algorithm.

Next is the scene with the narrowband PSK/QAM signals. The leftmost signal in this scene is BPSK, and the fourth, seventh, etc. signals are also BPSK. The second, fifth, eighth, etc. are QPSK, and the remainder are SQPSK. The bandwidth of the signals is monotonically decreasing as the signal index moves from left to right, so that every signal has a unique symbol rate, which means a unique cycle frequency. We choose the bit rate of the third BPSK signal, the seventh signal from the left, which is . The result is:

Now we turn to the scene containing an FH signal with intermittent wideband interferers. We choose the symbol rate of one of the interferers (situated near Hz), which is :

Next is the scene with a large number of signals of mixed bandwidth. We choose the symbol rate for the two signals between and , which is , and this rate is not a cycle frequency for any of the other signals:

Turning now to captured data, we consider the PCS scene and choose the chip rate for CDMA, which is 1.2288 MHz. This cycle frequency is exhibited by the CDMA signals at and MHz, as well as by the CDMA-EVDO signal at MHz:

When we switch the cycle frequency to the WCDMA chip rate of MHz, we obtain the following bands of interest:

which makes sense because the two signals near MHz are indeed WCDMA.

The broadcast FM signals can exhibit a cycle frequency of kHz due to their digital components:

As usual, if you find any errors in the post, please let me know in the comments. Question are welcome, as always. I’m also very interested in knowing about any papers in the literature that present solutions to this general spectral segmentation problem.

This is something I’ve struggled with a great deal myself. Developing an algorithm for real-world wideband signals of arbitrary sample rate is a major pain.

Most approaches I tried worked well at specific resolutions, but fail miserably at others.

I’ll try your approach and see how it compares to one I’ve made myself (recursively finds the biggest carrier until some SNR threshold is failed).

It would be interesting to see if this could be used to detect transponders – I have signals with multiple transponders – a PITA to analyse.

(Some of your LaTeX math is displaying as raw string)

Thanks! I found two in the second algorithm statement…

How do you avoid missing smaller carriers in the rough FFT at the start?

So an isolated carrier who’s power contribution will get lost compared to the noise in the FFT bin.

Or does that dictate the resolution of our FFTs.

I do miss weaker signals. At some point, the power of a weak signal is comparable to the variance (fluctuations) in the spectral estimates, and it is very hard to tell the difference. So the basic algorithm definitely favors stronger signals. And that gets me pretty far. Your comment reminds me of an extension to the algorithm that I didn’t talk about, but might in an update. That is “white-space detection”, where we apply the algorithm as stated, but at the end, the bands of interest are defined as the spaces in-between the detected bands. We can view these as

possiblynoise-only bands. How to really tell? One way is to then apply CSP (for example, the SSCA) to these putative white spaces to see if there are cyclostationary signals lurking there. Computationally expensive though.I agree it is a difficult problem when the generality we’re aiming for is so high. And the algorithmic approach I’ve sketched is quite useful for me, but is far from perfect. There are so many variations on RF scenes that it seems impossible to construct an algorithm that works flawlessly on all of them. In particular, since the algorithm works by processing multiple (data-adaptive) spectrum estimates, it is in the end energy detection, and weaker signals are missed.

Could you expand on the complications in the last section (matching for specific cyclic frequencies)?

I’m getting lots of artifacts on real data (with many signals), where there isn’t any carrier. Could be due to carrier spacing relating to the cyclic frequency [X(f+a/2) X*(f-a/2)] creating energy?

I’m using the FSM method.

It works well for my simulated data.

The complications arise from my definition of a band of interest as an arbitrary shaped bump between two flat spots. The flat spots correspond to noise-only intervals in the analysis band, and we can usually count on this because the noise is approximately flat across the entire band. So we can “drill down” until we hit the bedrock of the noise floor.

When using spectral correlation to find BOIs, for spectral regions that don’t exhibit the chosen cycle frequency, there is no “flat spot”, because the spectral correlation there is ideally zero. So it becomes harder to tell where the edges of a band of interest are.

Another problem is band-of-interest width estimation. Using the PSD, it is possible to estimate the band-of-interest width as the occupied bandwidth. Subsequent extraction of the band will generally capture almost all the signal energy, so downstream processing works. But with the spectral correlation function, recall that the width of the function can be quite a bit smaller than the occupied bandwidth. For example, for digital QAM/PSK with square-root raised-cosine pulses having roll-off and symbol rate , the occupied bandwidth is , but the width of the symbol-rate spectral correlation feature is just .

But what is the cycle frequency here, and how far apart are distinct signals? It is hard to reply without specifics. Are you using the coherence or the spectral correlation to try to detect BOIs? I also find artifacts are common if there are very narrowband signals present (like tones). These create havoc especially for the coherence function.

I’ve found that I can determine all symbol rates of interest present in a wideband signal using FFT(|x|) given enough data. I was using these as alpha. The nice thing about this is that I know what bandwidths I’m looking for (i.e. symbol rate plus some roll off).

Worst case scenario is for signals to have their 3dB bandwidths touching. Tones (i.e. satellite beacons) are causing havoc yes; but so are carriers of different bandwidths that just happen to be a/2 apart.

I’m straight up just using X(f+a/2) X*(f-a/2). So I suppose that makes it the spectral correlation. Should I be using coherence?

I’m intrigued by the use of the absolute value… contains various nonlinearities in one shot.

Ah. Suppose is a symbol rate for some signal in the data. If there are two signals with center frequencies separated by , then spectral correlation may be present due to some correlation between the two different signals. In other words, suppose the two signals are statistically independent (and therefore uncorrelated). Then if you pick frequencies in the region that is about halfway between the two carriers, and the cycle frequency , you should NOT see significant non-conjugate spectral correlation. However, if the two signals have some information in common, you very well might see spectral correlation. For example, if the two signals are transferring independent user data streams, but use the same periodically repeated training sequence (for, say, channel equalization) or some other in-common data, then you can see spectral correlation. So that would definitely confuse an automatic spectral segmentation algorithm using or a cycle detector.

Coherence, as opposed to spectral correlation, won’t help in this case, because there really is spectral correlation there, so the coherence will reflect that too. In general, I tend to use coherence as a detection statistic for cycle-frequency detection (as opposed to the cycle detectors) because of the relative ease of threshold computation. (A post on that is forthcoming, but see the comments on the SSCA post.)

I come back to this problem once in a while.

Two thoughts that I have (and I lack the knowledge to answer them).

Seems we humans would quite easily be able to separate between signals in spectrum. Is there a way to transform it to an AI / image processing problem?

Trying to determine BOIs – I wonder if there is any property of a signal such that two separate sub-bands can be identified as being part of one single signal?

Hey Sorin. Thanks for commenting. Yes, when I first started trying to construct an algorithm for automatic BOI estimation, I tried to formulate how I went about identifying BOIs when I stared at a PSD estimate. I came up with nothing profound: a BOI is a arbitrary-shaped bump between two flat spots. All of the subsequent work is trying to ensure that I can identify flat spots! Not as easy as it sounds once we look at real data sets from receivers with imperfections and spectral bands with strange goings-on, not to mention the dependence of the flatness on the parameters of the non-parametric PSD estimator.

One of the papers I referenced in the post tries to apply some image-processing tools to the BOI estimation problem. Not in an AI context though.

The second point you bring up, about trying to group two or more seemingly separate (in frequency) BOIs, is something I’ve tackled as well. The first time was in the context of textbook rectangular-pulse PSK/QAM and textbook MFSK. Each has significant sidelobes. Essentially I use symmetry to try to group bands together. Another context is a signal that is frequency hopping. You’d like to group (associate) all the individual BOIs that you might detect over time and frequency. I have ways of doing this that work well enough for simple FH models, but they are not foolproof, as they rely solely on spectrogram-like analysis (no CSP!).

Finally, I suppose if you can find the BOIs, interrogate their contents using CSP, then you might be able to group disparate BOIs together based on their shared cycle frequencies.

What about signals with more complex spectrum – like CPM/CPFSK that don’t look as nice as sqrtrc QAM/PSK? Or non line of sight signals that have a deep notch in the spectrum. How would you be able to decide if this is one signal or more, based on the shape alone.

I think the spectrum-based approach will fail for a lot of these kinds of cases. You can always come up with a propagation channel and ideal spectrum shape that when combined yields a measurable spectrum that looks exactly like the spectrum of two adjacent-channel signals. So, I admit failure! The multi-resolution spectrum-estimation approach I’ve outlined (and use every day!) cannot, by itself, resolve all ambiguities. Sigh. I suppose, however, you might do better if you can use additional information, such as the patterns of cycle frequencies that you see by (1) processing each “signal” separately, and (2) processing them after forcing them to be joined together. But that is going to cost some CPU cycles.

For signals with multi-modal spectra (like some CPM and traditional FSK signals), if the valleys between the peaks aren’t too shallow, we can still detect the correct full signal band with this kind of PSD-based approach. But once those valleys get small (near zero for a noiseless signal PSD), then my method will start to declare them as separate signals.

In fact, when I apply it to a textbook rectangular-pulse PSK signal, I often get an output consisting of a separate BOI for each of the sidelobes. Thankfully, I do not often encounter such spectra in captured data.

Hi Chad, I’m having a bit of trouble following your notation to try out this approach as compared to one I’m currently using. Do you have any example code (e.g., MATLAB)? Can this iterative algorithm work with a “baseline” PSD at a nominal high resolution and the compute the windowed spectra off that? I can have FPGA firmware accumulate up the PSD but it certainly appears that this type of algorithm requires software. An interesting low-cost addition is to “filter” the PSD with a non-linear denoising algorithm like Total Variance which can remove the fuzz but preserve the larger excursions.

Thanks for the comment Marty, and for visiting the CSP Blog.

I don’t have automatic spectral segmentation code for you, and I’m truly sorry for the difficult notation. Do you think there is an error using the notation or an inconsistency? The telescoping action makes it hard to write the algorithm in a readable way.

No, the algorithm requires access to the data itself to perform the successive PSD estimates as it telescopes.

I’m no hardware or FPGA expert, but it does seem to me this is a logic-heavy and iterative algorithm with an input-data-dependent number of iterations. I’ve certainly not convinced anybody to try to implement it in an FPGA.

Maybe your best option is to just try to take the core idea, of telescoping or multiresolution/iterative spectrum estimation and try to make your own version of the algorithm?

Can you provide a link to what you mean by the Total Variance Denoising algorithm?

Thanks for the clarification Chad, I’m having someone look to build something from your 2007 paper, or use it for inspiration for an approach that starts with an accumulated PSD. An FPGA is pretty efficient in running multiple FFTs and then accumulating the results for a fixed transform size. If I accumulate a couple hundred FFT frames the PSD noise drops to below 0.5 dB so I have a nice clean spectrum to start with.

A good place to look over some denoising algorithms is the work by Dr. Ivan Selesnick at NYU. http://eeweb.poly.edu/iselesni/research/index.html