% Illustrate Gibbs Phenomenon with Boxcar function
% GibbsBoxcar.m L Braile 03/24/06
% Calculate boxcar by FFT of a truncated (finite length) (sin x)/x
dt = 0.01; % sample interval (s)
n = 8;      % number of points in (sin x)/x [length of (sin x)/x)]
% Try n = 8, 16, 32, 64 128, to see that Gibbs effect is not 
% reduced by using a longer and longer (sin x)/x time series.
nf = 2048;    % number of frequencies if FFT
t = [dt:dt:n*dt]; % calculate n values of time
f = 15;       % frequency of (sin x)/x signal
x = sin(2*pi*t*f)./(2*pi*t*f); % calculate (sin x)/x
x1 = fliplr(x(1:(n-1)));  % make it symmetric
y = [x1 1 x];
t1 = [-(n)*dt:dt:n*dt-dt];
% ----------------------------------------------------------------
plot(t1,y,'-b','linewidth',1.5)  % Figure 1
set(gca,'fontsize',14,'linewidth',2)
xlabel('Time (s)','fontsize',14)
ylabel('Amplitude','fontsize',14)
title('N-Point (sin x)/x)','fontsize',14)
y2 = fftshift(y); % arrange order for FFT
figure
y3 = fliplr(y2);
% ----------------------------------------------------------------
plot(t1,y3,'-b','linewidth',1.5)  % Figure 2
set(gca,'fontsize',14,'linewidth',2)
xlabel('Time (s)','fontsize',14)
ylabel('Amplitude','fontsize',14)
title('N-Point (sin x)/x), Arranged for FFT','fontsize',14)
y3e = [y3(1:n) 0 zeros(1,nf-2*n) y3(n+2:2*n)]; % Add zeros to middle
                           % and make time series symmetric.
Y3 = fft(y3e,nf);  % calculate nf-point FFT of truncated (sin x)/x
rY3 = real(Y3);    % Take real part of Y3 (because of symmetry)
                   % all imaginary parts are zero.
figure
% ----------------------------------------------------------------
fnyq = 1/(2*dt);
df = fnyq/(nf/2);
f1 = [0:df:2*fnyq-df];
plot(f1,rY3,'r-','linewidth',1.5)         % Figure 3
set(gca,'fontsize',14,'linewidth',2)
xlabel('Frequency (Hz)','fontsize',14)
ylabel('Relative Amplitude','fontsize',14)
title('FFT of (sin x)/x)','fontsize',14)
rY3n = [rY3((nf/2)+2:nf) rY3(1:(nf/2)+1)];
figure
f2 = [-fnyq+df:df:fnyq];
% ----------------------------------------------------------------
plot(f2,rY3n,'r-','linewidth',1.5) % Figure 4
set(gca,'fontsize',14,'linewidth',2)
xlabel('Frequency (Hz)','fontsize',14)
ylabel('Relative Amplitude','fontsize',14)
title('FFT of (sin x)/x)','fontsize',14)
% ----------------------------------------------------------------