[SciPy-user] Fourier-sine transform with SciPy

Anne Archibald peridot.faceted@gmail....
Fri Jun 6 04:42:29 CDT 2008


2008/6/6 Lubos Vrbka <lists@vrbka.net>:

>> I don't think we have a function that computes discrete sine or cosine
>> transforms, but if f is real, you can get the sine transform you wrote
>> as the imaginary part of a complex Fourier transform.
> by this you mean using scipy.fftpack.rfft()? i must admit i am a bit
> lost in the various possibilities that are available.

Yes, there are rather a lot.

A complex fft - numpy.fft.fft - computes:

F_m = A \sum_{k=0}^{N-1} f_k (cos(2\pi m k/N)+i*sin(2\pi m k/N))

and an ifft computes:

f_k = B \sum_{m=0}^{N-1} F_m (cos(2\pi m k/N)-i*sin(2\pi m k/N))

So if you want

\sum_{k=0}^{N-1} f_k sin(2\pi m k/N)

all you have to do is take a complex FFT of your (real) data, then
take the imaginary part. This is somewhat wasteful, which is why FFTW
implements the sine and cosine transforms, but it's only a factor of
two or four (i.e. still far better than coding it by hand).

To save a little time, you can use a "real" FFT, which computes F_m as
above, but (1) takes advantage of the fact that f_k is known to be
real and (2) doesn't bother computing F_m for m>N/2 (or so), since
when f_k is real, F_m = conjugate(F_{N-m}) IIRC. So you may want to do
this and take the imaginary part.

Note that the sine transform only works on odd functions, so either
your input data should have f_k = f_{N-k} or you are representing only
half the values, i.e., N is twice the number of data points you have.
This too is somewhere a proper sine transform would help, if we had
it.

> scipy.fftpack.*fft* functions are clear (is this the netlib implementation?)
>
> scipy.ifft?? tells me (among others)
> File:             /usr/lib/python2.4/site-packages/numpy/fft/fftpack.py
>
> or should i use numpy.fft? what's the difference in fftpack in scipy and
> numpy?
>
> is there any place where these things are documented, i.e. what code
> does actually lie behind these functions?
>
> and just out of curiosity - from what i've seen during my search for the
>  information on google, fftw was (at least in the past) used in scipy.
> the scipy package on my debian box depends on fftw2, but the names refer
> to fftpack, that shouldn't be fftw but netlib implementation instead (if
> i am not wrong). could anybody clarify this for me, please? sorry if
> this was discussed before, but i wasn't able to find the relevant thread.

There are several reasons this is complicated. numpy's fft is designed
to use one of several external libraries depending on how numpy is
compiled. FFTPACK is always available, but is not especially fast.
FFTW is distributed under the GPL, so numpy cannot depend on it
without also being distributed under the GPL. But if it's an optional
plug-in, numpy can make use if it where it's available, since it's
much faster than FFTPACK (in spite of the fact that we use its API in
a not-very-efficient manner). numpy also supports several other
libraries, most notably Intel's proprietary and hardware-specific but
very fast implementation. At a python level, the user need not know
which backend is being used, and FFTPACK is mentioned a few times.

Anne


More information about the SciPy-user mailing list