Fourier Analysis – Derivative of Function Using Discrete Fourier Transform in MATLAB

fourier analysisMATLAB

I'm trying to find the derivative of a function $f(t)$ using the discrete fourier transform. My end goal is to do so numerically, but I suppose if I have the theory down then the rest should follow relatively easily.

From what I've found online, it works such that if $f(t) \rightarrow \hat{f}(\xi)$, then $f'(t) \rightarrow 2\pi i \xi\hat{f}(\xi)$

So if I have a dataset of a periodic signal, I thought that I could approximate its derivative by using a discrete fourier transform, multiplying it by $2 \pi i \xi$ and inverse fourier transforming it. However, it turns out that is is not exactly working out..

What I did was

t = linspace(0,4*pi,4096);
f = sin(t);

fftx = fft(f);
for l = 1:length(fftx)
    dffft(l) = 2*pi*1i*l*fftx(l);
end
df = ifft(dffft);

But the corresponding signal isn't just the cosine of t. Could anyone point out the stupid mistake?

Best Answer

There are a few issues with your code. The first is the use of linspace. It includes both endpoints of the interval, thus both $0$ and $4\pi$ appear. This is not what you want to happen with the discretization for the purpose of Fourier transform. E.g., linspace(0,2*pi,4) returns the vector $0,2\pi/3, 4\pi/3, 2\pi$ but what we actually want is $0,\pi/2, \pi , 3\pi/2$. Solution: introduce the step $dx=2\pi/N$ and create the vector a+[0:N-1]*dx.

Second, the correct version of $2\pi i \xi$ in the discrete setting is not obvious, due to multiple ways to enumerate the terms of Fourier series. In the code below this role is played by vector $k$. I adapted it from Finding Derivatives using Fourier Spectral Methods. 

a = 0;
b = 4*pi;
N = 4096;
dx = (b-a)/N;
t = a + dx*(0:N-1);
f = sin(t);
fftx = fft(f);
k = 2*pi/(b-a)*[0:N/2-1, 0, -N/2+1:-1];
dffft = 1i*k.*fftx;
df = ifft(dffft);
plot(t,df)

Why the strange $k$ vector? Part of it is scaling: the "wider" the function, the more narrow is its Fourier transform (it scales in the opposite direction). This is why we divide by $(b-a)$ there. To understand the weird [0:N/2-1, 0, -N/2+1:-1], think of what the components of fft are: they are the coefficients of the complex polynomial that matches $f$. For a simple example, take t=0:pi/4:(2*pi-pi/4); that is, $N=8$ here. Applying fft to 

f = 3+7*exp(1i*t)+8*exp(3i*t)+9*exp(5i*t)+10*exp(7i*t)

we get 8*[3 7 0 8 0 9 0 10] which is, as expected, the coefficients of exponentials in $f$ (multiplied by $8$, which isn't important). Multiplying the coefficients by 1i*[0:7] gives us the Fourier coefficients of $f'$, and ifft returns $f'$ itself. Sweet. As long as we work with polynomials of $e^{it}$, all this works fine.

Trouble arises when our $f$ has negative powers of $e^{it}$, for example $\sin t = \frac{1}{2}(e^{it}-e^{-it})$. Where does $e^{-it}$ go under the fft? Same place where $e^{7it}$ goes, since $e^{-it} \equiv e^{7it}$ on our grid. But... then our algorithm will multiply it by $7i$, not by $-i$! Most assuredly, $$(\sin t)' \ne \frac{1}{2}(ie^{it}-7ie^{-it})$$ This is the problem that [0:N/2-1, 0, -N/2+1:-1] solves. It simply recognizes that $e^{7it}$ should be treated as $e^{-it}$ for the purpose of taking the derivative: that is, it should be multiplied by $-i$, not by $7i$.

With $N=8$, the multiplying vector is [0 1 2 3 0 -3 -2 -1]. It

  • kills the constant term,
  • multiplies $e^{ikt}$ component by $k$ for $k=1,2,3$,
  • kills the $e^{4it}$ component (What else to do with it? We don't know if it's $e^{4it}$ or $e^{-4it}$. Note that $\sin 4t $ will have zero DFT at this scale).
  • multiplies $e^{ikt} = e^{i(k-8)t}$ component by $k-8$ for $k=5,6,7$.
Related Question