help-octave
[Top][All Lists]

## Re: spectral derivative using fft()

 From: Ben Abbott Subject: Re: spectral derivative using fft() Date: Mon, 17 Dec 2007 19:01:35 -0500

```
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

```