[Numpy-discussion] Bug or surprising undocumented behaviour in irfft
Anne Archibald
peridot.faceted at gmail.com
Wed Aug 29 18:31:55 EDT 2007
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.
First of all, the normalization:
In [65]: rfft(cos(2*pi*arange(8)/8.))
Out[65]:
array([ -3.44505240e-16 +0.00000000e+00j,
4.00000000e+00 -1.34392280e-15j,
1.22460635e-16 -0.00000000e+00j,
-1.16443313e-16 -8.54080261e-16j, 9.95839695e-17 +0.00000000e+00j])
In [66]: rfft(cos(2*4*pi*arange(8)/8.))
Out[66]: array([ 0.+0.j, 0.+0.j, 0.-0.j, 0.+0.j, 8.+0.j])
So a cosine signal gives 0.5*N if its frequency F is 0<F<=N/2, but N
if its frequency is N/2+1 (or zero). This is fine; it's the way
Fourier transforms work.
Now, suppose we take a signal whose Fourier transform is [0,0,0,0,1]:
In [67]: n=8; irfft([0,0,0,0,1],n)[0]*n
Out[67]: 1.0
In [68]: n=16; irfft([0,0,0,0,1],n)[0]*n
Out[68]: 2.0
In [69]: n=32; irfft([0,0,0,0,1],n)[0]*n
Out[69]: 2.0
In [70]: n=64; irfft([0,0,0,0,1],n)[0]*n
Out[70]: 2.0
The value at zero - a non-interpolated point - changes when you interpolate!
Similarly, if we are reducing the number of harmonics:
In [71]: n=8; irfft([0,0,1,0,0],n)[0]*n
Out[71]: 2.0
In [72]: n=4; irfft([0,0,1,0,0],n)[0]*n
Out[72]: 1.0
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?
Thanks,
Anne M. Archibald
More information about the NumPy-Discussion
mailing list