%% Filtering in Matlab tutorial 1 - Ellen Lau, adapted from example on Mathworks website, Nov. 11 2006
close all
clear all
clf
clc


%%This is a frequency analysis example from the matlab website to get better understanding of doing fft (fast fourier transform)
%%No filtering, just getting comfortable with sinusoids and ffts and effects of noise

%"A common use of Fourier transforms is to find the frequency components of a signal buried in a noisy time domain signal. 
%Consider data sampled at 1000 Hz. Form a signal containing a 50 Hz sinusoid of amplitude 0.7 and 120 Hz sinusoid 
%of amplitude 1 and corrupt it with some zero-mean random noise:"

%%I designed this tutorial with a cut-and-paste approach in mind: read each comment and then cut-and-paste commands below into command window

%%%%%%%%%%%%%%%%%%%%


%%First enter the parameters of the problem: we're sampling 1000 times per time unit (probably a second) = Fs (sampling frequency)
%%so each time sample covers 1/1000 of a second = T 
%%the length of the entire signal that we have to look at is 1000 samples
%%we make a time vector that has an entry for each time sample in the whole signal time--so goes from 0 to the (signal length - 1) 
%%(because we're starting to count from zero instead of 1) and has total time samples number of entries (which should equal the length
%%which are listed in time units (so here, fraction of a second)


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


%%take a look at the sample time and the time vector we just entered, just so you know you've got a handle on them


t
T


% Now we're going to create an idealized signal, which is composed of two parts, a 50 Hz sinusoid at .7 amplitude and
% and a 120 Hz sinusoid at 1 amplitude, as a function of time. Create the signal and take a look at its components (in the first
% windows) and the summed whole signal. 


figure;

x1 = 0.7*sin(2*pi*50*t);
x2 = sin(2*pi*120*t);

subplot(3,1,1),plot(Fs*t(1:50),x1(1:50)),title('50 hz sin wave');
subplot(3,1,2),plot(Fs*t(1:50),x2(1:50)),title('120 hz sin wave');

x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);     %%create whole signal from component sinusoids

subplot(3,1,3),plot(Fs*t(1:50),x(1:50)),title('whole pure signal x(t)');

% Take a look at the frequency amplitude spectrum of the signal you just created. Remember the signal is purely a 50 Hz sinusoid
% And a 120 hz sinusoid, so in theory should only have amplitude at those two frequencies. However, since we are using a finite
% time window, we'll see small amounts of frequency leakage at other frequencies


NFFT = 2^nextpow2(L); % Next power of 2 up from length of the signal--need this for the fast fourier transform algorithm 
Y = fft(x,NFFT)/L;    % calculate the fast fourier transform and scale it to length
f = Fs/2*linspace(0,1,NFFT/2);   %%linspace creates a linearly spaced vector between 0 and 1 with NFFT/2 points evenly distributed
                                 %%this step gives you the x axis for your frequency spectrum. go up only to NFFT/2, not NFFT, b/c
                                 %%spectrum is redundant in second half
figure;
% Plot single-sided amplitude spectrum.
bar(f,2*abs(Y(1:NFFT/2))) 
title('Single-Sided Amplitude Spectrum of x(t)')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')


% Note that the frequency spectrum is basically as we expect, with most of the signal amplitude concentrated at 50 and 120, and
% just a bit of leakage at other frequencies; note also that the amplitudes at those two frequencies correspond to the amplitudes
% definied in our sine wave equations