[Matrix-SIG] FFT derivative

[William Park]

On Thu, Nov 04, 1999 at 05:00:43PM -0500, Jean-Bernard ADDOR wrote:
> Hey Numerical people!
> 
> I wrote a python code for derivation in Fourrier space, and it gives very
> bad results. The idea is to multiply the transform of the data by -ik.
> 
> Do you have any hint? Do you know about any code in any language which do
> that, and I could just translate into python? Are any books with detailled
> example of the use of the FFT to compute derivatives?

Do you want h'(t) or H'(f)?  For your reference,
    h'(t)		<==>  j2\pi f H(f)
    -j2\pi t h(t)	<==>  H('f)
so the discrete version can be found by substituting
    f = n/(NT),  t = kT
where T is time sampling interval.

> 
> 
> 	Jean-Bernard
> 
> def FFTderivative(a, order=1, axis=-1):
>     from Numeric import arange
>     from FFT import fft, inverse_fft
>     a = Numeric.asarray(a) # works on arbitrary python sequence
>     l = a.shape[axis]/2
>     m = a.shape[axis] - 2*l
>     if axis != -1: a = Numeric.swapaxes(a, axis, -1)
>     b = fft(concatenate((a, a[...,::-1]), -1)) # JD trick, could be made
> nicer
>     c =
> (1j*concatenate((arange(a.shape[-1]),arange(-a.shape[-1],0))))**order
>     #c = (1j*concatenate((arange(l+m),arange(-l,0))))**order # was before
> fourier space doubling
>     # from the FT theory it should be -1j instead of 1j, why?

The usual engineering convention is to use e^{-j2\pi ft} to go from time
domain to frequency domain.  But, there are people who still use
e^{j2\pi ft} or, even worse, e^{j\omega t}.

>     d = b*c
>     r = inverse_fft(d)[...,:a.shape[-1]]
>     if a.typecode() != 'D' and a.typecode() != 'F': # not complex type
>         r = r.real
>     if axis != -1: r = Numeric.swapaxes(r, axis, -1)
>     return r
> 
> 
> 
> _______________________________________________
> Matrix-SIG maillist  -  Matrix-SIG@python.org
> http://www.python.org/mailman/listinfo/matrix-sig