Frontmatter

The objective today is to get you up to speed to be able to handle the new homework assignment that is out.

Note about exercises

You are expected to follow along this notebook and recreate the same plots in MATLAB. Yes the code is presented in the textbook; but let's do it again by hand to understand what is happening.

That this notebook is not in MATLAB langauge shouldn't be a problem since this is an example straight out of the textbook and the code is already there for you.

Example 4.4 in CSSB

This example achieves two objectives:

• How to deal with unevenly sampled data

• How to generate the power spectral density

There are two files that we will use:

• hr_pre.mat: Heart rate observed in the normal state

• hr_med.mat: Heart rate observed in the meditative state

Look for these files in the Ed resources pane. The problem here is that the data was sampled unevenly. By this we mean that a pair of values were recorded whenever a heart beat was detected: (timestamp, heart_rate).

So we will need to manually convert the data to evenly sampled values.

Let us load these files up into a variable ws to simulate the MATLAB workspace

As usual when we start any sort of analysis, the first step is to plot the data and see what it looks like.

Let us confirm that these are indeed unevenly sampled. The diff command in MATLAB is used to take the numerical derivative. So if we have an evenly (also called linearly) spaced time vector, then its derivative should be constant (remember derivative of a linear function is a constant). So to check if your time vector is evenly/linearly sampled; it suffice to take the diff of it and plot it.

To generate a uniformly sampled sequence we will need two things

• A vector of evenly spaced time samples

• A method to interpolate the hr_med and hr_pre variables at those chosen time instances.

In MATLAB this is achieved using the interp1 command. Look up MATLAB documentation on how to use it.

Switching to even sampling

Convert MAT to Julia Array

To switch to even sampling, we will have to perform interpolation. This is done for us by MATLAB's interp1 command. This command takes as input the current (unevenly sampled) time vector, the values corersponding to this time vector and a new (evenly sampled) time vector for which we want inerpolated values. To create this evenly sampled "query" vector ... first we start by finding the minimum and maximum timestamps.

We want to sample between them at 100 Hz. This notebook by default only shows the very last output of a code block. So what you see below is the evenly sampled time vector for meditative state, nt_med.

Next we create the interpolated vectors and mean center them. We mean center them because the definition of the power spectral density as the Fourier transform of the autorcorrelation function already includes mean-centering (recall formula for auto-correlation). But we are going to do the direct computation method (from the signals Fourier transform); therefore we need to do it manually.

Calculating the PSD

Now that we have switched our dataset to be evenly sampeld values (via linear interpolation) let us calculate the Power Spectral Density using the direct method.

First we calculate the FFT and generate the frequency vector. The following MATLAB function is the equivalent of the built in fftfreq function used by this notebook.

fftfreq for MATLAB

function [freqs] = fftfreqs(n, fs)

	if mod(n, 2) == 0
    	freqs = [0:floor(n/2)-1, -floor(n/2):-1]  * fs/n;
	else
    	freqs = [0:floor((n-1)/2), floor(-(n-1)/2):-1]  * fs/n;
	end
end

The above function takes as input the length of the DFT to be performed and the frequency at which the signal whose DFT is being taken was sampled.

NOTE: The textbook creates the frequency vectors in MATLAB manually. See the example code in the textbook to see how it is done.

Averaging

Recall averaging reduces noise. The same is true even when constructing periodograms or spectrograms.

Read Lecture 09 if you haven't already as well as Example 4.5 in CSSB.

Example 4.6

In this example we are going to smoothen the above images by the use of averaging. Specifcally we use the welch method for averaging. MATLAB provides an implementation with its DSP or Signal Processing Toolbox and so is behind a paywall. Instead we can use the following replacement routine from the course website.

welch routine from book/notes

function [PS,f] = welch(x,nfft,noverlap,fs)
	%Function to calculate averaged spectrum
	%[ps,f] = welch(x,nfft,noverlap,fs);
	%  Output arguments
	%		sp	spectrogram
	%		f	frequency vector for plotting
	%  Input arguments
	%		x data
	%		nfft window size
	%      	noverlap number of overlaping points in adjacent segements
	%	    fs sample frequency
	%	Uses Hanning window
	%
	N = length(x);                     % Get data length
	half_segment = round(nfft/2);      % Half segment length
	if nargin < 4                      % Test inputs
		fs = 2*pi;                     % Default fs
	end
	if nargin < 3 || isempty(noverlap) == 1
		noverlap = half_segment;       % Set default overlap at 50%
	end
	nfft = round(nfft);                % Make sure nfft is an interger
	noverlap = round(noverlap);        % Make sure noverlap is an interger
	%
	f = (1:half_segment)* fs/(nfft);   % Calculate frequency vector
	increment = nfft - noverlap;       % Caluclate window shift
	nu_avgs =  round(N/increment);     % Determine the number of segments
	%
	for k = 1:nu_avgs 			       % Calculate spectra for each data point
	   first_point = 1 + (k - 1) * increment;
	   if (first_point + nfft -1) > N
		   first_point = N - nfft + 1;          % Prevent overflow
	   end
		data = x(first_point:first_point+nfft-1); % Get data segment
		if k == 1
		   PS = abs((fft(data)).^2);            % Calculate Power spectrum Eq. 4.16
	   else
		   PS = PS + abs((fft(data)).^2);       % Calculate Power spectrum
	   end
	end
	PS = PS(2:half_segment+1)/(nu_avgs*nfft/2);  % Remove redundant points, no avg
end

Check list I

Make sure your submission has all the above pairs of plots for full credit

Folding code