close all; clear; clc;

%%
N = 201; t = linspace(0,1,N+1); t(end) = []; dt = t(2)-t(1); % Time axis
v = linspace(0,1-1/N,N) + (mod(N,2)/N-1)/2; v = v/dt; % Frequency axis
vh = linspace(0, (1-mod(N,2)/N)/2, ceil((N+1)/2)); vh = vh/dt;

rng(42);
% Make an interesting signal
signal = sin(2*pi*[t(t<1/3)*20 t(t>=1/3 & t< 1/2)*60 t(t>=1/2)*40]);
signal = signal.*(sin(2*pi*t*3).^2) + 0.2*randn(size(signal));

d1 = 0.05; d2 = 0.4; % Widths of the window
w = @(c,d) exp( -(4*log(2).*((t-c)./d).^2 )); % Gaussian window

ts = linspace(0,1,N);
for ti = 1:numel(ts)
    Y = fft( signal.*w(ts(ti),d1) )/N; YY1(ti,:) = Y(1:ceil((N+1)/2));
    Y = fft( signal.*w(ts(ti),d2) )/N; YY2(ti,:) = Y(1:ceil((N+1)/2));
end

% Plot
cm = lines;
figure(1); clf

s1 = subplot(1,3,1); hold on;
plot(t,signal,'k');
plot(t,w(1/3,d1)-2,'color',cm(1,:),'LineWidth',1.5);
plot(t,w(2/3,d2)-2,'color',cm(2,:),'LineWidth',1.5);
axis tight square; box on; grid on;
xlabel('t'); ylabel('y');
s2 = subplot(1,3,2); hold on;
colormap(viridis);
imagesc(vh,ts,abs(YY1));
s2.YDir = 'normal';
axis tight square; box on; grid on;
xlim([0 80]);
xlabel('v'); ylabel('t');

s3 = subplot(1,3,3); hold on;
colormap(viridis);
imagesc(vh,ts,abs(YY2));
s3.YDir = 'normal';
axis tight square; box on; grid on;
xlim([0 80]);
xlabel('v'); ylabel('t');
