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



%%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


%% Now let's try adding some noise to our pure waveform--adding or subtracting a random amount of noise to each point. here we're using a gaussian
%% distribution for the noise centered around 0 with a stdev of 2


y = x + 2*randn(size(t));     % Sinusoids plus noise


%% plot the waveform with noise and notice how much it corrupts the signal compared to your earlier plot of the pure signal

figure;
plot(Fs*t(1:50),y(1:50))
title('Signal Corrupted with Zero-Mean Random Noise, y(t)')
xlabel('time (milliseconds)')


%% finally, take a look at the frequency spectrum of the signal corrupted by noise


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);


figure;
% Plot single-sided amplitude spectrum.
bar(f,2*abs(Y(1:NFFT/2))) 
title('Single-Sided Amplitude Spectrum of y(t)')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')


%% notice that although the noise did a lot of damage in the time domain, in the frequency domain (as depicted by the 
%% amplitude spectrogram) the key properties of the original signal are still retrievable despite the noise.  That's because
%% the noise was basically 'white noise'--equal amounts at all frequencies. In this case, it would be harder to get rid of the
%% noise using a frequency filter, although if you knew the range where you expected to see the signal--say, under 150 Hz--you could low-pass
%% at least the random noise above this cutoff. However, the benefit of this case, all the non-randomness in the frequency distribution
%% can be correctly attributed to the signal. In other cases, like noise from nearby electrical sources, noise will have a lot of amplitude
%% at a particular frequencies, so if such sources of noise are possible, you can not assume that the frequency domain will give you a good
%% representation of the signal only.  However, if you are aware of the frequencies of the noise, you can construct a filter more easily to 
%% get rid of it.