[SciPy-User] fft not giving Hermitian output for real input
Aronne Merrelli
aronne.merrelli@gmail....
Sat Feb 25 11:24:11 CST 2012
On Thu, Feb 23, 2012 at 8:15 PM, Jeff Alstott <jeffalstott@gmail.com> wrote:
> I have a particular data set, d0, which is real:
>
> In [72]: d0
> Out[72]:
> array([ 0.00907105, 0.0372916 , 0.01402867, ..., -0.04779497,
> -0.07054817, -0.0436582 ])
>
> It doesn't produce Hermitian output from numpy.fft.fft.
>
>
> In [87]: y0 = fft.fft(d0)
> In [88]: y0
> Out[88]:
> array([ 7.77156117e-14 +0.00000000e+00j,
> 6.89226454e-13 -1.56319402e-13j,
> 1.72140080e-13 +0.00000000e+00j, ...,
> -7.95807864e-13 -1.13686838e-13j,
> -3.41060513e-13 +1.13686838e-13j, 1.25055521e-12
> -3.41060513e-13j])
> In [89]: y0[1].conj()==y0[-1]
> Out[89]: False
> In [90]: y0[1].conj()==y0[-2]
> Out[90]: False
>
>
I might be missing something, but this seems like floating point rounding
error to me. I don't know what your d0 really looks like, but with some
random numbers,
In [1]: d0 = np.random.normal(0,1,1024)
In [2]: y0 = np.fft.fft(d0)
In [3]: np.allclose( y0[1:].conj(), y0[-1:0:-1] )
Out[3]: True
In [4]: np.max( (y0[1:].conj() - y0[-1:0:-1]).imag )
Out[4]: 4.3520742565306136e-13
In [5]: np.max( (y0[1:].conj() - y0[-1:0:-1]).real )
Out[5]: 2.957634137601417e-13
I tried this in MATLAB and the two halves are exactly equal. I would assume
this means that the MATLAB FFT implementation does something extra to
eliminate the rounding error for the special case of real input, whereas
the NumPy FFT does not take that step. I'm pretty sure the NumPy FFT is
just calling FFTPACK, so you'd need to check the implementation details
there.
Cheers,
Aronne
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.scipy.org/pipermail/scipy-user/attachments/20120225/40ffee20/attachment.html
More information about the SciPy-User
mailing list