close all; clear; clc;

%%
rng(45);

N = 2^10; t = linspace(0,1,N+1); t(end) = []; dt = t(2)-t(1);
y = sin(2*pi*t.*2) + randn(size(t)); % Signal + noise
v = linspace(0,1-1/N,N) + (mod(N,2)/N-1)/2; v = v/dt;

% Decimate signal
naive = y(1:64:N);
aaliased = decimate(y,64);
N2 = N/64;
v2 = linspace(0,1-1/N2,N2) + (mod(N2,2)/N2-1)/2; v2 = v2/(64*dt);

%% Plot
t2 = t(1:64:N);

cm = lines;
figure(1); clf;
subplot(1,2,1); hold on;
plot(t,y,'Color',[cm(1,:) 0.6]);
plot(t2,naive,'o-','Color',cm(2,:),'Linewidth',1,'MarkerSize',4);
plot(t2,aaliased,'o-','Color',cm(3,:),'Linewidth',1,'MarkerSize',4);
grid on; box on; axis square;
xlabel('time / s'); ylabel('signal');

subplot(1,2,2);
hold on;
h(1) = stem(v,abs(fftshift(fft(y)))/N,'MarkerSize',4);
h(2) = stem(v2,abs(fftshift(fft(naive)))*64/N,'MarkerSize',4);
h(3) = stem(v2,abs(fftshift(fft(aaliased)))*64/N,'MarkerSize',4);
grid on; box on; axis square;

xlim([-20 20]); ylim([0 0.6]);
xlabel('frequency / Hz'); ylabel('fft( signal )');
grid on
