[Numpy-discussion] Bug or surprising undocumented behaviour in irfft
Wed Aug 29 21:24:50 CDT 2007
On 29/08/2007, Charles R Harris <email@example.com> wrote:
> > Is this also appropriate for the other FFTs? (inverse real, complex,
> > hermitian, what have you) I have written a quick hack (attached) that
> > should do just that rescaling, but I don't know that it's a good idea,
> > as implemented. Really, for a complex IFFT it's extremely peculiar to
> > add the padding where we do (between frequency -1 and frequency zero);
> > it would make more sense to pad at the high frequencies (which are in
> > the middle of the array). Forward FFTs, though, can reasonably be
> > padded at the end, and it doesn't make much sense to rescale the last
> > data point.
> It all depends on the data and what you intend. Much of my experience is
> with Michaelson interferometers and in that case the interferogram is
> essentially an autocorrelation, so it is desirable to keep its center at
> sample zero and let the left side wrap around, so ideally you fill in the
> middle as you suggest. You can also pad at the end if you don't put the
> center at zero, but then you need to phase shift the spectrum in a way that
> corresponds to rotating the center to index zero and padding in the middle.
> I expect you would want to do the same thing for complex transforms if they
> are of real data and do the nyquist divided by two thingy. If the high
> frequencies in a complex transform are actually high frequencies and not
> aliases of negative frequencies, then you will want to just append zeros.
> That case also occurs, I have designed decimating complex filters that
> produce output like that, they were like single sideband in the radi o
So is it a fair summary to say that for irfft, it is fairly clear that
one should adjust the Nyquist coefficient, but for the other varieties
of FFT, the padding done by numpy is just one of many possible
Should numpy be modified so that irfft adjusts the Nyquist
coefficient? Should this happen only for irfft?
> I usually multiply the forward transform by the sample interval, in secs or
> cm, and the unscaled inverse transform by the frequency sample interval, in
> Hz or cm^-1. That treats both the forward and inverse fft like
> approximations to the integral transforms and makes the units those of
> spectral density. If you think trapezoidal rule, then you will also see
> factors of .5 at the ends, but that is a sort of apodization that is
> consistent with how Fourier series converge at discontinuities. In the
> normal case where no interpolation is done the product of the sample
> intervals is 1/N, so it reduces to the usual convention. Note that in your
> example the sampling interval decreases when you do the interpolation, so if
> you did another forward transform it would be scaled down to account for the
> extra points in the data.
That's a convenient normalization.
Do you know if there's a current package to associate units with numpy
arrays? For my purposes it would usually be sufficient to have arrays
of quantities with uniform units. Conversions need only be
multiplicative (I don't care about Celsius-to-Fahrenheit style
conversions) and need not even be automatic, though of course that
would be convenient. Right now I use Frink for that sort of thing, but
it would have saved me from making a number of minor mistakes in
several pieces of python code I've written.
More information about the Numpy-discussion