matlabfft

Slight error with created inverse DFT function in MATLAB

I am trying to write my own inverse DFT function in MATLAB as I need a subset of these Fourier mode basis in real space for an application. However, when I check my inverse DFT routine after extracting the Fourier coefficients against Matlab's ifft function and the original waveform there seems to be a slight error. The code is seen from below.

% 1D fft and inverse fft
clear; clc; close all;
L=1; % total domain size
n=100;
x=linspace(0,L,n); %define x values
dx=x(1)-x(2);
F=1*cos(2*pi*5*x); %declare the function
figure(1)
plot(x,F,'linewidth',2,'Displayname','Original');
Fhat=fftshift(fft(F));%compute fft and shift to center
if(mod(n,2)==0) % even split
    dk=2*pi/L.*[-n/2:1:1:n/2-1];
else %odd split
    dk=2*pi/L*[-(n-1)/2:1:(n-1)/2];
end
figure % plot FFt results
plot(dk/(2*pi),abs(Fhat),'displayname','Frequency')
% Now here is my inverse fft
f=zeros(1,n);
for i=1:n
    f=f+Fhat(i).*exp(1j*dk(i).*x)
end
f=f*1/n;
figure(1); hold on;
plot(x,real(f),'b-.','linewidth',2,'displayname','From my own invDFT')
plot(x,ifft(ifftshift(Fhat)),'M--','linewidth',2,'displayname','From matlab ifft')
xlabel('x','fontsize',16)
ylabel('y','fontsize',16)
legend()

Below is an image showing the discrepancy between the waveforms. Comparison of methods

I'm not sure if this issue is a floating point problem or if my frequency bins are incorrect. This is my first time using fft. Any insight is greatly appreciated.

Solution

  • TL;DR

    Your code is calculating the inverse as

    enter image description here

    instead of

    enter image description here

    But why?

    Well, first of all your signal is periodic of period L, so when you use a vector x going all the way from 0 to L you overlap onto the next period of the signal, which will induce confusion (and aliasing).

    This said, the DFT (and its inverse) do not explicitly use frequencies, they use bins indexed on the length n of the signal:

    enter image description here

    replicating those bins is going to lower the risk of confusion, in your case:

    clear; clc; close all;
L=1; % total domain size
n=100;
x=linspace(0,L,n); %define x values
F=1*cos(2*pi*5*x); %declare the function
figure(1)
plot(x,F,'linewidth',2,'Displayname','Original');
Fhat=fftshift(fft(F));%compute fft and shift to center
if(mod(n,2)==0) % even split
    dk=2*pi/n.*[-n/2:1:1:n/2-1]; % <-- using n
else %odd split
    dk=2*pi/n*[-(n-1)/2:1:(n-1)/2]; % <-- using n
end
figure % plot FFt results
% plot(dk/(2*pi),abs(Fhat),'displayname','Frequency')
% Now here is my inverse fft
f=zeros(1,n);
for i=1:n
    f=f+Fhat(i).*exp(1j*dk(i).*(0:(n-1))) % <-- critical difference
end
f=f*1/n;
figure(1); hold on;
plot(x,real(f),'b-.','linewidth',2,'displayname','From my own invDFT')
plot(x,ifft(ifftshift(Fhat)),'M--','linewidth',2,'displayname','From matlab ifft')
xlabel('x','fontsize',16)
ylabel('y','fontsize',16)
legend()

    The critical difference is that the exponent j-1 in the previous equation does not go all the way to n but only to n-1, contrary to what your x/L was doing. It had the same effect as resampling all your frequencies by a factor n/(n-1), hence the slight shift of frequency you observed.

    The result is now coherent:

    enter image description here