Merge branch 'Mark-Kramer/master'

Chapter5/Chapter 5 Solutions/README.md

@@ -1,3 +1,3 @@
 Chapter 5 Solutions
 
-Each solution is a function that is either self-contained or makes a call to various helper functions within this folder. However, please ensure that Chapter 5’s files (particularly the .mat files) are within this folder or within your MATLAB working directory. If not, then the various commands load('Ch3-EEG-#.mat') will fail. Question 8 also calls functions within the Chronic toolbox, which is freely distributed online. 
+Each solution is a function that is either self-contained or makes a call to various helper functions within this folder. However, please ensure that Chapter 5’s files (particularly the .mat files) are within this folder or within your MATLAB working directory. If not, then the various commands load('Ch5-EEG-#.mat') will fail. Question 8 also calls functions within the Chronic toolbox, which is freely distributed online. 

Chapter6/Chapter 6 Solutions/Chapter_6_Question_1.m

@@ -0,0 +1,104 @@
+function Chapter_6_Question_1()
+ % Question 1 - Analyze Ch6-EEG-2.mat
+
+ %Load Data
+ load('Ch6-EEG-2.mat')
+ x = EEG - mean(EEG);
+ T = t(end);
+ N = length(t);
+ sample_interval = t(2) - t(1);
+ sample_freq = 1/sample_interval;
+ nyquist_freq = sample_freq/2;
+ freq_resolution = 1/T;
+ faxis = 0:freq_resolution:nyquist_freq;
+
+ %Visualize Data
+ figure()
+ plot(t, x, 'k', 'LineWidth', 2)
+ xlabel('Time (seconds)')
+ ylabel('Data')
+ title('Single Trial EEG Data')
+ set(gca, 'FontSize', 14)
+
+ %Get Fourier Transform of Data
+ xf = fft(x);
+ Sxx = 2*sample_interval^2/T * (xf.*conj(xf));
+ Sxx = Sxx(1:N/2+1); %ignore negative freqs
+ Sxx(1) = 0; %set DC freq to zero power to correct for MATLAB rounding
+
+ %Plot Fourier Transform of Data
+ figure()
+ plot(faxis, 10*log10(Sxx), 'r', 'LineWidth', 2)
+ xlim([0 80])
+ xlabel('Frequency (Hz)')
+ ylabel('Power (dB)')
+ title('Fourier Transform Of EEG Data')
+ set(gca, 'FontSize', 14)
+
+ %Construct a Low Pass Filter for +/- 20 Hz.
+ n = 100;
+ Wn = 20/nyquist_freq;
+ b = transpose(fir1(n, Wn, 'low'));
+
+ %Fourier Transform Filter
+ bf = fft(b, N);
+ Mb = bf.*conj(bf);
+ faxis_filter = fftshift((-nyquist_freq:freq_resolution:nyquist_freq-freq_resolution));
+ [~, isorted] = sort(faxis_filter, 'descend');
+
+ %Visualize FFT of FIR Filter
+ figure()
+ plot(faxis_filter(isorted), Mb(isorted), 'b', 'LineWidth', 2)
+ xlabel('Frequency (Hz)')
+ ylabel('Magnitude')
+ title('FFT of FIR Filter')
+ xlim([-30, 30])
+ set(gca, 'FontSize', 14)
+
+ %Apply FIR Filter
+ xnew = filtfilt(b, 1, x);
+
+ %Get Spectrum of New Data
+ xfnew = fft(xnew);
+ Sxxnew = 2*sample_interval^2/T * (xfnew.*conj(xfnew));
+ Sxxnew = Sxxnew(1:N/2+1); %ignore negative frequencies
+ Sxxnew(1) = 0; %set DC freq to zero power to correct for MATLAB rounding
+
+ %Visualize New Data
+ figure()
+ subplot(2, 1, 1)
+ plot(t, xnew, 'k', 'LineWidth', 2)
+ xlabel('Time (seconds)')
+ ylabel('Data')
+ title('Filtered Data')
+ set(gca, 'FontSize', 14)
+ subplot(2, 1, 2)
+ plot(faxis, 10*log10(Sxxnew), 'r', 'LineWidth', 2)
+ xlabel('Frequency (Hz)')
+ ylabel('Power (dB)')
+ xlim([0 80])
+ title('Filtered EEG Data FFT')
+ set(gca, 'FontSize', 14)
+
+ % a) Unfortunately the data appears to be a bit too noisy for me to make
+ % out any structure in the data. Thus, I can't observe any clear rhythms. 
+
+ % b) The fourier transform shows a lot of flucutating noise/power in the
+ % frequencies greater than 20 Hz, but there is no clear peak. While there
+ % is a peak at 10 Hz, this rhythm is not visible in the EEG data. Thus, it
+ % seems like the EEG contains a plethora of frequencies above 20 or 30 Hz,
+ % but none seem to jump out, indicating that they might all just be white
+ % noise. That lone 10 Hz peak is what we will focus on, since we are
+ % assuming nothing neural would create power in all frequencies above 20 to
+ % 30 Hz in such a way as seen here.
+
+ % c) We construct a low pass filter to isolate the frequencies +/- 20 Hz.
+ % This will eliminate the higher frequency (seemingly) white noise and
+ % leave us with our 10 Hz sinusoid. After creating the filter, we visualize
+ % it in the frequency domain to ensure it's doing what we expect. Then, we
+ % use built-in MATLAB functions to apply the filter and create filtered EEG
+ % data. Then we Fourier Transform the filtered data to clearly see the 10
+ % Hz sinusoid that makes up our signal. The sinusoid is also clearly
+ % visible in the filtered EEG.
+
+end

