% Testing the fast Fourier transform
%
% Samuli Siltanen November 2014

% Construct evaluation points for the function
M = 512;
xvec = ([0:(M-1)]/M)*2*pi;
Dx   = xvec(2)-xvec(1);

% Evaluate the function
g = MyFun3(xvec);
%g = exp(i*5*xvec);

% Compute the FFT
Fg = fft(g);

% Reorganize the FFT elements
Fg = fftshift(Fg);

% Put some of the frequency components to zero
nvec        = [[-M/2:-1],[0:(M/2-1)]];
Fg_filtered = Fg;
Fg_filtered(abs(nvec)>2) = 0;
g_filtered  = ifft(fftshift(Fg_filtered));

figure(1)
clf

% Plot the original function and the FFT
subplot(2,1,1)
plot(xvec,g,'k')
hold on
plot(xvec,g_filtered,'r')
axis([0 2*pi -1.5 1.5])
subplot(2,1,2)
plot(nvec,abs(Fg),'k')
hold on 
plot(nvec,abs(Fg_filtered),'r')
xlim([-M/2 M/2-1])

