[Numpy-discussion] zoom FFT with numpy?

[Charles R Harris]

On 3/15/07, Warren Focke <focke at slac.stanford.edu> wrote:
>
>
>
> On Wed, 14 Mar 2007, Charles R Harris wrote:
>
> > On 3/14/07, Ray Schumacher <subscriber100 at rjs.org> wrote:
> > >
> > > 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)
> >
> >
> > That's usually fine. You might want to zero fill to get more samples
> through
> > the band. It would help if you gave the sample frequency in Hz too.
> Anyway,
> > unless time is important, I would just zero fill by a factor of 4-8 and
> > transform. You can get the same effect with a chirp-z transform, but
> again
> > this is depends on how much programming work you want to do. If you just
> > have a couple of lines in the band that you want to locate you could
> also
> > try maximum entropy methods.
>
> Or do
>
> 1) explicit multipy by the transform matrix but only use a few rows.
>
> but evaluate the transform matrix on a frequency grid finer than the
> 1/(2048*tsamp) that is actually independent.  Then fit sin(f)/f to the
> power spectrum.


You can actually zero fill by a factor of two, then build optimum least
squares interpolators for bandlimited signals using a reasonable number of
sample around each frequency interval. The result can be fitted to single
precision accuracy with a 9'th degree polynomial and an ordinary zero solver
used on the derivative. Works well, I did this some 20 years ago as part of
a package for fourier spectroscopy, but it is probably more work than
warranted for the present case.

Either of these should give better interpolation than
> linear.  I've seen this done (and pass peer review) to determine pulsar
> frequencies.  I also remain unconvinced that interpolation provides a
> better result, but that can be determined by analyzing fake data with a
> known frequency.
>
> If you're trying to determine the significance of the result, the fit
> should somehow take into account the fact that the interpolated data
> points are not real degrees of freedom.  But your power estimates are
> already not independent since you've applied a Hann window.  Probably
> should also fit to the line response of a Hann window rather than
> sin(f)/f.


Sinc interpolation will work fine for the windowed spectrum as it contains
the same range of frequencies as the original. Where you can gain something
is explicitly interpolating the unwindowed spectrum with the Hann, or
stronger, window. Because the window functions fall off much faster than the
sinc you don't need to use so many points in the convolution.

Chuck
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20070315/edc721f1/attachment.html>