Chapter6/Chapter 6 Solutions/Chapter_6_Question_2.m

@@ -0,0 +1,84 @@
+function Chapter_6_Question_2()
+ % Question 2 - Explore Hanning Filters
+
+ %Load Data
+ load('Ch6-EEG-1.mat')
+ x = EEG(1, :) - mean(EEG(1, :));
+ T = t(end);
+ N = length(t);
+ sample_interval = t(2) - t(1);
+ sample_freq = 1/sample_interval;
+ nyquist_freq = sample_freq/2;
+ freq_resolution = 1/T;
+ faxis = fftshift(-nyquist_freq:freq_resolution:nyquist_freq-freq_resolution);
+
+ %Find Indices For +/- 60 Hz
+ ind = find((abs(abs(single(faxis))-60)) <= 0.1*freq_resolution);
+
+ %Describe Three Hann Filters (in Frequency Domain)
+ win1 = 7;
+ win2 = 15;
+ win3 = 30;
+ hann_filter1 = ones(1, N);
+ hann_filter1(ind(1)-win1:ind(1)+win1) = 1-transpose(hann(2*win1+1));
+ hann_filter1(ind(2)-win1:ind(2)+win1) = 1-transpose(hann(2*win1+1));
+ hann_filter2 = ones(1, N);
+ hann_filter2(ind(1)-win2:ind(1)+win2) = 1-transpose(hann(2*win2+1));
+ hann_filter2(ind(2)-win2:ind(2)+win2) = 1-transpose(hann(2*win2+1));
+ hann_filter3 = ones(1, N);
+ hann_filter3(ind(1)-win3:ind(1)+win3) = 1-transpose(hann(2*win3+1));
+ hann_filter3(ind(2)-win3:ind(2)+win3) = 1-transpose(hann(2*win3+1));
+
+ %Get Time Domain Hann Filters
+ lag_axis = (-N/2+1:N/2)*sample_interval;
+ impulse = zeros(1, N);
+ impulse(N/2) = 1;
+ i_hann_filter1 = ifft(hann_filter1);
+ i_hann_filter2 = ifft(hann_filter2);
+ i_hann_filter3 = ifft(hann_filter3);
+ impulse_response_t1 = zeros(N, 1);
+ impulse_response_t2 = zeros(N, 1);
+ impulse_response_t3 = zeros(N, 1);
+ for i = 1:N
+ impulse_response_t1(i) = sum(circshift(i_hann_filter1, [1, i-1]).*impulse);
+ impulse_response_t2(i) = sum(circshift(i_hann_filter2, [1, i-1]).*impulse);
+ impulse_response_t3(i) = sum(circshift(i_hann_filter3, [1, i-1]).*impulse);
+ end
+
+ %Visualize Frequency and Time Domain of Filters
+ figure()
+ subplot(2, 1, 1)
+ hold on
+ plot(faxis, hann_filter1, 'b', 'LineWidth', 2)
+ plot(faxis, hann_filter2, 'g', 'LineWidth', 2)
+ plot(faxis, hann_filter3, 'r', 'LineWidth', 2)
+ hold off
+ xlabel('Frequency (Hz)')
+ ylabel('Magnitude')
+ xlim([-80 80])
+ ylim([0 1])
+ legend({'Win = 7', 'Win = 15', 'Win = 30'})
+ title('Fourier Representation of Hanning Filter')
+ set(gca, 'FontSize', 14)
+ subplot(2, 1, 2)
+ hold on
+ plot(lag_axis, impulse_response_t1, 'b', 'LineWidth', 2)
+ plot(lag_axis, impulse_response_t2, 'g', 'LineWidth', 2)
+ plot(lag_axis, impulse_response_t3, 'r', 'LineWidth', 2)
+ hold off
+ xlabel('Time (Seconds)')
+ ylabel('Impulse Response')
+ xlim([-0.5 0.5])
+ ylim([-.03 .03])
+ legend({'Win = 7', 'Win = 15', 'Win = 30'})
+ title('Time Representation of Hanning Filter')
+ set(gca, 'FontSize', 14)
+
+ % In the frequency domain, increasing the width of the window causes more
+ % frequencies to be affected by the filter. More specifically, the filter
+ % more gradually changes from 0 to 1 and vice versa. In the time domain,
+ % this causes an increase in the precision with which the filter works.
+ % Larger windows are tightly centered around the impulse, whereas small
+ % windows extend over more time.
+
+end

