[Numpy-discussion] zoom FFT with numpy?

[Anne Archibald]

On 15/03/07, Ray Schumacher <subscriber100 at rjs.org> wrote:

> The desired band is rather narrow, as the goal is to determine the f of a
> peak that always occurs in a narrow band of about 1kHz around 7kHz
>
>  >2) frequency shift, {low pass}, and downsample

By this I would take it to mean, multiply by a complex sinusoid at 7
kHz, optionally low-pass filter, dispose of every high frequency in
the result by downsampling (boxcar averaging will be an ultra-simple
way to include a crummy low-pass filter) then FFT. (Or some equivalent
purely-real operation, doesn't much matter.) Basically implement a
heterodyne receiver in software.

>  > Much depends on the length of the data run and what the signal to noise is.
>
> data is stable for 2-4k samps, and s/n is low, ~5-10

It seems like this would benefit from some filtering... you could try
numpy.convolve to apply a FIR filter, but scipy.signal has a whole
heap of tools for efficiently designing and applying FIR and IIR
filters.

If you've got a good filter design arrangement, I suppose you could
just band-pass filter the data then downsample and FFT; the aliasing
would drop the frequency you want into the band-pass in a predictable
way. Saves generating and multiplying the  complex sinusoids.

> What I had been doing is a 2048 N full real_FFT with a Hann window, and
> further analyzing the side lobe/bin energy (via linear interp) to try to
> more precisely determine the f within the peak's bin. (How legitimately
> valuable that is I'm not sure... IANAM)

This can be a good idea; if what you're observing really is an
isolated sinusoidal signal, you can get frequency accuracy limited
only by your signal-to-noise ratio. Off the top of my head I don't
remember how you do the magic, but Ransom et al. 2002 ("Fourier
Techniques for Very Long Astrophysical Time Series Analysis",
published in ApJ but more accessible at
http://lanl.arxiv.org/abs/astro-ph/0204349 ; hi Scott!) describe the
process without undue pain. They also discuss the issue of computing
only part of an FFT and claim that it's rarely worth the trouble.

> It's awfully hard to beat the FFTPACK library's full FFT by hand-rolling
> anything else, I'd guess. I was hoping for some all-numpy procedure that I
> could use weave with; that might stand a chance of being faster.

Well, mixing it down then FFTing should be pretty fast. Particularly
if the problem is that you have to do it a zillion times, rather than
that your FFT is huge, so that you can keep the LO array around. And
scipy.signal should help.

Anne M. Archibald