%-- SOLUTIONS to Day 2 Exercises
%-- Peter S. Shawhan

display('%----- Part 1a/1b')
%-- Read in the template
template = load('day2template.txt');
%-- Report the length of the template in seconds
template_secs = length(template) / 16384
%-- Generate an array of times (in seconds) of the samples in the template
ttemplate = (0:length(template)-1)/16384;
%-- Plot the template vs. time
plot(ttemplate,template)
title('(1b) Template')
grid on
%-- By counting cycles per second on the plot, I estimate the starting
%--     frequency to be around 40 Hz
input('(1b) Judge starting frequency by counting cycles per sec; Press RETURN to continue');

%-- Zoom in on the end of the template to check the ending frequency
plot(ttemplate(13000:16384),template(13000:16384));
title('(1b) Zoomed in on end of template')
grid on
%-- I judge the ending frequency to be around 200 Hz
input('(1b) Judge ending frequency by counting cycles per sec; Press RETURN to continue');

display('%----- Part 1c')
%-- Load in the data time series
data = load('day2data.txt');
%-- Report the length of the data in seconds
data_secs = length(data) / 16384
%-- Generate an array of times (in seconds) of the samples in the data
tdata = (0:length(data)-1)/16384;
%-- Plot it
plot(tdata,data);
title('(1c) Original data time series')
%-- No signal evident in the plot
input('(1c) No signal evident; Press RETURN to continue');

display('%----- Part 1d')
%-- Plot just the first 10000 data points
plot(tdata(1:10000),data(1:10000));
title('(1d) First 10000 points of original time series')
input('(1d) Low-frequency noise is visible; Press RETURN to continue');

display('%----- Part 1e')
%-- Calculate the PSD.  Note how we tell pwelch what the sampling frequency is.
pwelch(data,[],[],[],16384)
input('(1e) PSD of the data; Press RETURN to continue');

%-------------------------------
display('%----- Part 2a')
%-- Apply the simple difference filter
filtdata = data - [0 data(1:length(data)-1)];

display('%----- Part 2b')
%-- Plot it
plot(tdata(1:10000),filtdata(1:10000));
title('(2b) First 10000 points of time series after simple filter')
%-- The filtered data looks more like white noise
input('(2b) No low-frequency noise evident any more; Press RETURN to continue');

display('%----- Part 2c')
pwelch(filtdata,[],[],[],16384);
%-- The filter suppressed low frequencies and increased high frequencies some.
%-- One can work out the transfer function of the filter analytically;
%-- it has magnitude 2*sin(pi*f/16384) (and also introduces a phase shift,
%--     but that doesn't affect the power spectrum)
input('(2c) Low freq noise suppressed, high freq increased some; Press RETURN to continue');

%-------------------------------
display('%----- Part 3a')
%-- Design the filter, and apply it.  Note that you have to pass '4' as the
%-- first argument to get an 8th-order bandpass filter - a quirk of Matlab
[b,a] = butter(4,[35/8192 250/8192]);
%-- Filter the data
filtdata = filter(b,a,data);

display('%----- Part 3b')
plot(tdata,filtdata);
title('(3b) Time series after Butterworth bandpass filter')
%-- A signal appears in the filtered data (larger than the noise) at t~6
input('(3b) A signal is visible; Press RETURN to continue');

%-------------------------------
display('%----- Part 4a')
%-- Create an empty array to store all the correlation results.  If you don't
%-- do this, the loop spends most of its time extending the results array!
corr = zeros(1,length(data)-length(template));
%-- Now loop over time offsets
for ioff=0:length(data)-length(template);
  %-- There are a few different ways to take the dot product
  %- This version takes the sum of an element-by-element product vector:
%%  corr(ioff+1) = sum( data(ioff+1:ioff+length(template)).*template' );
  %- This version uses the matrix product of a row vector with a column vector:
  corr(ioff+1) = data(ioff+1:ioff+length(template))*template;
end

display('%----- Part 4b')
plot(corr);
title('(4b) Correlation, i.e. output from matched filter in time domain')
%-- Now a clear spike is seen!  Maybe also a smaller one near the right
%--    side of the plot...
input('(4b) A clear signal spike is seen; Press RETURN to continue');

%-------------------------------
% Part 5

display('%----- Part 5a')
rms_orig = sqrt( mean(data.^2))

display('%----- Part 5b')
%-- Generate an array with integers spaced by 16
sel = (1:16:length(data));
dsdata = data(sel);
%-- Another way to do that:
% dsdata = downsample(data,16);

%-- Calculate the rms
rms_simply_decimated = sqrt( mean(dsdata.^2))
%-- This rms is essentially the same as before, because taking a subsample of
%-- the original data points does not change their characteristic spread

display('%----- Part 5c')
[b,a] = butter(6,512/8192,'low');
            %-- The 'low' argument is optional -- makes lowpass by default
filtdata = filter(b,a,data);
%-- Select every 16th sample
dsdata = filtdata(sel);
plot(tdata(sel),dsdata);
title('(5c) Downsampled from original time series by factor of 16')
input('(5c) Plot of properly down-sampled data; Press RETURN to continue');

display('%----- Part 5d')
rms_properly_downsampled = sqrt( mean(dsdata.^2))
%-- The new RMS is somewhat smaller than the original, because we filtered out
%-- the high-frequency parts of the original time series

display('%----- Part 5e')
%-- Note that we give pwelch the NEW sampling rate
pwelch(dsdata,[],[],[],1024);
input('(5e) PSD of downsampled data, fNyquist=512 Hz; Press RETURN to continue');

display('%----- Part 5f')
%-- To down-sample the template, it shouldn't be necessary to filter first
%-- because there shouldn't be any signal power at high frequency.  I will
%-- go ahead and filter, but I would expect the results to be the same if
%-- I didn't filter.
filttemplate = filter(b,a,template);
tsel = (1:16:length(template));
dstemplate = filttemplate(tsel);

%-- Do the matched filter on the downsampled data and template
for ioff=0:length(dsdata)-length(dstemplate);
  dscorr(ioff+1) = sum( dsdata(ioff+1:ioff+length(dstemplate)).*dstemplate' );
end

plot(dscorr);
title('(5f) Correlation after downsampling by a factor of 16')
input('(5f) Output from time-domain matched filter after down-sampling; Press RETURN to continue');
