CODEWITHSUNDEEP
matlabsignal-processingfftdft
Mikhail Genkin
Mikhail Genkin

Reputation: 3460

Discrepancy with differentiation in Fourier space

I am solving differential equations in Fourier space using Matlab. However, I got a problem: after differentiating my real signal I get complex answer (which is incorrect).

Consider an example with diffirentiation over x (multiplication by ik in fourier space):

 a=rand(6,1).';
 fr=fftshift(-3:1:2);
 ifft(1i*fr.*fft(a))

The output is complex. I figured out why it happens: our frequency spectrum is -3,-2,-1,0,1,2. So, there is no pair for the highest frequency (we have -3, but don't have 3). I am wondering how to fix it.

If we think about it, there is technically non-zero contribution from the highest frequency. If the fourier amplitude on the frequency -3 is c0, that means, that in reality we have amplitude c0/2 on frequencies -3 and 3, so that after differentiation we get:

(c0/2)*i*(-k)*exp(-ikx)+(c0/2)*i*(k)*exp(ikx)=kc0*sin(kx)

I am curious, how to implement correct differentiation. My problem is 2D, so that I use fft2 and ifft2. But problem is of the same origin.

Thanks

Upvotes: 1

Views: 208

Answers (3)

Luis Mendo
Luis Mendo

Reputation: 112749

You need to take three things into account:

  • Differentiating in time corresponds to multiplying the DFT by 1-exp(-1j*2*pi/N*fr), where N is the signal period and fr = 0:N-1 are the frequency samples. This stems from the time-shift property of the DFT (see for example here).
  • This differentiation should be circular (see above link), because the DFT inherently assumes that the time signal is periodic. So the first differentiated sample in the time domain will be a(1)-a(N), the second will be a(2)-a(1), etc.
  • You may get a very small imaginary part due to floating-point accuracy.

So, the code should be:

a = rand(6,1).';
N = numel(a);
fr = 0:N-1;
a_diff_fr = ifft((1-exp(-1j*2*pi/N*fr)).*fft(a));

Check:

>> a_diff_fr % imag part should be small
a_diff_fr =
  Columns 1 through 5
  -0.5490 - 0.0000i   0.3169 - 0.0000i  -0.5662 + 0.0000i   0.6851 + 0.0000i  -0.5155 - 0.0000i
  Column 6
   0.6287 + 0.0000i

>> real(a_diff_fr) % real part only
ans =
   -0.5490    0.3169   -0.5662    0.6851   -0.5155    0.6287

>> a([1 2:N])-a([N 1:N-1]) % circular differentiation
ans =
   -0.5490    0.3169   -0.5662    0.6851   -0.5155    0.6287

Upvotes: 2

Mikhail Genkin
Mikhail Genkin

Reputation: 3460

Luis Mendo suggested correct solution to my problem. I also figured out how to fix my approach: the think is, that sinusoid component is always zero on the high frequency, only cosinusoidal harmonic is visible:

signal=sin(2.*(linspace(0,2*pi*(1-1/4),4)));
q=fftshift(fft(signal))./4

Here q is zero. But if do it with cos(2x) signal:

signal=cos(2.*(linspace(0,2*pi*(1-1/4),4)));
q=fftshift(fft(signal))./4

Here q(1)=1. So, in my approach I have to set imaginary part of the highest harmonic to 0 after multiplication by ik, as long as sinusoidal harmonic is invisible. As Matt suggested, one can simply use 'symmetric' option in ifft routine

Upvotes: 0

Matt
Matt

Reputation: 2802

I think you're just missing a 2 pi in there. Here an example using x = 2 cos(2*pi*t)

Tp = 10; % sample length
deltaTime = Tp / 200; % time step
time = 0:deltaTime:Tp; % time
x = 2*cos(2*pi*time); % function
plot(time, x)
fMax = 1/deltaTime/2; % maximum frequency
fMin = 1/Tp; % lowest observable frequency
xfft = fft(x) ./ (length(time)/ 2); % fft scaled to original amplitude, in case you want to plot it.
freq = -1*fMax:fMin:fMax; % frequencies of the fft
xd_fft = xfft .* fftshift(freq) * 1i*2*pi; % note the extra 2 pi in here.
xd = ifft(xd_fft, 'symmetric') * (length(xd_fft)/ 2); % reverse the scaling and take the ifft.
xd2 = -4*pi*sin(2*pi*time);
plot(time, xd2, '.')

The equation that I have for doing this is:

x(t) = Xe^(iwt) and

x'(t) = iwXe^(iwt)

Recall that w = 2*pi*f.

If I knew how to use greek symbols on SO I would. But I think you get what I'm saying.

Upvotes: 1

Related Questions