clear all, close all, clc

% create a simple signal with 2 frequencies
dt = .002; % sampling interval
Fs = 1/dt; % sampling rate 500 Hz
%shannon nqyuist theorem Fs/2 = 500Hz and beyond cannot be fully
%captured resolved by Fs = 1000 Hz

t = 0:dt:2-dt; % second 
f1 = 50; %50 Hz
f2 = 120; % 120 Hz
fclean = sin(2*pi*f1*t) + sin(2*pi*f2*t);
f = fclean + 2.5*randn(size(t)); % add large gaussian white noise

% plot clean signal vs noisy signal together
figure; set(gcf, 'Position', [1500 200 2000 2000])
subplot(3, 1, 1)
plot(t, f, 'c', 'LineWidth', 3), hold on
plot(t, fclean, 'k', 'LineWidth', 2.5)
l1 = legend('Noisy', 'Clean'); set(l1, 'FontSize', 32)
ylim([-10 10]); set (gca, 'FontSize', 32)

%"could we tell real clean signal and denoise?" 

% compute FFT
n = length(t);
fhat = fft(f, n); % one liner :) How convenient! sample size n
% f --> f_hat (each entry is complex number. Magnitude and phase)
% Magnitue : How important that freq is 
% Phase : more sin or cos?
freq = 1/(dt*n)*(0:n-1); % bin number to -> actual frequency in Hz mapping

% figure;
% stem(freq, abs(fhat));
% xlabel('Frequency in Hz')
% ylabel('abs of each f_hat value')
% 
% figure;
% stem([0:n-1], abs(fhat));
% xlabel('Bin Number')
% ylabel('abs of each f_hat value')


PSD = fhat.*conj(fhat)/n; % just getting the magnitude

% single side
L = 1:floor(n/2); % just plot the first half of freqs
subplot(3,1,3); set(gca, 'FontSize', 16)
plot(freq(L), PSD(L), 'c', 'LineWidth', 3), hold on
% set(gca, 'FontSize', 32)
% 
% 
% % c x c_bar = getting magnitude of 
% % c = a+ib 
% % c_bar = a - ib
% % |c|^2 = a^2 + b^2 
% 
% % Quite clear to see have peaks!! 
% % Can we filter??
% 
threshold = 100;
indicies = PSD >= threshold; % vectors of 0 or 1, (0 if less, 1 if greater)
%PSDclean = PSD.* indicies;
fhat_filtered = indicies.*fhat;
f_filtered = ifft(fhat_filtered);

subplot(3,1,2); set(gca, 'FontSize', 16)
plot(t, f_filtered, '-', 'Color', [.5 .1 0], 'LineWidth', 2.5)
l1 = legend('Filtered'); set(l1, 'FontSize', 32)
ylim([-10 10]); set(gca, 'FontSize', 32)




