spectral derivative using fft()

[Ben Abbott]

On Dec 17, 2007, at 5:39 PM, Gastón Araguás wrote:

> hi,
>
> I want to implement an spectral derivative using fft() and ifft().
> Theoretically it can be done scaling the fft coefficients by the
> factor  (2 pi i w)^n, where n is the derivate order, and applying the
> ifft() to them.
> It is, the derivative of a function f_j(x_j) can be obtained by
> derivation of it's spectral representation, so if we have
>
> f_j = sum_{w=-M} ^{M} e^{2 pi i w j} ftilde(w)
>
> with j = 0, ..., 2M, then the f'_j can be obtained by
>
> f'_j = ( 2 pi i w )*sum_{w=-M} ^{M} e^{2 pi i w j} ftilde(w)
>
> in terms of the fft() and ifft() functions it should be:
> f'_j = ifft( ( 2 pi i w ) * fft( f_j ) )
> with w = -M, ..., M
>
> but it doesn't gives the correct result. I suspect that the problem is
> in the order of the returned fft() coefficients.
> I found that the coef. ftilde(w) are stored in this way:
> [ftilde(0) ftilde(1) ...   ftilde(M-1/2) ftilde(M+1/2) ftilde(-M+1/2)
> ... ftilde(-1)] = fft(f_j)
> or, defining N = 2M+1 it holds
> [ftilde(0) ftilde(1) ...   ftilde(N/2-1) ftilde(N/2) ftilde(-N/2+1)
> ... ftilde(-1)] = fft(f_j)
> is that ok?
>
> Here is my test code (it works ok with this simple function y =
> cos(2*pi*x), but not for example with y = x + cos(2*pi*x) ) :
> N = 2^7;
> h = 1/N;
> x = [0:N-1]*h;
> y = cos(2*pi*x);
> figure; plot (real(y));
> Y = fft( y , N );
> M = (N-1)/2;
> W = (2*pi*i)*[0:M-1/2 0 -M+1/2:-1];
> Y2 = W.*Y;
> y2 = ifft( Y2 );
> figure; plot (real(y2));
> pause;
>
> i hope someone can help, thanks in advanced
> Gastón

At 1st glance, I noticed a problem. You used the following to  
represent the derivative of "f" (I've modified it slightly).

	fp = ifft ( ( 2i * pi  * m ) * fft( f ) );

which is correct for continuous signals with infinite limits (Fourier  
Transform). However, for FFT's the signals are discrete & periodic  
(Discrete Fourier Transforms). So the derivative must be periodic as  
well. I haven't gone through all the details, but the result should  
look something like

	fp = 1i *  ifft ( sin ( 2 * pi * m ) .* fft ( f ) );

I'm likely off by some scale factor, but you get the idea, yes?

Ben