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

%-- I generally like to put constants into descriptive variables.  It is also
%-- very handy if I decide to change the parameters later.

samprate = 16384;
nyquist = samprate/2;

%-------------------------------
% Part 1

display('%----- Part 1a')
nsec = 1.0;
ndata = nsec*samprate;
t = (0:ndata-1)/samprate;
%-- Note that the array has exactly ndata samples.  Other ways to do that:
%  t = 0:(1/ndata):(ndata-1)/ndata;
%  t = 0:1/ndata:1-1/ndata;
%  t = linspace(0,(ndata-1)/ndata,ndata);

display('%----- Part 1b')
%-- Generate sine waves - remember the factors of 2pi !
twosines = 5*sin(2*pi*154.7*t) + 2*sin(2*pi*820*t);

display('%----- Part 1c')
noise = 1.4*randn(1,ndata);

display('%----- Part 1d')
timeseries = twosines + noise;
plot(t,timeseries);
title('Time series for part 1');
xlabel('Time (s)');
ylabel('Signal (arbitrary units)');
input('Whole time series; Press RETURN to continue');
%-- Here are a few ways to zoom
set(gca,'XLim',[0 0.05]);
xlim([0 0.05]);
input('Zoomed in on the first 0.05 sec; Press RETURN to continue');

display('%----- Part 1e')
timeseries_mean = mean(timeseries)
%-- There are a few different ways to calculate the RMS - see below.
%-- Note that by "rms" we usually mean "rms deviation from the mean", so the
%-- "correct" thing to do is to subtract the mean, although in this case
%-- it does not matter much because the mean is essentially zero.
rms_deviation_from_zero = sqrt(mean(timeseries.^2))
rms_deviation_from_mean = sqrt(mean((timeseries-timeseries_mean).^2))
rms_using_std_function = std(timeseries,1)     % The 1 here is meaningful --
                                               % do 'help std' to see why

input('Press RETURN to continue');

%-------------------------------
% Part 2

display('%----- Part 2a')
tsfft = fft(timeseries);

display('%----- Part 2b')
%-- For a complex number, 'abs' gives its magnitude
periodo = abs(tsfft).^2;
semilogy(periodo);
input('Whole periodogram; Press RETURN to continue');

%-- Now zoom in
xlim([1 1000]);
input('Low freq zoom, showing spectral leakage around 74.7 Hz but not around 820 Hz (because 820 Hz is an exactly periodic freq); Press RETURN to continue');

%-- Zoom in again
range = (16384+1-1000:16384);
semilogy(range,periodo(range));
input('High end of array, corresponding to negative freqs; Press RETURN to continue');

display('%----- Part 2c')
%-- Do 'help fft' to read about Matlab's convention for the frequency series
%-- that stores an fft.
help fft
input('Press RETURN to continue');
%-- The frequency runs from 0 to just under the sampling rate, NOT from
%-- negative Nyquist to positive Nyquist.  Mathematically, those cover
%-- equivalent frequency ranges, just with wrapping around.  Here are a
%-- few ways to generate the frequencies for the frequency series:
fpos = (0:ndata-1) / nsec;
display(sprintf('\nAll-positive frequency series runs from %g to %g Hz',...
       fpos(1), fpos(length(fpos)) ));
semilogy(fpos,periodo);
title('Part 2c: FFT vs. all-positive frequency series');
xlabel('Frequency (Hz)');
ylabel('Fourier transform of signal (arbitrary units)');
input('Periodogram vs. all-positive frequencies; Press RETURN to continue');

%-- Above the Nyquist frequency, subtract off the sampling frequency
toohigh = find(fpos>=nyquist);
f=fpos; f(toohigh) = f(toohigh) - samprate;
display(sprintf(...
'\nSymmetric frequency series runs from %g to %g Hz, then from %g to %g Hz',...
       fpos(1), max(f), min(f), f(length(f)) ));
semilogy(f,periodo);
title('Part 2c: FFT vs. symmetric frequency series');
xlabel('Frequency (Hz)');
ylabel('Fourier transform of signal (arbitrary units)');
input('Periodogram vs. symmetric frequencies; Press RETURN to continue');

%-- Zoom in
range = find(abs(f)<1000);
semilogy(f(range),periodo(range));
title('Part 2c: FFT vs. symmetric frequency series, zoomed');
xlabel('Frequency (Hz)');
ylabel('Fourier transform of signal (arbitrary units)');
input('Zoom of periodogram vs. symmetric frequencies; Press RETURN to continue');

