Algorithm for automatic channel detection

Question

I'm currently working on a spare-time project to perform automatic modulation classification (AMC) on radio signals (more precisely, I'm interested in L-band satellite channels), using SDR. I would like it to be able to discover channels in the power spectrum, along with their frequencies and bandwidths, so I can direct the application of AMC on the output of this step.

My first (naive) approach was rather simple: after every N samples (say 1024) apply a window function, perform the FFT on the last N, compute an estimation of the power spectrum and apply some exponential smoothing to reduce the noise. Then, in the smoothed power spectrum, find the maximum and minimum signal levels, calculate some threshold value based on a weighted mean of both levels and use this threshold to determine which frequency bins belong to a channel.

This works well in my unit tests (one QPSK channel + gaussian noise). However, in real-life scenarios I either get a few channels or a lot of false-positives. Of course I can fix this by fine-tuning the weights in the threshold calculation, but then it wouldn't be automatic anymore.

I've been doing some research on Google but maybe I'm not using the right search keywords, or there is no real interest on this subject (which would be strange, as frequency scanners must perform this task somehow).

How could I find the appropriate values for the mean weights? Maybe there is a better approach than calculating a threshold for the whole spectrum?

Answer 1

Your smoothing approach seems a bit counter-productive: Why would noise in a sufficiently long DFT form something like a "sharp" shape? You're most likely suppressing narrowband carriers that way.

Then: There's really a lot of signal detectors, many simply based on energy detection in spectrum (just like your approach).

However, the better an estimator has to be, the more info on the signal you're looking for you'll need to have. Are you looking for 20 kHz wide narrowband channels, or for dozens of Megahertzes of high-rate QPSK satellite downlink? Do you know anything about the channel/pulse shape? Autocorrelation?

This is a bit of a wide field, so I'd propose you look at something that already works:

gr-inspector of Sebastian Müller is a proven and pretty awesome automatic channel detector, and can, for some types of transmissions, also infer modulation parameters.

See a demo video (dubstep warning!); what you seem to be sketching looks something like this in that:

But: that's just one of the things it can do.

More importantly: Energy-band detection based on a DFT is really just one of many things you can do to detect signals. It wouldn't, for example, detect spread-spectrum transmissions, like GPS, for example. For that, you'd need knowledge of the spreading technique or at least a large autocorrelation-based detector.

Answer 2

After considering Marcus Müller's suggestion, I developed an algorithm that still relies on a global energy threshold but also on an estimation of the noise floor. It can be improved in many ways, but as its simplest realization already provides acceptable results with real-world captures (again, L-band captures at 250 ksps) I believe it could be a good ground for other people to start with.

The central idea of the algorithm is to calculate the energy threshold based on a continuous estimation of the noise power, updating it with every update of the spectrogram. This estimation keeps track of the maximum and minimum levels attained by each FFT bin after during all FFT runs, using them to estimate the PSD in that bin and discard outliers (i.e. channels, spurs...). The algorithm is to be executed every fixed number of samples. More formally:

Parameters of the algorithm:

int N        /* spectrogram length (FFT bins) */
complex x[N] /* last N samples */
float alpha  /* spectrogram smoothing factor between 0 (infinite smoothing,
                spectrogram will never be updated) and 1 (no smoothing at all) */
float beta   /* level decay factor between 0 (levels will never decay) and 1
                (levels will be equal to the spectrogram). */
float gamma  /* smooth factor applied to updates of the current noise estimation
                between 0 (no updates allowed) and 1 /* no smoothing */
float SNR    /* detection threshold above noise, in dB */

Recommended values for alpha, beta and gamma are 1e-2, 1e-3 and .5 respectively.

Variables:

