%see Matlab documentation
%doc fft
N = 1000;
Fs = N;                 % Sampling frequency
T = 1/N;                   % Sample time
L = N;                     % Length of signal
t = (0:L-1)*T;             % Time vector

% Define the true signal
% Sum of a 50 Hz sinusoid and a 120 Hz sinusoid
x = sin(2*pi*50*t) + 2*sin(2*pi*120*t); 
figure(1)
plot(Fs*t(1:100),x(1:100),'r','linewidth',3)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
hold on
pause(5)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%introduce noise in the signal
noise = 2*randn(size(t));
y = x + noise;     % Sinusoids plus noise
plot(Fs*t(1:100),y(1:100),'k','linewidth',2)
title('Signal Corrupted with Zero-Mean Random Noise','fontsize',16)
xlabel('time (milliseconds)','fontsize',16)
legend('true signal','noisy signal')
pause(10)

%perform the Fourier transform fft
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2); %frequencies axis 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Plot single-sided amplitude spectrum (positive values)
figure(2)
plot(f,2*abs(Y(1:NFFT/2)),'.','markersize',20) 
title('Single-Sided Amplitude Spectrum of y(t)','fontsize',16)
xlabel('Frequency (Hz)','fontsize',16)
ylabel('|Y(f)|','fontsize',16)