display('%----- Part 2d')
%-- Remember the factor of ndata normalizing the frequency domain sum
time_domain_energy = sum( timeseries.^2 )
freq_domain_energy = sum( abs(tsfft).^2 ) / ndata
periodogram_sum = sum( periodo ) / ndata
%-- Those are equal!

display('%----- Part 2e')
%-- The Hann window is a "raised cosine" window
win = (1-cos(2*pi*t))/2;
windowed_timeseries = win.*timeseries;
wintsfft = fft(windowed_timeseries);
plot(t,windowed_timeseries);
title('Part 2e: Windowed time series');
xlabel('Time (s)');
ylabel('Signal (arbitrary units)');
input('Windowed data; Press RETURN to continue');

periodo = abs(wintsfft).^2;
semilogy(periodo);
input('Whole periodogram after windowing; Press RETURN to continue');

%-- Now zoom in
semilogy(1:1000,periodo(1:1000));
input('Low freq zoom, showing less spectral leakage (and much more similar) now; Press RETURN to continue');

%-------------------------------
% Part 3

display('%----- Part 3a')
%-- If you don't assign the output fom the pwelch function, it produces a nice
%-- formatted plot (though not with correct freq scale unless you tell it)
pwelch(timeseries)
input('pwelch produces a one-sided PSD with log y scale; Press RETURN to continue');

pwelch(timeseries,[],[],[],samprate)
input('pwelch with correct frequency scale (in kHz); Press RETURN to continue');

display('%----- Part 3b')
%-- Zoom in
set(gca,'XLim',[0 1.6]);
input('Less spectral leakage than part 2b, but coarser freq resolution due to averaging; Press RETURN to continue');

%-------------------------------
% Part 4

display('%----- Part 4a')
%-- Apply the simple difference filter
filtts = timeseries - [0 timeseries(1:length(timeseries)-1)];
%-- Plot it
plot(t,filtts);
title('Part 4a: Filtered time series')
%-- The filtered data looks more like white noise
input('Filtering suppressed the low-frequency content; Press RETURN to continue');

display('%----- Part 4b');
%-- Consider a sinusoidal input, x_k = A * exp(2*pi*i*f*k*Dt) where Dt is the
%-- time step.  It is not too hard to show that
%-- y_k = A * exp(2*pi*i*f*k*Dt) * 2 * sin(pi*f*Dt) * exp(-pi*i*f*Dt) .
%-- Thus the transfer function is y_k/x_k = 2 * sin(pi*f*Dt) * exp(-pi*i*f*Dt)
%-- i.e. magnitude sin(pi*f*Dt) and phase -pi*f*Dt

display('%----- Part 4c');
pwelch(filtts,[],[],[],samprate);
%-- The filter suppressed low frequencies and increased high frequencies some.
input('Low freq component suppressed, high freq increased some; Press RETURN to continue');

%-- Also plot amplitude spectral density on a linear scale
[pwdata,pwf] = pwelch(filtts,[],[],[],samprate);
plot(pwf,sqrt(pwdata));
title('Part 4c: Amplitude spectral density');
input('Previously flat noise floor now looks like a quarter sine wave; Press RETURN to continue');

display('%----- Part 4d');
%-- Design the filter, and apply it.  Note that you have to pass '2' as the
%-- first argument to get an 4th-order bandpass filter - a quirk of Matlab
[b,a] = butter(2,[65/nyquist 85/nyquist]);
%-- Filter the original time series
filtts = filter(b,a,timeseries);
%-- Plot it
plot(t,filtts);
title('Part 4d: Time series after Butterworth filter');
input('See mostly the 74.7 Hz component; Press RETURN to continue');

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

display('%----- Part 5a');
downts = timeseries(1:16:length(timeseries));
rms_original_timeseries = std(timeseries,1)
rms_every_16th_sample = std(downts,1)

display('%----- Part 5b');
%-- Design a 4th-order Butterworth lowpass filter, and apply it.
[b,a] = butter(4,512/nyquist,'low');
%-- Filter the original time series
filtts = filter(b,a,timeseries);
downts = filtts(1:16:length(filtts));
rms_anti_aliased = std(downts,1)
display('Lower rms, because the second sinusoidal peak and high-frequency white noise were filtered out');
