close all; clear; clc;

%% Sub and super-sampled signal
N = 13; t = linspace(0,1,N+1); t(end) = []; dt = t(2)-t(1);
fs = 1/dt; % Sampling frequency
%f = [1 2.2 3 4]' + randn(4,1);
f = [0.79 2.11 3.05 4.11]'; % Signal frequency components
y1 = @(x) sum( sin(2*pi*x.*f) ,1);
y2 = @(x) sum( sin(2*pi*x.*(fs-f)) ,1);

t_dense = linspace(0,1,1e3);

% Plot 1
figure(1); clf;
subplot(1,2,1);
plot(t_dense,y1(t_dense),'k'); hold on;
plot(t_dense,y2(t_dense),'k:'); hold on;
plot(t,y1(t),'o');
plot(t,y2(t),'o');
grid on; box on; axis square;
xlabel('time / s'); ylabel('signal');

%% Nyquist reflection surface
clear YY
N = 51*2+1; t = linspace(0,1,N+1); t(end) = [];
dt = t(2)-t(1); fs = 1/dt;
y = @(t,f) cos(2*pi*t.*f) + cos(2*pi*t.*(f-fs/4));

fqs = linspace(fs/4, 5*fs/4, ceil((N+1)/2)); % Signal frequencies
for fi = 1:numel(fqs)
    Y = fft( y(t,fqs(fi)) )/N; Y = Y(1:ceil((N+1)/2)); % Half spectrum
    Y(2:ceil(N/2)) = Y(2:ceil(N/2))*2; % Adjust amplitude of middle samples
    YY(fi,:) = abs(Y); % f x v matrix
end

vh = linspace(0, (1-mod(N,2)/N)/2, ceil((N+1)/2)); vh = vh/dt;

% Plot 2
a = subplot(1,2,2);
imagesc(vh,fqs,log10(YY)); grid on;
axis tight square;
colormap(viridis);
xlabel('FT frequency bin / Hz'); ylabel('signal frequency / Hz');
a.YDir = 'normal';