close all; clear; clc;

%%
N = 400; 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

rng('shuffle'); seed = round(rand*101); rng(seed); fprintf('seed = %d\n',seed);
%rng(42);

signal = sin(2*pi*t*5); % Hz baseband signal
carrier = sin(2*pi*t*80); % 80 Hz carrier

noise = 1.5*randnd(-2,1,N) + 0.5;

unmodulated = signal + noise;
modulated = signal.*carrier + noise;
demodulated = modulated.*carrier;

%% 6-pole IIR Butterworth filter with cutoff at 10 HZ
[b,a] = butter(6,10*dt*2);
lpfiltered = filtfilt(b,a,demodulated); % Zero-phase filtering
[h,~] = freqz(b,a,N,'whole'); % Calculate the filter frequency response
filterwindow = fftshift(abs(h)); % Shape of the filter window

%% Plot
w = hann(N)'; w = w./mean(w); % Hann window
cm = lines;
figure(1); clf;
s1 = subplot(1,2,1); hold on;

plot(t,w.*(unmodulated)+40,'color',cm(1,:));
plot(t,w.*(modulated)+30,'color',cm(2,:));
plot(t,w.*(demodulated)+20,'color',cm(4,:));
plot(t,w.*(lpfiltered)+10,'color',cm(5,:));
plot(t,w.*(signal),'k');
grid on; box on; axis square;
s1.YTickLabel = [];
xlabel('time / s'); ylabel('signal \times hann window');
ylim([-5 50]);

s2 = subplot(1,2,2); hold on;
plot(v,abs(fftshift(fft( w.*(unmodulated) )))/N+4,'color',cm(1,:));
plot(v,abs(fftshift(fft( w.*(modulated) )))/N+3,'color',cm(2,:));
plot(v,abs(fftshift(fft( w.*(demodulated) )))/N+2,'color',cm(4,:));
plot(v,abs(fftshift(fft( w.*(lpfiltered) )))/N+1,'color',cm(5,:));
plot(v,filterwindow*0.8+1,'k:');
plot(v,abs(fftshift(fft( w.*(signal) )))/N,'k');
xlim([-200 200]); ylim([-0.5 5.5]);
grid on; box on; axis square;
s2.YTickLabel = [];
xlabel('frequency / Hz'); ylabel('fft( signal )');