Chapter6/Chapter 6 Solutions/Chapter_6_Question_3.m

@@ -0,0 +1,71 @@
+function Chapter_6_Question_3()
+ % Question 3 - Trial 1 Autocovariance Before/After Filters
+
+ %Load Data
+ load('Ch6-EEG-1.mat')
+ x = EEG(1, :) - mean(EEG(1, :));
+ T = t(end);
+ N = length(t);
+ sample_interval = t(2) - t(1);
+ sample_freq = 1/sample_interval;
+ nyquist_freq = sample_freq/2;
+ freq_resolution = 1/T;
+ faxis = fftshift(-nyquist_freq:freq_resolution:nyquist_freq-freq_resolution);
+
+ %Find Frequencies Near 60 Hz
+ indices = find((abs(abs(single(faxis))-60)) <= 1);
+
+ %Create And Apply Rectangular Filter
+ rectangular_filter = ones(1, N);
+ rectangular_filter(indices) = 0;
+ xf = fft(x);
+ xfnew = xf.*rectangular_filter;
+ x_rect = ifft(xfnew);
+
+ %Create And Apply FIR Filter
+ n = 100;
+ Wn = 30/nyquist_freq;
+ b = transpose(fir1(n, Wn, 'low'));
+ x_FIR = filtfilt(b, 1, x);
+
+ %Visualize Data Before/After Filter(s)
+ figure()
+ hold on
+ plot(t, x, 'k', 'LineWidth', 2)
+ plot(t, x_rect, 'b', 'LineWidth', 2)
+ plot(t, x_FIR, 'r', 'LineWidth', 2)
+ hold off
+ xlabel('Time (seconds)')
+ ylabel('Data')
+ legend({'Raw Data', 'Data + Rect Filter', 'Data + FIR Filter'})
+ title('Data Before/After Filter(s)')
+ set(gca, 'FontSize', 14)
+
+ %Get Biased Autocovariance
+ [x_autocov, lags] = xcorr(x, 200, 'biased');
+ [x_rect_autocov, ~] = xcorr(x_rect, 200, 'biased');
+ [x_FIR_autocov, ~] = xcorr(x_FIR, 200, 'biased');
+
+ %Plot Biased Autocovariance
+ figure()
+ hold on
+ plot(lags.*sample_interval, x_autocov, 'k', 'LineWidth', 2)
+ plot(lags.*sample_interval, x_rect_autocov, 'b', 'LineWidth', 2)
+ plot(lags.*sample_interval, x_FIR_autocov, 'r', 'LineWidth', 2)
+ hold off
+ xlabel('Time (seconds)')
+ ylabel('Autocovariance')
+ legend({'Raw Data', 'Data + Rect Filter', 'Data + FIR Filter'})
+ title('Autocovariance Before/After Filter(s)')
+ set(gca, 'FontSize', 14)
+
+ % The raw data show obvious 60 Hz correlation, which is the noise we are
+ % trying to get rid of. The rectangular-filtered data has massively reduced
+ % this, but if you zoom in on the y-axis you can see minor fluctuations
+ % that create a 60 Hz rhythm. However, it does pick up the roughly 20 Hz
+ % rhythm that underlies the evoked response, since we see a few undulations
+ % that occur once every 0.05 seconds or so (20 Hz). The FIR-filtered data
+ % shows the same .05 undulations, but without the added small 60 Hz
+ % fluctuations.
+
+end

Chapter6/Chapter 6 Solutions/Chapter_6_Question_4.m

