[Matrix-SIG] FFT derivative

[Warren B. Focke]

On Thu, 4 Nov 1999, Jean-Bernard ADDOR wrote:
> 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?

Try this:

--------------------------
import FFT, Numeric, umath

def deriv(series, num=None, axis=-1):
    fspace = FFT.fft(series, num, axis)
    fspace = fderiv(fspace, num, axis)
    return(FFT.inverse_fft(fspace, num, axis))

def fderiv(fspace, num=None, axis=-1):
    if num == None:
        num = fspace.shape[axis]
    n2 = (num - 1) / 2
    iomega = ((Numeric.arange(num) + n2) % num) - n2
    iomega = iomega.astype('D') * 2.0 * Numeric.pi * 1.0j / num
    dims = len(fspace.shape)
    if dims > 1:
        newshape = [1] * dims
        newshape[axis] = num
        iomega.shape = newshape
    dspace = fspace * iomega
    # if the length is even, zero out the value at nyquist
    nyq = num / 2
    if (nyq*2 == num):
        if dims > 1:
            sly = [slice(None)] * dims
            sly[axis] = nyq
        else:
            sly = nyq
        dspace[sly] = 0.0
    return(dspace)
------------------------------

Note that the results will still look "bad" if the sequence you are taking
the derivative of is not continuous at high order across its endpoints.

>     # from the FT theory it should be -1j instead of 1j, why?

Probably because some people put the minus sign on the forward transform,
and others on the reverse, and FFTPACK does the opposite of what you
expected.


Warren Focke