
% 1. Create an artificial data set
deltat = .001; % time step
t = 0:deltat:1.023; % time vector
n = length(t); % we have 1024 data points
f1 = 50; f2 = 120;  % we use two frequencies: f1 and f2
% our artificial time series: two sine waves and random noise
y = sin(2*pi*f1*t)+2*sin(2*pi*f2*t)+0.5*randn(size(t));

% 2. calculate FFT and power spectrum
Y = fft(y);
Py = Y.*conj(Y); % |Y|^2

% 3. plot signal and power spectrum
Fs = 1/deltat; % sampling frequency
f = [0:n/2-1]/n*Fs; % frequency intervals for plotting
% only up to the Nyquist frequency (1/(2* deltat) = Fs/2)
Py(n/2+1:n)=[]; % chop off Fourier coeffs at and above Nyquist frequency

figure 
subplot(1,2,1)
plot(t,y)
set(gca,'FontSize',12)
xlabel('Time [s]')
ylabel('y(t)')
title('Time domain') 
axis([0 1 -4 4])
%
subplot(1,2,2)
semilogy(f,Py)
set(gca,'FontSize',12)
xlabel('Frequency [1/s]')
ylabel('|Y|^2')
title('Frequency domain')
axis([0 500 10^(-2) 10^(6)])
