close all; clear; clc;

%%
N = 61; % Number of sampled points
t_end = 5;
t = linspace(0,t_end,N+1); t(end) = []; % time axis / s

f = 2; % Frequency / Hz
tc = 1;
y = @(x) sin(2*pi*f*x) .* exp(-x./tc); % Signal
% Exact Fourier transform of the signal
%Y = @(v) (2*f*pi) ./ (4*f.^2*pi^2 - (1i-v*2*pi).^2);
Y = @(v) 2*pi*f*tc ./ (4*f^2*pi^2*tc^2 - (1i-2*pi*tc*v).^2);

%% Construct frequency axis
dt = t(2)-t(1); % N/max(t) Time sampling
v = linspace(0,1-1/N,N) + (mod(N,2)/N-1)/2; % Frequency indices
v = v/dt; % Frequency axis / Hz

%% Calculate fourier transforms
fast = fftshift(fft(y(t)))/N;
naive = 1/N*sum(bsxfun(@times,y(t),exp(2 * pi * 1i .* v' .* (t-t(1)))),2)';
s = ( tc.*(1 - exp(-t_end/tc)) ) / t_end;
exact = Y(v).*s; % Scale to account for amplitude drop 

%% Plot
% Dense axis for nice plots
t_dense = linspace(0,2*pi,1e3); v_dense = linspace(min(v),max(v),1e5);

figure(1); clf;
subplot(2,1,1);
h = plot(t_dense,y(t_dense),'k-',t,y(t),'o');
grid on; box on; xlabel('t'); ylabel('y'); xlim([min(t) max(t)]);
title('Signal'); legend(h,'Signal','Sampled points');

subplot(2,1,2);
h = plot(v,abs(fast),'x',v,abs(naive),'o',v,abs(exact),'d'); hold on;
%plot(v_dense,abs(Y(v_dense))/t_end*tc/2,'k-');
grid on; box on; xlabel('f'); ylabel('FT\{ y \}');
title('Fourier Transform'); legend(h,'FFT','Naive','Exact');
