%% FFTGIBBS.M   Apr 23, 2007:
%% To show Gibbs phenomenon, using FFTs to sum Fourier series.
%% First it shows three curves superposed, for 10, 30 and 100 terms;
%% then four graphs on a page, for 2, 3, 4 and 6 terms.
%
clear; clg; hold off;
% Set up +/- square pulse
 nn = 2^13;
 x(1:nn) = ones(1,nn);
 x(nn/2+1:nn) = -ones(1,nn/2);

% take its FT and blot out high frequencies: NOTE wrap-round, so that
% I blot out all but the first j and last j points.
 y = fft(x);
 j=10;
 y(1+j:nn-j) = zeros(1, nn-2*j);

% then take the Inverse fft and plot the result
 z = ifft(y);
 plot(real(z), 'y')
 hold on;

% repeat with j=30
 y = fft(x);
 j=30;
 y(1+j:nn-j) = zeros(1, nn-2*j);
 z = ifft(y);
 plot(real(z), 'c')

 % repeat with j=100
 y = fft(x);
 j=100;
 y(1+j:nn-j) = zeros(1, nn-2*j);
 z = ifft(y);
 plot(real(z), 'r')
 hold off;

pause;

 % do the plots again, on separate graphs this time
 nn = 2^13;
 x(1:nn) = ones(1,nn);
 x(nn/2+1:nn) = -ones(1,nn/2);
%% vals = [3, 10, 30, 100];
 vals = [2, 3, 4, 6];
 z = fft(x);

 for k=1:4
 y = z;
 j=vals(k);
 y(1+j:nn-j) = zeros(1, nn-2*j);
 y = ifft(y);
 subplot(2,2,k), plot(real(y), 'y'); title([int2str(j) ' terms']);
 end;