Real interbinning

(This is another of my highly-technical "note to self" signal processing posts. I'll put up something less arcane soon.)


The Fourier transform is great for finding periodic signals: you take an FFT and a periodic signal looks like a peak in the output. Well, in an ideal world, that is; you only really get such a neat and tidy peak if the periodic signal is exactly at a Fourier frequency, which happens when it makes exactly an integral number of turns over the course of the data set. If the signal is somewhere between two Fourier frequencies, the power is spread over several Fourier output values. While it's possible to interpolate a very accurate value based on the ~32 nearest values, this can be expensive, and there's a shortcut called "interbinning"; it doesn't reconstruct the phase, but you just take a scaled average of the two neighbouring bins and get a decent approximation to the value at the midpoint between two independent Fourier frequencies.

My problem, for this post, is that the theory is all nicely worked out when going from the time domain to the frequency domain, but I want to do something analogous while going from the frequency domain back to the time domain, if that's possible. (I haven't done a literature search, or even looked very carefully at the frequency-domain interbinning papers; I thought this would be a good exercise for me.)

Just to be specific about the process I'm trying to analyze: I have the Fourier series of some pulse, assumed a delta function, and I truncate it and do an inverse real FFT, obtaining n points. Now I want to efficiently estimate the amplitude of the original pulse, in the presence of noise (though not much noise, since I need a highly significant detection to be worth following up).

The signal, then, has the profile given by a shifted Dirichlet kernel, sin(pi (n+1) x)/sin(pi x) (note the different definition of n I use from what's in the link). So when the peak is midway between two sampling points, the amplitude of the two neighboring points is a fraction 1/(n tan(pi/2n)) of the height of the main peak (note that there's a little fiddle here because, if we're doing a real inverse FFT, only the real part of the top Fourier coefficient is used; this value takes that into account). This is a gradually increasing function converging to 2/pi, or about 0.637. So if we average two adjacent values and divide by this value, we'll get the height the peak would have had had it been exactly on a sampling point.

Of course, we have noise on every sampling point, and scaling up the average like this will amplify the noise. But the averaging reduced it by a factor root two, so the loss is not as bad as it seems at first. In fact, as the plot up above shows, for all but the smallest sizes the efficiency is always about 2sqrt(2)/pi or 90% (this lowered efficiency essentially comes from using only the power in two samples rather than all of them). A ten percent loss of efficiency is not good, certainly, but we can to some extent compensate for it: simply lower our threshold for interbinned candidates by ten percent, then at the end discard any that don't improve when "polished". Since this approach is necessary anyway for the frequency sampling and interpolation, and since not interbinning roughly doubles our runtime, this is probably reasonable as long as we don't get utterly deluged with candidates as a result.

That said, the ten percent loss is for signals exactly halfway between sampling points. What about other offsets? The plot on the left shows that the worst-case loss of sensitivity (in the case of large n) is about 15%. This is still not great, but as above, I think it is probably manageable. And for a factor of two in runtime, I think it's probably worth it.

Upon reflection, of course I don't have to guess whether it will be worth it: since I know my thresholds are Gaussian, I can compute how many excess false positives I get by reducing the threshold in this way. The answer will depend on the threshold I set, of course. But if I'm running down around one false positive per FFT, that's something like one in 10⁷, and at those levels, there are some hundred bogus false positives due to the lowered threshold for every "real" false positive. That's annoying but not necessarily a problem. With those numbers, though, it would be nice to short-circuit the "polishing" step. If we split the polishing up into steps: inverse real FFT with lots of bins, then search in period with proper interpolation, we can short-circuit at each step if the power doesn't improve enough to pass the real threshold. As long as it's only some hundred candidates, 1024-point FFTs are a very cheap price to pay. Realistically, beams with zillions of RFI candidates are going to be more of a problem, since there's a risk they'll completely bog down the processing, polishing thousands of candidates from a beam that's too contaminated to be worth using. Short-circuiting the polishing can help with these too: you just declare by fiat that not more than, say, a thousand candidates will be reported from any FFT, so once you hit a thousand candidates, you start raising the threshold.