@@ -0,0 +1,73 @@
+function Chapter_6_Question_4()
+ % Question 4 - Changing Order of FIR Filters
+
+ %Load the EEG data to define useful parameters.
+ load('Ch6-EEG-1.mat') 
+ x = EEG(1,:) - mean(EEG(1,:));
+ sample_interval = t(2) - t(1);
+ N = length(x);
+ sample_freq = 1/sample_interval;
+ nyquist_freq = sample_freq/2;
+ freq_resolution = 1/(N*sample_interval);
+ faxis = fftshift((-nyquist_freq:freq_resolution:nyquist_freq-freq_resolution));
+
+ %Define the filter orders, cutoff frequency, and colors.
+ n = [1000, 500, 100, 50];
+ Wn = 30/nyquist_freq; 
+ clr = {'r', 'g', 'b', 'm'};
+ filtered_data = zeros(length(n), N);
+
+ %For each filter, design and visualize the filter.
+ for k = 1:length(n)
+
+ %Design filter.
+ b = fir1(n(k), Wn, 'low');
+
+ %Visualize filter.
+ subplot(2,1,1)
+ hold on
+ plot(b+0.05*k, 'Color', clr{k})
+ hold off
+ axis tight
+ axis off
+
+ bf = fft(b,N); %NOTE: zero pad to same length as x.
+ subplot(2,1,2)
+ hold on
+ plot(faxis,bf.*conj(bf), 'Color', clr{k})
+ hold off
+ axis tight
+ xlim([-55 55])
+ xlabel('Frequency [Hz]')
+
+ %Store Filtered Data (note: we cannot use filtfilt)
+ x1 = filter(b, 1, x); %filter once
+ x1 = x1(N:-1:1); %reverse
+ x2 = filter(b, 1, x1); %filter again
+ filtered_data(k, :) = x2(N:-1:1); %reverse again
+
+ end
+ subplot(2,1,1)
+ legend({'n=1000'; 'n=500'; 'n=100'; 'n=50'})
+
+ %Plot Filtered Data
+ figure()
+ subplot(5, 1, 1)
+ plot(t, x, 'k', 'LineWidth', 2)
+ xlabel('Time (seconds)')
+ ylabel('Data')
+ title('Raw Data')
+ set(gca, 'FontSize', 14)
+ for i = 2:5
+ subplot(5, 1, i)
+ plot(t, filtered_data(i-1, :), clr{i-1}, 'LineWidth', 2)
+ xlabel('Time (seconds)')
+ ylabel('Data')
+ title(['Order: ', num2str(n(i-1))])
+ set(gca, 'FontSize', 14)
+ end
+
+ % Increasing the order of the FIR filter causes the EEG signal to smooth
+ % out. 
+
+end

Chapter6/Chapter 6 Solutions/Chapter_6_Question_5.m

@@ -0,0 +1,51 @@
+function Chapter_6_Question_5()
+ % Question 5 - Bandstop FIR Filter
+
+ %Load the EEG data to define useful parameters.
+ load('Ch6-EEG-1.mat') 
+ x = EEG(1,:) - mean(EEG(1,:));
+ sample_interval = t(2) - t(1);
+ N = length(x);
+ sample_freq = 1/sample_interval;
+ nyquist_freq = sample_freq/2;
+ freq_resolution = 1/(N*sample_interval);
+ faxis = fftshift((-nyquist_freq:freq_resolution:nyquist_freq-freq_resolution));
+ [~, isorted] = sort(faxis, 'descend');
+
+ %Create Bandstop Filter (eliminate 50 to 70 Hz)
+ n = 100;
+ Wn = [50/nyquist_freq, 70/nyquist_freq]; 
+ b = transpose(fir1(n, Wn, 'stop'));
+
+ %Fourier Transform Bandstop Filter
+ bf = fft(b, N);
+ Mb = bf.*conj(bf);
+
+ %Visualize Fourier Transform of Bandstop Filter
+ figure()
+ plot(faxis(isorted), Mb(isorted), 'b', 'LineWidth', 2)
+ xlabel('Frequency (Hz)')
+ ylabel('Magnitude')
+ title('FFT of Bandstop FIR Filter')
+ xlim([-100, 100])
+ set(gca, 'FontSize', 14)
+
+ %Filter Data
+ xnew = filtfilt(b, 1, x);
+
+ %Visualize Data Before/After Filter
+ figure()
+ subplot(2, 1, 1)
+ plot(t, x, 'k', 'LineWidth', 2)
+ xlabel('Time (seconds)')
+ ylabel('Data')
+ title('Raw Data')
+ set(gca, 'FontSize', 14)
+ subplot(2, 1, 2)
+ plot(t, xnew, 'b', 'LineWidth', 2)
+ xlabel('Time (seconds)')
+ ylabel('Data')
+ title('Data after Bandstop Filter')
+ set(gca, 'FontSize', 14)
+
+end