469,926 Members | 2,517 Online
Bytes | Developer Community
New Post

Home Posts Topics Members FAQ

Post your question to a community of 469,926 developers. It's quick & easy.

? on scipy.fftpack

I've been debugging a simulation I wrote a while ago, and it seems that
the problem is in the fft module itself. I'm trying a simple test by
just feeding the function a basic real gaussian. Obviously, I should
get back the same real gaussian, but what I get is not even close. Can
anyone help me? Thanks in advance.

<code>
from scipy import *
import pylab

fft, ifft = fftpack.rfft, fftpack.irfft

def main():
sigma = 50.
x = arange(-100, 100, 1, 'f')
g = exp(-x**2 / (2 * sigma)) / sqrt(2 * pi)
testFFT(x, g)

def testFFT(x, f):
pylab.plot(x, f)
pylab.title('Initial gaussian')
pylab.show()
pylab.clf()

f = fft(f)
pylab.plot(x, f)
pylab.title('After first FFT')
pylab.show()
pylab.clf()

f = fft(f)
pylab.plot(x, f)
pylab.title('After first iFFt')
pylab.show()
pylab.clf()

if __name__ == '__main__': main()
</code>

Jun 28 '06 #1
3 3194
Mike Duffy wrote:
I've been debugging a simulation I wrote a while ago, and it seems that
the problem is in the fft module itself. I'm trying a simple test by
just feeding the function a basic real gaussian. Obviously, I should
get back the same real gaussian, but what I get is not even close. Can
anyone help me? Thanks in advance.


You will probably want to ask scipy questions on scipy-user. There aren't many
scipy people here.

http://www.scipy.org/Mailing_Lists

I haven't run your code, yet, but one of the things you are running into is the
FFT packing convention for FFTs on real functions. Please read the docstring:

Type: function
Base Class: <type 'function'>
String Form: <function rfft at 0x204fbf0>
Namespace: Interactive
File:
/Library/Frameworks/Python.framework/Versions/2.4/lib/python2.4/site-packages/scipy-0.5.0.1980-py2.4-mac
osx-10.4-ppc.egg/scipy/fftpack/basic.py
Definition: rfft(x, n=None, axis=-1, overwrite_x=0)
Docstring:
rfft(x, n=None, axis=-1, overwrite_x=0) -> y

Return discrete Fourier transform of real sequence x.

The returned real arrays contains
[y(0),Re(y(1)),Im(y(1)),...,Re(y(n/2))] if n is even
[y(0),Re(y(1)),Im(y(1)),...,Re(y(n/2)),Im(y(n/2))] if n is odd
where
y(j) = sum[k=0..n-1] x[k] * exp(-sqrt(-1)*j*k* 2*pi/n)
j = 0..n-1
Note that y(-j) = y(n-j).

Optional input:
n
Defines the length of the Fourier transform. If n is not
specified then n=x.shape[axis] is set. If n<x.shape[axis],
x is truncated. If n>x.shape[axis], x is zero-padded.
axis
The transform is applied along the given axis of the input
array (or the newly constructed array if n argument was used).
overwrite_x
If set to true, the contents of x can be destroyed.

Notes:
y == rfft(irfft(y)) within numerical accuracy.

--
Robert Kern

"I have come to believe that the whole world is an enigma, a harmless enigma
that is made terrible by our own mad attempt to interpret it as though it had
an underlying truth."
-- Umberto Eco

Jun 28 '06 #2

Robert Kern wrote:
You will probably want to ask scipy questions on scipy-user. There aren't many
scipy people here.

http://www.scipy.org/Mailing_Lists

I haven't run your code, yet, but one of the things you are running into is the
FFT packing convention for FFTs on real functions. Please read the docstring:


Ok, thanks a lot. I was unaware of that mailing list, I will certainly
go there next. I have read the documentation, but I'm not sure what
packing convention you are referring to.

Jun 28 '06 #3
Oops, sorry. I see what you mean. I was reading the docs for the
regular (complex) fft function, since that was what I was initially
having the bug with. I was just using the rfft to simplify things, but
I guess that was ironic. Anyway, I appreciate your help and don't mean
to burden you. Again, I will suubmit this to the scipy list.

Jun 28 '06 #4

This discussion thread is closed

Replies have been disabled for this discussion.

Similar topics

3 posts views Thread by avik | last post: by
3 posts views Thread by hawkesed | last post: by
1 post views Thread by tkpmep | last post: by
reply views Thread by Julien Fiore | last post: by
2 posts views Thread by Aage Andersen | last post: by
1 post views Thread by tkpmep | last post: by
2 posts views Thread by robert | last post: by
2 posts views Thread by Peter Maas | last post: by
2 posts views Thread by Frank Moyles | last post: by
By using this site, you agree to our Privacy Policy and Terms of Use.