7
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ Hm, the 1/N I see, that was a bit stupid indeed. Why should I subtract 1 though? Your fix still doesn't completely work though. The sign of the derivative is incorrect (I'm seeing minus a cosine), and the amplitude also seems off. Also, should I take the absolute value of the ifft? As the result of ifft is still complex. $\endgroup$ – user3183724 Apr 5 '14 at 18:06
  • $\begingroup$ Right, right, fair enough. That's a bit silly of me then. $\endgroup$ – user3183724 Apr 5 '14 at 18:11
10
$\begingroup$

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$.
$\endgroup$
  • $\begingroup$ I almost understand, thank you for the effort. The only thing that is a little vague to me is the k = 2*pi/(b-a)*[0:N/2-1, 0, -N/2+1:-1]; part. I don´t really follow where that comes from, nor do I understand it in the link you´ve provided. It´s 2pi times the length of the interval, then multiplied by an array that is not obvious to me. $\endgroup$ – user3183724 Apr 5 '14 at 21:26
  • $\begingroup$ @user3183724 I expanded the answer. $\endgroup$ – user127096 Apr 5 '14 at 22:34
  • 1
    $\begingroup$ Thank you. I still find it rather complex, but I understand now. I should have been able to come up with it on my own however, as it was part of an assignment, but I really do not see how. Ah well, that's something for me to worry about! $\endgroup$ – user3183724 Apr 6 '14 at 13:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.