close all; clear; clc;

%%
N = 21; t = linspace(0,1,N+1); t(end) = []; % time axis / s
y1 = @(x) sin(2*pi*3*x); y2 = @(x) sin(2*pi*4.5*x);

% 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

% Window functions and names
windows = {@(N) boxcar(N)', @(N) hann(N)', @(N) flattopwin(N)'};
names = {'Boxcar','Hann','Flattop'};

dtft = @(y,v) 1/N*sum(bsxfun(@times,y,exp(2*pi*1i.*v'.* (t-t(1)))),2)';
v_dense = linspace(min(v),max(v),1e3);
t_dense = linspace(min(t),max(t),1e3);

%% Plot
figure(1); clf;
for j = 1:3
    window = windows{j}; % Window
    w = window(N); w = w./mean(w); % Scale to conserve signal power
    w_dense = window(numel(t_dense));
    w_dense = w_dense./mean(w_dense);
    
    Y1 = fftshift(fft(y1(t).*w))/N;
    Y2 = fftshift(fft(y2(t).*w))/N;
    Y1_dense = dtft(y1(t).*w,v_dense); 
    Y2_dense = dtft(y2(t).*w,v_dense); 
    
    s1 = subplot(2,3,j); hold on;
    plot(t_dense,y1(t_dense).*w_dense,'k-');
    plot(t,y1(t).*w,'o');
    %plot(t_dense,y2(t_dense).*w_dense,'-'); plot(t,y2(t).*w,'k.');
    if j ~= 1
        plot(t_dense,w_dense,'k:');
        plot(t_dense,-w_dense,'k:');
    end
    grid on; box on; axis square; xlabel('t'); xlim([min(t) max(t)]);
    title(names{j});
    if j == 1
        ylabel('y*w');
    else
        s1.YTickLabel = '';
    end
    ylim([-1 1]*5);
    drawnow;
    
    s2 = subplot(2,3,j+3); hold on;
    plot(v_dense,abs(Y1_dense),'k-');
    plot(v,abs(Y1),'x');
    plot(v_dense,abs(Y2_dense),'k:');
    plot(v,abs(Y2),'o');
    grid on; box on; axis square; xlabel('v / Hz');
    if j == 1
        ylabel('Y');
    else
        s2.YTickLabel = '';
    end
    ylim([0 0.6]); xlim([-5 5]+3);
    drawnow;
end