%% Filtering in Matlab tutorial 2 - Ellen Lau, Nov. 7 2006

%% In this tutorial we examine the properties of several different filters
%% in the matlab package. without knowing too much technical information
%% about them before hand, we can discover most of what we need to
%% know about their application and effects on our data by looking at
%% 1)their frequency response function, 2) their effect on our data in the
%% frequency domain, 3) their effect on our data in the time domain. 

%% Let's get back our 50 + 120 hz component wave from the last problem:

Fs = 1000;                    % Sampling frequency
T = 1/Fs;                     % Sample time
L = 1000;                     % Length of signal
t = (0:L-1)*T;                % Time vector
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)');

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


%% The way filtering works in Matlab is, you create the filter and get back
%% the coefficients, and then you do a separate command to do the actual
%% filtering of the waveform.  

%% Let's look at how this would work in a Butterworth filter, which is
%% known as one of the more gradual cutoff filters

[b,a] = butter(9,100/500,'low');

%% given a sampling rate of 1000, this command returns a 9th order lowpass digital Butterworth filter with normalized cutoff frequency 100 Hz. 
%% It returns the filter coefficients in length n+1 row vectors b and a, with coefficients in descending powers of z.

%% The cutoff frequency is that frequency where the magnitude response of the
%% filter is sqrt(1/2), which means it's a *power* cutoff

%% The order basically means, how many surrounding points are used to
%% compute the current point

%% We can look at the frequency response of this filter by using the freqz command

figure;
freqz(b,a,128,1000)

%% The top plot is change in amplitude / change in frequency. can see that
%% the drop-off is very gradual. the bottom is phase, we won't worry about
%% this too much yet, but basically this shows us that there's a rapid
%% increase in the amount of phase shift with increase in frequency, and
%% bad for us (if we want to keep our 50 Hz component but get rid of our
%% 120 Hz component), there's a significant amount of phase shift at 50 Hz

%% Let's look at how this filter impacts our 50 + 120 pure signal

xf = filter(b,a,x);

Y = fft(xf,NFFT)/L;    % calculate the fast fourier transform and scale it to length
f = Fs/2*linspace(0,1,NFFT/2); 

% Plot single-sided amplitude spectrum.

figure;
subplot(2,1,1),plot(Fs*t(1:50),xf(1:50)),title('x(t) + 100 Hz 9th order Butterworth filter');
subplot(2,1,2),bar(f,2*abs(Y(1:NFFT/2))),title('Single-Sided Amplitude Spectrum of x(t) + 100 Hz 9th order Butterworth filter'),xlabel('Frequency (Hz)'),ylabel('|Y(f)|');

%% We can see that although there's some mess, we're mostly left with our
%% 50 Hz signal (same 2.5 peaks in this period as with the waveform we
%% looked at earlier of the pure 50 Hz component), but notice the phase
%% shift when compared to the pure one. In the frequency domain we see that
%% the 120 hz power has mostly been attenuated, but there's still a not
%% insubstantial amount left.

%% Now you may want to play around with the parameters of this filter and
%% see how it changes your signal.

%% Now let's look at another kind of filter in Matlab: elliptic or 'Cauer'
%% filters. They tend to have steeper rolloffs than a Butterworth, and they
%% give you more control over the output of the filter--you decide how much
%% ripple you want in the passed frequencies and how much attenuation you
%% want in the stopped band.

[b,a] = ellip(6,3,50,100/500);

% this filter is a sixth order lowpass filter with a passband edge
% frequency of 100 Hz (normalized to db is .6), 3 dB of ripple in the
% passband, and 50 dB of attenuation in the stopband

figure;
freqz(b,a,512,1000)

%% note that the cutoff is basically as good as the butterworth, even
%% though the order of the filter is lower, which gives you slightly less
%% phase shift. elliptical filters tend to give better performance at lower
%% orders. 


%% Let's look at how this elliptical filter impacts our 50 + 120 pure signal

xf = filter(b,a,x);

Y = fft(xf,NFFT)/L;    % calculate the fast fourier transform and scale it to length
f = Fs/2*linspace(0,1,NFFT/2); 

% Plot single-sided amplitude spectrum.

figure;
subplot(2,1,1),plot(Fs*t(1:50),xf(1:50)),title('x(t) + 100 Hz 6th order elliptical filter');
subplot(2,1,2),bar(f,2*abs(Y(1:NFFT/2))),title('Single-Sided Amplitude Spectrum of x(t) + 100 Hz 6th order elliptical filter'),xlabel('Frequency (Hz)'),ylabel('|Y(f)|');


%% the waveform is still phase-shifted about the same amount as the
%% butterworth, which is bad, but you have a smoother waveform, due to the
%% improved attenuation of the 120 hz component visible in the spectrogram.
%% it also attenuates the passed frequency a bit more than the
%% butterworth however.

%%if we do the same filter with a lower order, we get less phase shift, but
%%more of the 120 hz component creeps back in, and you get a bit more still
%%attenuation in the passed frequency

[b,a] = ellip(3,3,50,100/500);
figure;
freqz(b,a,512,1000)
xf = filter(b,a,x);

Y = fft(xf,NFFT)/L;    % calculate the fast fourier transform and scale it to length
f = Fs/2*linspace(0,1,NFFT/2); 
figure;
subplot(2,1,1),plot(Fs*t(1:50),xf(1:50)),title('x(t) + 100 Hz 3rd order elliptical filter');
subplot(2,1,2),bar(f,2*abs(Y(1:NFFT/2))),title('Single-Sided Amplitude Spectrum of x(t) + 100 Hz 3rd order elliptical filter'),xlabel('Frequency (Hz)'),ylabel('|Y(f)|');


%% another possibility is to make the half-power cutoff even lower, say at 80, so you get
%% greater attenuation in the 120 hz frequency, without a huge cost to the
%% phase

[b,a] = ellip(3,4,50,80/500);
figure;
freqz(b,a,512,1000)
xf = filter(b,a,x);

Y = fft(xf,NFFT)/L;    % calculate the fast fourier transform and scale it to length
f = Fs/2*linspace(0,1,NFFT/2); 
figure;
subplot(2,1,1),plot(Fs*t(1:50),xf(1:50)),title('x(t) + 80 Hz 3rd order elliptical filter');
subplot(2,1,2),bar(f,2*abs(Y(1:NFFT/2))),title('Single-Sided Amplitude Spectrum of x(t) + 80 Hz 3rd order elliptical filter'),xlabel('Frequency (Hz)'),ylabel('|Y(f)|');

%% however, you have a bit more still attenuation of the passed frequency



%% the above examples illustrate the inescapable tradeoff of precision in
%% time (reduced phase shift) and precision in frequency (relative power of
%% frequencies that are passed or stopped). this is why it is important to
%% assess beforehand in which domain precision is more important in your
%% analysis