complex dft[N]  /* Fourier transform of the last N samples */
float spect[N]  /* smoothed spectrogram */
float spmin[N]  /* lower levels for spectrogram bins */
float spmax[N]  /* upper levels for spectrogram bins */
int runs        /* FFT run counter (initially zero) */
float N0        /* Current noise level estimation */
float new_N0    /* New noise level estimation */
float min_pwr   /* Minimum power density found */
int min_pwr_bin /* FFT bin where min_pwr is */
int valid       /* Number of valid bins for noise estimation */
int min_runs    /* Minimum number of runs required to detect channels */

Algorithm:

min_runs = max(2. / alpha, 1. / beta)
dft = FFT(x);
++runs;
if (runs == 1) then /* First FFT run */
    spect = dft * conj(dft) /* |FFT(x)|^2 */
    spmin = spect /* Copy spect to spmin */
    spmax = spect /* Copy spect to spmax */
    N0    = min(spect); /* First noise estimation */
else
    /* Smooth spectrogram w.r.t the previous run */
    spect += alpha * (dft * conj(dft) - spect)

    /* Update levels. This has to be performed element-wise  */
    new_N0 = 0
    valid = 0
    min_pwr = INFINITY
    min_pwr_bin = -1
    for (int i = 0; i < N; ++i)
        /* Update current lower levels or raise them */
        if (spect[i] < spmin[i]) then
            spmin[i] = spect[i]
        else
            spmin[i] += beta * (spect[i] - spmin[i]);
        end

        /* Update current upper levels or decrease them */
        if (spect[i] > spmax[i]) then
            spmax[i] = spect[i]
        else
            spmax[i] += beta * (spect[i] - spmax[i]);
        end

        if (runs > min_runs) then
            /* Use previous N0 estimation to detect outliers */
            if (spmin[i] < N0 or N0 < spmax[i]) then
                new_N0 += spect[i]
                ++valid
            end

            /* Update current minimum power */
            if (spect[i] < min_pwr) then
                min_pwr = spect[i]
                min_pwr_bin = i
            end
        end
    end

    /*
     * Check whether levels have stabilized and update noise
     * estimation accordingly
     */
    if (runs > min_runs) then
        /*
         * This is a key step: if the number of valid bins is
         * 0 this means that our previous estimation was
         * absolutely wrong. We reset it with a cruder estimation
         * based on where the minimum value of the current
         * spectrogram was found
         */

        if (valid == 0) then
            N0 = .5 * (spmin[min_pwr_bin] + spmax[min_pwr_bin])
        else
            N0 += gamma * (new_N0 / valid - N0)
        end

        /*
         * Detect channels based on this threshold (trivial,
         * not detailed)
         */

        detect_channels(spect, 10^(SNR / 10) * N0)
    end
end

Even though this algorithm makes the strong assumption that the noise floor is flat (which is false in most cases as in real-world radios, tuner output passes through a low-pass filter whose response is not flat), it works even if this condition doesn't hold. These are some of the algorithm results for different values of alpha, N = 4096 and SNR = 3 dB. Noise estimation is marked in yellow, channel threshold in green, upper levels in red, spectrogram in white and lower levels in cyan. I also provide an evolution of the N0 estimation after every FFT run:

Results for alpha = 1e-1:

Results for alpha = 1e-2. Note how the number of valid bins has been reduced as the spectrogram got clearer:

Results for alpha = 1e-3. In this case, the levels are so tight and the noise floor so obviously non-flat that there are novalid bins from one FFT run to another. In this case we fall back to the crude estimation of looking for the bin with the lowest power density:

The min_runs calculation is critical. To prevent the noise level to updrift (this is, to follow a channel and not the noise floor) we must wait at least 2. / alpha FFT runs before trusting the signal levels. This value was found experimentally: in my previous implementations, I was intuitively using 1. / alpha which failed miserably for alpha = 1e-3:

I haven't tested this yet on other scenarios (like burst transmissions) where this algorithm may not perform as well as with continuous channels because of the persistence of min/max levels, and it may fail to detect burst transmissions as outliers. However, given the kind of channels I'm working with, this is not a priority for now.