% Sampling frequency 

Fs = 1000; 


% Time vector of 1 second 

t = 0:1/Fs:1; 


% Create a sine wave of 200 Hz.

x = sin(2*pi*t*200) + sin(2*pi*t*300) + 2*randn(1,length(t)); 


% Use next highest power of 2 greater than or equal to length(x) to
% calculate FFT. (technical step, not strictly necessary, just speeds up the fft operation)
nfft= 2^(nextpow2(length(x))); 


% Take fft, pad with zeros so that length(fftx) is equal to nfft  (again, not totally necessary) 

fftx = fft(x,nfft); 


% Calculate the number of unique points
NumUniquePts = ceil((nfft+1)/2); 


% FFT is symmetric, throw away second half 

fftx = fftx(1:NumUniquePts); 


% Take the magnitude of fft of x and scale the fft so that it is not a function of the length of x

mx = abs(fftx)/length(x); 


% Take the square of the magnitude of fft of x. 

mx = mx.^2; 


% Since we dropped half the FFT, we multiply mx by 2 to keep the same energy.
% The DC component and Nyquist component, if it exists, are unique and should not be multiplied by 2.


if rem(nfft, 2) % odd nfft excludes Nyquist point
  mx(2:end) = mx(2:end)*2;
else
  mx(2:end -1) = mx(2:end -1)*2;
end


% This is an evenly spaced frequency vector with NumUniquePts points. 

f = (0:NumUniquePts-1)*Fs/nfft; 


% Generate the plot, title and labels. 
subplot(2,1,1), plot(t,x); xlabel('time (seconds)'); ylabel('signal');
subplot(2,1,2), plot(f,mx); xlabel('Frequency (Hz)');  ylabel('Power'); 