[Numpy-discussion] Bug or surprising undocumented behaviour in irfft

[Charles R Harris]

On 8/29/07, Charles R Harris <charlesr.harris at gmail.com> wrote:
>
> Anne,
>
> On 8/29/07, Anne Archibald <peridot.faceted at gmail.com> wrote:
> >
> > Hi,
> >
> > numpy's Fourier transforms have the handy feature of being able to
> > upsample and downsample signals; for example the documentation cites
> > irfft(rfft(A),16*len(A)) as a way to get a Fourier interpolation of A.
> > However, there is a peculiarity with the way numpy handles the
> > highest-frequency coefficient.
>
>
> <snip>
>
> The upshot is, if I correctly understand what is going on, that the
> > last coefficient needs to be treated somewhat differently from the
> > others; when one pads with zeros in order to upsample the signal, one
> > should multiply the last coefficient by   0.5. Should this be done in
> > numpy's upsampling code? Should it at least be documented?
>
>
> What is going on is that the coefficient at the Nyquist  frequency appears
> once in the unextended array, but twice when the array is extended with
> zeros because of the Hermitean symmetry. That should probably be fixed in
> the upsampling code.
>

The inverse irfft also scales by dividing by the new transform size instead
of the original size, so the result needs to be scaled up for the
interpolation to work. It is easy to go wrong with fft's when the correct
sampling/frequency scales aren't carried with the data. I always do that
myself so that the results are independent of transform size/interpolation
and expressed in some standard units.


In [9]: a = array([1, 0, 0, 0], dtype=double)

In [10]: b = rfft(a)

In [11]: b[2] *= .5

In [12]: irfft(b,8)
Out[12]:
array([ 0.5      ,  0.3017767,  0.       , -0.0517767,  0.       ,
       -0.0517767,  0.       ,  0.3017767])

In [13]: 2*irfft(b,8)
Out[13]:
array([ 1.        ,  0.60355339,  0.        , -0.10355339,  0.        ,
       -0.10355339,  0.        ,  0.60355339])

I don't know where that should be fixed.

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