%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% clean up the workspace
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all ; close all ; clc;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% define a sum of sinusoids
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N = 1000;
t = linspace(0,1,N)';
[A1,w1] = deal(3,20);
[A2,w2] = deal(1,60);
x = A1 * cos(2*pi * w1 * t) + A2 * cos(2*pi * w2 * t);

figure();
plot(t , x , 'k');
ylim(10 * [-1 +1]);
grid on;
xlabel('t');
ylabel('x');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% add noise to the sinusoid
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
randn('seed' , 0);
sigma = 4;
y = x + sigma * randn(N,1);

figure();
plot(t , y , 'k');
ylim(10 * [-1 +1]);
grid on;
xlabel('t');
ylabel('y');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute the DFT (note: the FFT is much more efficient)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x_dft = nan(N,1);
y_dft = nan(N,1);
n = (0:(N-1))';
for k = 0:(N-1)
    x_dft(k+1) = (1/N) * sum(x .* exp(-2*pi*1i * k*n/N));
    y_dft(k+1) = (1/N) * sum(y .* exp(-2*pi*1i * k*n/N));
end

figure();
plot(0:(N-1) , abs(x_dft) , 'k');
xlim([0 100]);
grid on;
xlabel('f');
ylabel('X');

figure();
plot(0:(N-1) , abs(y_dft) , 'k');
xlim([0 100]);
grid on;
xlabel('f');
ylabel('Y');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% filter in the frequency domain
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
z_dft = y_dft;
z_dft(z_dft < 0.3) = 0;

figure();
plot(0:(N-1) , abs(z_dft) , 'k');
xlim([0 100]);
grid on;
xlabel('f');
ylabel('Z');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute the inverse DFT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
z = nan(N,1);
for k = 0:(N-1)
    z(k+1) = real(sum(z_dft .* exp(2*pi*1i * k*n/N)));
end

figure();
plot(t , z , 'k');
ylim(10 * [-1 +1]);
grid on;
xlabel('t');
ylabel('z');