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.
