diff --git a/HW_1_PNS_EE379K-385V_Neural_Eng/HW_1_PNS_sp26_EE379K_385J_Neural_Eng.pdf b/HW_1_PNS_EE379K-385V_Neural_Eng/HW_1_PNS_sp26_EE379K_385J_Neural_Eng.pdf new file mode 100644 index 0000000..ca27108 Binary files /dev/null and b/HW_1_PNS_EE379K-385V_Neural_Eng/HW_1_PNS_sp26_EE379K_385J_Neural_Eng.pdf differ diff --git a/HW_1_PNS_EE379K-385V_Neural_Eng/c1_dataVis.m b/HW_1_PNS_EE379K-385V_Neural_Eng/c1_dataVis.m index 604a4aa..900e552 100644 --- a/HW_1_PNS_EE379K-385V_Neural_Eng/c1_dataVis.m +++ b/HW_1_PNS_EE379K-385V_Neural_Eng/c1_dataVis.m @@ -7,51 +7,192 @@ clc load('data.mat'); %% Example: Plot the raw signal -signal=Flex.signal; -labels=Flex.trigger; -TRIG = gettrigger(labels,0.5); -TRIGend = gettrigger(-labels,-0.5); +flexSignal=Flex.signal; +pinchSignal = Pinch.signal; +VFSignal = VF.signal; +flexLabels=Flex.trigger; +pinchLabels=Pinch.trigger; +VFLabels = VF.trigger; -figure('units','normalized','Position',[0.1,0.1,0.7,0.4]) -plot((1:length(signal))./fs,zscore(signal)); -hold on; -plot((1:length(signal))./fs,zscore(labels),'y'); -stem(TRIG./fs,ones(length(TRIG),1)*max(zscore(labels)),'Color','g'); -stem(TRIGend./fs,ones(length(TRIG),1)*max(zscore(labels)),'Color','r'); -grid on; grid minor; -xlim([0,length(signal)./fs]) -xlabel('Time (s)') -ylabel('Amplitude (uV)') -title('Raw VF signal with labels for stimulation periods') +%% Plot the raw signals +% Grouping variables into cell arrays to loop through them cleanly +signals = {flexSignal, pinchSignal, VFSignal}; +labels_all = {flexLabels, pinchLabels, VFLabels}; +titles = {'Raw Flex Signal with Stimulation Pattern (Yellow)', ... + 'Raw Pinch Signal with Stimulation Pattern (Yellow)', ... + 'Raw VF Signal with Stimulation Pattern (Yellow)'}; + +% Create one large figure for all three subplots +figure('units','normalized','Position',[0.1, 0.1, 0.7, 0.8]) + +for i = 1:3 + sig = signals{i}; + lbls = labels_all{i}; + + % Find trigger start and end points using your custom function + TRIG = gettrigger(lbls, 0.5); + TRIGend = gettrigger(-lbls, -0.5); + + % Create a subplot (3 rows, 1 column, current index i) + subplot(3, 1, i); + + % Plot normalized signal and labels + plot((1:length(sig))./fs, zscore(sig)); + hold on; + plot((1:length(sig))./fs, zscore(lbls), 'y', 'LineWidth', 1.5); + + % Plot stem markers for triggers (added conditional checks just in case a signal has no triggers) + if ~isempty(TRIG) + stem(TRIG./fs, ones(length(TRIG),1)*max(zscore(lbls)), 'Color', 'g'); + end + if ~isempty(TRIGend) + stem(TRIGend./fs, ones(length(TRIGend),1)*max(zscore(lbls)), 'Color', 'r'); + end + + % Formatting + grid on; grid minor; + xlim([0, length(sig)./fs]); + xlabel('Time (s)'); + % Note: Because you used zscore(), the amplitude is no longer in uV, but in standard deviations + ylabel('Amplitude (stdDevs)'); + title(titles{i}); +end %% Example: PSD estimates -figure('units','normalized','Position',[0.1,0.1,0.5,0.5]) -[rows_act,cols_act,values_act] = find(labels>0); -[rows_rest1,cols_rest,values_rest] = find(labels==0); -notOfInterest = signal(rows_rest1); -signalOfInterest=signal(rows_act); -h = spectrum.welch; % creates the Welch spectrum estimator -SOIf=psd(h,signalOfInterest,'Fs',fs); % calculates and plot the one sided PSD -plot(SOIf); % Plot the one-sided PSD. -temp =get(gca); -temp.Children(1).Color = 'b'; +figure('Name', 'Raw PSD Estimates', 'units','normalized','Position',[0.1,0.1,0.5,0.5]) +h = spectrum.welch; -% %% Bandpass Filtering +% --- REST (All signals combined) --- +% Find indices where labels are 0 for each signal type +[flex_rows_rest, ~, ~] = find(flexLabels == 0); +[pinch_rows_rest, ~, ~] = find(pinchLabels == 0); +[vf_rows_rest, ~, ~] = find(VFLabels == 0); + +% Extract the actual signal data during those rest periods +flex_restData = flexSignal(flex_rows_rest); +pinch_restData = pinchSignal(pinch_rows_rest); +vf_restData = VFSignal(vf_rows_rest); + +% Concatenate all rest data into one large vector for a robust estimate +all_restData = [flex_restData; pinch_restData; vf_restData]; + +% Calculate and Plot Rest PSD in Black +SOIf_rest = psd(h, all_restData, 'Fs', fs); +plot(SOIf_rest.Frequencies, 10*log10(SOIf_rest.Data), 'k', 'LineWidth', 1.5); +hold on; + +% --- FLEX --- +[flex_rows_act, ~, ~] = find(flexLabels>0); +flex_signalOfInterest = flexSignal(flex_rows_act); +SOIf_flex = psd(h, flex_signalOfInterest, 'Fs', fs); +plot(SOIf_flex.Frequencies, 10*log10(SOIf_flex.Data), 'g'); % Plot in Green + +% --- PINCH --- +[pinch_rows_act, ~, ~] = find(pinchLabels>0); +pinch_signalOfInterest = pinchSignal(pinch_rows_act); +SOIf_pinch = psd(h, pinch_signalOfInterest, 'Fs', fs); +plot(SOIf_pinch.Frequencies, 10*log10(SOIf_pinch.Data), 'r'); % Plot in Red + +% --- VF --- +[VF_rows_act, ~, ~] = find(VFLabels>0); +VF_signalOfInterest = VFSignal(VF_rows_act); +SOIf_VF = psd(h, VF_signalOfInterest, 'Fs', fs); +plot(SOIf_VF.Frequencies, 10*log10(SOIf_VF.Data), 'b'); % Plot in Blue + +% --- FORMATTING --- +grid on; +xlabel('Frequency (Hz)'); +ylabel('Power/Frequency (dB/Hz)'); +legend('Rest', 'Flex', 'Pinch', 'VF'); +title('Power Spectral Density Estimates (Raw Signals)'); + +%% Bandpass Filtering % -% fc1 = 0; % first cutoff frequency in Hz -% fc2 = 100; % second cutoff frequency in Hz +fc1 = 800; % first cutoff frequency in Hz +fc2 = 2200; % second cutoff frequency in Hz % % % normalize the frequencies -% Wp = [fc1 fc2]*2/fs; +Wp = [fc1 fc2]*2/fs; % Build a Butterworth bandpass filter of 4th order % check the "butter" function in matlab +n = 2; +[b, a] = butter(n, Wp, 'bandpass'); % Filter data of both classes with a non-causal filter % Hint: use "filtfilt" function in MATLAB % filteredSignal = ; -%% +filteredFlex = filtfilt(b, a, flexSignal); +filteredPinch = filtfilt(b, a, pinchSignal); +filteredVF = filtfilt(b, a, VFSignal); + +%% Compare VF Signal Before and After Filtering (Time Domain) + +% Group the raw and filtered VF signals to loop through them cleanly +vf_signals = {VFSignal, filteredVF}; +vf_labels = {VFLabels, VFLabels}; % Labels remain the exact same +vf_titles = {'Raw VF Signal with Stimulation Pattern (Yellow)', ... + 'Filtered VF Signal (800-2200 Hz) with Stimulation Pattern (Yellow)'}; + +% Create one figure for the two subplots (slightly shorter since it's only 2 plots) +figure('Name', 'VF Filter Comparison (Time Domain)', 'units', 'normalized', 'Position', [0.1, 0.1, 0.7, 0.6]) + +for i = 1:2 + sig = vf_signals{i}; + lbls = vf_labels{i}; + + % Find trigger start and end points using your custom function + TRIG = gettrigger(lbls, 0.5); + TRIGend = gettrigger(-lbls, -0.5); + + % Create a subplot (2 rows, 1 column, current index i) + subplot(2, 1, i); + + % Plot normalized signal and labels + plot((1:length(sig))./fs, zscore(sig)); + hold on; + plot((1:length(sig))./fs, zscore(lbls), 'y', 'LineWidth', 1.5); + + % Plot stem markers for triggers + if ~isempty(TRIG) + stem(TRIG./fs, ones(length(TRIG),1)*max(zscore(lbls)), 'Color', 'g'); + end + if ~isempty(TRIGend) + stem(TRIGend./fs, ones(length(TRIGend),1)*max(zscore(lbls)), 'Color', 'r'); + end + + % Formatting + grid on; grid minor; + xlim([0, length(sig)./fs]); + xlabel('Time (s)'); + ylabel('Amplitude (stdDevs)'); + title(vf_titles{i}); +end +%% Compare VF Signal Before and After Filtering (Frequency Domain / PSD) +figure('Name', 'VF PSD Comparison', 'units', 'normalized', 'Position', [0.2, 0.2, 0.5, 0.4]) + +% --- RAW VF --- +[VF_rows_act, ~, ~] = find(VFLabels>0); +VF_raw_signalOfInterest = VFSignal(VF_rows_act); +h = spectrum.welch; +SOIf_VF_raw = psd(h, VF_raw_signalOfInterest, 'Fs', fs); +plot(SOIf_VF_raw.Frequencies, 10*log10(SOIf_VF_raw.Data), 'b'); % Plot Raw in Blue +hold on; + +% --- FILTERED VF --- +% Reuse the exact same trigger rows since the timing hasn't changed +VF_filt_signalOfInterest = filteredVF(VF_rows_act); +SOIf_VF_filt = psd(h, VF_filt_signalOfInterest, 'Fs', fs); +plot(SOIf_VF_filt.Frequencies, 10*log10(SOIf_VF_filt.Data), 'r'); % Plot Filtered in Red + +% --- FORMATTING --- +grid on; grid minor; +title('PSD of VF Signal Before and After Filtering'); +xlabel('Frequency (Hz)'); +ylabel('Power/Frequency (dB/Hz)'); +legend('Raw VF', 'Filtered VF (800 - 2200 Hz)'); +% xlim([0, 3000]); % Zoom in to see the filter cutoff roll-off clearly \ No newline at end of file diff --git a/HW_1_PNS_EE379K-385V_Neural_Eng/c2_discriminability.m b/HW_1_PNS_EE379K-385V_Neural_Eng/c2_discriminability.m new file mode 100644 index 0000000..83cbb11 --- /dev/null +++ b/HW_1_PNS_EE379K-385V_Neural_Eng/c2_discriminability.m @@ -0,0 +1,114 @@ +%% Setup and Feature Extraction +% Ensure c1_dataVis.m has been run so filtered signals and labels exist. + +% Parameters from Part 2.2b (Feature Selection) +WSize_sec = 0.1; +Olap_pct = 0; + +% Group data for easy looping +signals_list = {filteredFlex, filteredPinch, filteredVF}; +labels_list = {flexLabels, pinchLabels, VFLabels}; +names = {'Flex', 'Pinch', 'VF'}; + +% Initialize struct to store the extracted features for each class +results = struct(); + +fprintf('Extracting features (WSize=%.2fs, Olap=%.2f)...\n', WSize_sec, Olap_pct); + +for k = 1:3 + sig = signals_list{k}; + lbl = labels_list{k}; + + % --- Windowing Logic --- + WSize_samp = floor(WSize_sec * fs); + nOlap = floor(Olap_pct * WSize_samp); + hop = WSize_samp - nOlap; + nx = length(sig); + len = fix((nx - (WSize_samp - hop)) / hop); + + % Preallocate + MAV_vec = zeros(1, len); + VAR_vec = zeros(1, len); + feat_lbl = zeros(1, len); + + % Get triggers for labeling + Rise = gettrigger(lbl, 0.5); + Fall = gettrigger(-lbl, -0.5); + + for i = 1:len + idx_start = (i-1)*hop + 1; + idx_end = idx_start + WSize_samp - 1; + + segment = sig(idx_start:idx_end); + + % Calculate Features + MAV_vec(i) = mean(abs(segment)); + VAR_vec(i) = var(segment); + + % Labeling: 1 if window is strictly inside stimulation, 0 otherwise + is_stim = any(idx_start >= Rise & idx_end <= Fall); + feat_lbl(i) = double(is_stim); + end + + % --- Separate Rest vs Stimulus Data --- + % Store the features in the results struct + results(k).Name = names{k}; + results(k).MAV_Stim = MAV_vec(feat_lbl == 1); + results(k).MAV_Rest = MAV_vec(feat_lbl == 0); + results(k).VAR_Stim = VAR_vec(feat_lbl == 1); + results(k).VAR_Rest = VAR_vec(feat_lbl == 0); +end + +%% Calculate Fisher's Discriminant Ratio (FDR) +% Formula: J = (mu1 - mu2)^2 / (var1^2 + var2^2) +% Note: Using Variance (sigma^2) directly in denominator + +fprintf('\n=== Fisher''s Discriminant Ratio (FDR) Results ===\n'); +fprintf('Higher value = Better discriminability\n'); +fprintf('----------------------------------------------------------\n'); +fprintf('%-25s | %-12s | %-12s\n', 'Comparison', 'FDR (MAV)', 'FDR (VAR)'); +fprintf('----------------------------------------------------------\n'); + +% 1. Rest vs Stimulus (Internal comparison for each file) +for k = 1:3 + % Data for Stimulus vs Rest + data_stim = results(k); + + % --- MAV --- + mu1 = mean(data_stim.MAV_Stim); var1 = var(data_stim.MAV_Stim); + mu2 = mean(data_stim.MAV_Rest); var2 = var(data_stim.MAV_Rest); + fdr_mav = ((mu1 - mu2)^2) / (var1 + var2); + + % --- VAR --- + mu1 = mean(data_stim.VAR_Stim); var1 = var(data_stim.VAR_Stim); + mu2 = mean(data_stim.VAR_Rest); var2 = var(data_stim.VAR_Rest); + fdr_var = ((mu1 - mu2)^2) / (var1 + var2); + + fprintf('%-25s | %-12.4f | %-12.4f\n', [names{k} ' vs Rest'], fdr_mav, fdr_var); +end + +fprintf('----------------------------------------------------------\n'); + +% 2. Stimulus vs Stimulus (Compare 'Stim' vector of one to 'Stim' vector of another) +pairs = [1 2; 1 3; 2 3]; % Flex-Pinch, Flex-VF, Pinch-VF + +for p = 1:3 + idx1 = pairs(p, 1); + idx2 = pairs(p, 2); + + name1 = names{idx1}; + name2 = names{idx2}; + + % --- MAV --- + mu1 = mean(results(idx1).MAV_Stim); var1 = var(results(idx1).MAV_Stim); + mu2 = mean(results(idx2).MAV_Stim); var2 = var(results(idx2).MAV_Stim); + fdr_mav = ((mu1 - mu2)^2) / (var1 + var2); + + % --- VAR --- + mu1 = mean(results(idx1).VAR_Stim); var1 = var(results(idx1).VAR_Stim); + mu2 = mean(results(idx2).VAR_Stim); var2 = var(results(idx2).VAR_Stim); + fdr_var = ((mu1 - mu2)^2) / (var1 + var2); + + fprintf('%-25s | %-12.4f | %-12.4f\n', [name1 ' vs ' name2], fdr_mav, fdr_var); +end +fprintf('----------------------------------------------------------\n'); \ No newline at end of file diff --git a/HW_1_PNS_EE379K-385V_Neural_Eng/c2_featureExtraction.asv b/HW_1_PNS_EE379K-385V_Neural_Eng/c2_featureExtraction.asv new file mode 100644 index 0000000..8f142d4 --- /dev/null +++ b/HW_1_PNS_EE379K-385V_Neural_Eng/c2_featureExtraction.asv @@ -0,0 +1,35 @@ + +filteredSignal = VF_filt_signalOfInterest; % bandapass filtered signal +label = VFLabels; % labels of stimulus locations + +WSize = 50; % window size in s +Olap = 0; % overlap percentage + +%% Extracting Features over overlapping windows + +WSize = floor(WSize*fs); % length of each data frame, 30ms +nOlap = floor(Olap*WSize); % overlap of successive frames, half of WSize +hop = WSize-nOlap; % amount to advance for next data frame +nx = length(signal); % length of input vector +len = fix((nx - (WSize-hop))/hop); %length of output vector = total frames + +% preallocate outputs for speed +[MAV_feature, VAR_feature, featureLabels] = deal(zeros(1,len)); + +Rise1 = gettrigger(label,0.5); % gets the starting points of stimulations +Fall1 = gettrigger(-label,-0.5); % gets the ending points of stimulations + +for i = 1:len + segment = filteredSignal(((i-1)*hop+1):((i-1)*hop+WSize)); + MAV_feature(i) = ; + VAR_feature(i) = ; + + % re-build the label vector to match it with the feature vector + featureLabels(i) = sum(arrayfun(@(t) ((i-1)*hop+1) >= Rise1(t) && ((i-1)*hop+WSize) <= Fall1(t), 1:length(Rise1))); +end + +%% Plotting the features +% Note: when plotting the features, scale the featureLabels to the max of +% the feature values for proper visualization + + diff --git a/HW_1_PNS_EE379K-385V_Neural_Eng/c2_featureExtraction.m b/HW_1_PNS_EE379K-385V_Neural_Eng/c2_featureExtraction.m index b960386..a9cbbce 100644 --- a/HW_1_PNS_EE379K-385V_Neural_Eng/c2_featureExtraction.m +++ b/HW_1_PNS_EE379K-385V_Neural_Eng/c2_featureExtraction.m @@ -1,35 +1,123 @@ +%% Parameter Sweep Setup +% Note: Ensure you have run c1_dataVis.m first to get filteredFlex/fs/flexLabels +filteredSignal = filteredPinch; % Using Flex signal (High SNR) for demonstration +label = pinchLabels; % Labels of stimulus locations -filteredSignal = ; % bandapass filtered signal -label = ; % labels of stimulus locations +% Sweep parameters +WSize_values = [0.05, 0.1, 0.3]; % Window sizes in seconds +Olap_values = [0, 0.25, 0.75]; % Overlap percentages -WSize = ; % window size in s -Olap = ; % overlap percentage +% Get trigger points once (indices in the raw signal) +Rise1 = gettrigger(label, 0.5); % Start of stimulation +Fall1 = gettrigger(-label, -0.5); % End of stimulation -%% Extracting Features over overlapping windows - -WSize = floor(WSize*fs); % length of each data frame, 30ms -nOlap = floor(Olap*WSize); % overlap of successive frames, half of WSize -hop = WSize-nOlap; % amount to advance for next data frame -nx = length(signal); % length of input vector -len = fix((nx - (WSize-hop))/hop); %length of output vector = total frames - -% preallocate outputs for speed -[MAV_feature, VAR_feature, featureLabels] = deal(zeros(1,len)); - -Rise1 = gettrigger(label,0.5); % gets the starting points of stimulations -Fall1 = gettrigger(-label,-0.5); % gets the ending points of stimulations - -for i = 1:len - segment = filteredSignal(((i-1)*hop+1):((i-1)*hop+WSize)); - MAV_feature(i) = ; - VAR_feature(i) = ; - - % re-build the label vector to match it with the feature vector - featureLabels(i) = sum(arrayfun(@(t) ((i-1)*hop+1) >= Rise1(t) && ((i-1)*hop+WSize) <= Fall1(t), 1:length(Rise1))); -end - -%% Plotting the features -% Note: when plotting the features, scale the featureLabels to the max of -% the feature values for proper visualization +% Create Figure +figure('Name', 'Feature Extraction Sweep with MAV & VAR SNR', 'units', 'normalized', 'Position', [0, 0, 1, 1]); +plot_idx = 1; % Counter for subplot index +%% Loop through all combinations +for w = 1:length(WSize_values) + for o = 1:length(Olap_values) + + % Current parameters + current_WSize_s = WSize_values(w); + current_Olap_pct = Olap_values(o); + + % Windowing calculations + WSize_samp = floor(current_WSize_s * fs); % Window size in samples + nOlap = floor(current_Olap_pct * WSize_samp); % Overlap in samples + hop = WSize_samp - nOlap; % Hop size + nx = length(filteredSignal); + len = fix((nx - (WSize_samp - hop)) / hop); % Total number of frames + + % Preallocate + MAV_feature = zeros(1, len); + VAR_feature = zeros(1, len); + featureLabels = zeros(1, len); + + % Feature Extraction Loop + for i = 1:len + % Extract segment + idx_start = (i-1)*hop + 1; + idx_end = idx_start + WSize_samp - 1; + + % Check bounds + if idx_end > nx + break; + end + + segment = filteredSignal(idx_start:idx_end); + + % Calculate Features + MAV_feature(i) = mean(abs(segment)); + VAR_feature(i) = var(segment); + + % Re-build label vector + % Strict: Window must be fully inside stimulation to count as '1' + is_stim = any(idx_start >= Rise1 & idx_end <= Fall1); + featureLabels(i) = double(is_stim); + end + + %% Calculate SNR + % Separate Stimulus and Rest values using logical indexing + stim_indices = (featureLabels == 1); + rest_indices = (featureLabels == 0); + + % --- MAV SNR --- + mean_mav_stim = mean(MAV_feature(stim_indices)); + mean_mav_rest = mean(MAV_feature(rest_indices)); + if mean_mav_rest > 0 + mav_snr = 20 * log10(mean_mav_stim / mean_mav_rest); + else + mav_snr = 0; + end + + % --- VAR SNR --- + mean_var_stim = mean(VAR_feature(stim_indices)); + mean_var_rest = mean(VAR_feature(rest_indices)); + if mean_var_rest > 0 + var_snr = 20 * log10(mean_var_stim / mean_var_rest); + else + var_snr = 0; + end + + %% Plotting MAV (Left Column) + subplot(9, 2, plot_idx); + plot(MAV_feature, 'b', 'LineWidth', 0.5); hold on; + plot(featureLabels * max(MAV_feature), 'r', 'LineWidth', 1); + + % Formatting Title with SNR + title_str = sprintf('MAV (W=%.2g, Olap=%.2g) | SNR=%.2f dB', ... + current_WSize_s, current_Olap_pct, mav_snr); + grid on; + title(title_str, 'FontSize', 8); + + if mod(plot_idx, 2) == 1 + ylabel('MAV'); + end + xlim([1, len]); + set(gca, 'XTickLabel', []); + + plot_idx = plot_idx + 1; + + %% Plotting VAR (Right Column) + subplot(9, 2, plot_idx); + plot(VAR_feature, 'k', 'LineWidth', 0.5); hold on; + plot(featureLabels * max(VAR_feature), 'r', 'LineWidth', 1); + + % Formatting VAR with SNR + title_str_var = sprintf('VAR (W=%.2g, Olap=%.2g) | SNR=%.2f dB', ... + current_WSize_s, current_Olap_pct, var_snr); + grid on; + title(title_str_var, 'FontSize', 8); + + if mod(plot_idx, 2) == 0 + ylabel('VAR'); + end + xlim([1, len]); + set(gca, 'XTickLabel', []); + + plot_idx = plot_idx + 1; + end +end \ No newline at end of file diff --git a/HW_1_PNS_EE379K-385V_Neural_Eng/c3_classification.m b/HW_1_PNS_EE379K-385V_Neural_Eng/c3_classification.m index e12fb04..b96b80e 100644 --- a/HW_1_PNS_EE379K-385V_Neural_Eng/c3_classification.m +++ b/HW_1_PNS_EE379K-385V_Neural_Eng/c3_classification.m @@ -1,107 +1,178 @@ -close all;clc; -%% -% Inputs: -% -------- -% MAVClass1: the features of the VF case (stimulus and rest features) -% MAVClass2: the features of the Pinch case (stimulus and rest features) -% TriggerClass1: labels for VF features (stimulus or rest label) -% TriggerClass2: labels for Pinch features (stimulus or rest label) +%% c3_classification_complete.m +% This script performs 10-fold cross-validation for Part 2.3 of the assignment. +% It answers: +% 1. Classifiability of Stimuli vs Rest +% 2. Classifiability of Stimuli vs Stimuli +% 3. Comparison of MAV vs VAR features +% 4. Evaluation of Confusion Matrices -% Build the datasets -MAV_class1 = MAVClass1(find(TriggerClass1==1)); -MAV_rest1 = MAVClass1(find(TriggerClass1==0)); +clearvars -except filteredFlex filteredPinch filteredVF flexLabels pinchLabels VFLabels fs; +clc; -VAR_class1 = VARClass1(find(TriggerClass1==1)); -VAR_rest1 = VARClass1(find(TriggerClass1==0)); +% Check if data is loaded +if ~exist('filteredFlex', 'var') + error('Error: Filtered signals not found. Please run c1_dataVis.m first.'); +end -MAV_class2 = MAVClass2(find(TriggerClass2==1)); -MAV_rest2 = MAVClass2(find(TriggerClass2==0)); +%% 1. Feature Extraction (WSize=100ms, Olap=0) +fprintf('1. Extracting Features (100ms Window, 0%% Overlap)...\n'); -VAR_class2 = VARClass2(find(TriggerClass2==1)); -VAR_rest2 = VARClass2(find(TriggerClass2==0)); +WSize_sec = 0.1; +Olap_pct = 0; +WSize = floor(WSize_sec * fs); +nOlap = floor(Olap_pct * WSize); +hop = WSize - nOlap; -% Concantenate the rest classes -MAV_rest = [MAV_rest1 MAV_rest2]; -VAR_rest = [VAR_rest1 VAR_rest2]; +% Organize data for looping +% Index 1=VF, 2=Flex, 3=Pinch +sigs = {filteredVF, filteredFlex, filteredPinch}; +lbls = {VFLabels, flexLabels, pinchLabels}; +names = {'VF', 'Flex', 'Pinch'}; +feats = struct(); % Structure to hold features -%% -% Class1 vs Rest dataset -MAV_Data_Class1vsRest = [MAV_class1 MAV_rest]; -MAV_Labels_Class1vsRest = [ones(1,length(MAV_class1)) 2*ones(1,length(MAV_rest))]; +for k = 1:3 + sig = sigs{k}; + lab = lbls{k}; + + nx = length(sig); + len = fix((nx - (WSize - hop)) / hop); + + MAV_vec = zeros(1, len); + VAR_vec = zeros(1, len); + LBL_vec = zeros(1, len); + + Rise = gettrigger(lab, 0.5); + Fall = gettrigger(-lab, -0.5); + + for i = 1:len + idx_start = (i-1)*hop + 1; + idx_end = idx_start + WSize - 1; + segment = sig(idx_start:idx_end); + + MAV_vec(i) = mean(abs(segment)); + VAR_vec(i) = var(segment); + + % Label: 1 if window is strictly inside stimulation + is_stim = any(idx_start >= Rise & idx_end <= Fall); + LBL_vec(i) = double(is_stim); + end + + feats(k).MAV = MAV_vec; + feats(k).VAR = VAR_vec; + feats(k).LBL = LBL_vec; + feats(k).Name = names{k}; +end -VAR_Data_Class1vsRest = [VAR_class1 VAR_rest]; -VAR_Labels_Class1vsRest = MAV_Labels_Class1vsRest; +%% 2. Define Comparisons +% We need to run classification for these specific pairs: +comparisons = { + 'VF vs Rest', 1, 0; % 0 denotes "Rest" class + 'Flex vs Rest', 2, 0; + 'Pinch vs Rest', 3, 0; + 'Flex vs Pinch', 2, 3; + 'Flex vs VF', 2, 1; + 'Pinch vs VF', 3, 1; +}; -% Class2 vs Rest dataset -MAV_Data_Class2vsRest = [MAV_class2 MAV_rest]; -MAV_Labels_Class2vsRest = [ones(1,length(MAV_class2)) 2*ones(1,length(MAV_rest))]; +%% 3. Classification Loop (10-Fold CV) +fprintf('\n2. Running 10-Fold Cross Validation...\n'); +fprintf('----------------------------------------------------------------\n'); +fprintf('%-20s | %-12s | %-12s | %-15s\n', 'Comparison', 'Acc (MAV)', 'Acc (VAR)', 'Best Feature'); +fprintf('----------------------------------------------------------------\n'); -VAR_Data_Class2vsRest = [VAR_class2 VAR_rest]; -VAR_Labels_Class2vsRest = MAV_Labels_Class2vsRest; +for c = 1:size(comparisons, 1) + comp_name = comparisons{c, 1}; + idx1 = comparisons{c, 2}; + idx2 = comparisons{c, 3}; + + % --- Prepare Data for Class 1 --- + % Get Stimulus features (Label == 1) + f1_MAV = feats(idx1).MAV(feats(idx1).LBL == 1); + f1_VAR = feats(idx1).VAR(feats(idx1).LBL == 1); + + % --- Prepare Data for Class 2 (or Rest) --- + if idx2 == 0 + % If comparing vs Rest, get Rest features (Label == 0) from the SAME signal + f2_MAV = feats(idx1).MAV(feats(idx1).LBL == 0); + f2_VAR = feats(idx1).VAR(feats(idx1).LBL == 0); + label_names = {feats(idx1).Name, 'Rest'}; + else + % If comparing vs another Stimulus, get Stimulus features (Label == 1) + f2_MAV = feats(idx2).MAV(feats(idx2).LBL == 1); + f2_VAR = feats(idx2).VAR(feats(idx2).LBL == 1); + label_names = {feats(idx1).Name, feats(idx2).Name}; + end + + % Combine Data + X_MAV = [f1_MAV, f2_MAV]'; % Transpose to column vector + X_VAR = [f1_VAR, f2_VAR]'; + + % Create Labels (1 for Class 1, 2 for Class 2) + Y = [ones(length(f1_MAV), 1); 2 * ones(length(f2_MAV), 1)]; + + % --- 10-Fold Cross Validation --- + k = 10; + cv = cvpartition(Y, 'KFold', k); % Random split (answering Q4!) + + acc_mav = 0; + acc_var = 0; + conf_mav = zeros(2,2); % Accumulate confusion matrix + + for i = 1:k + train_idx = cv.training(i); + test_idx = cv.test(i); + + % MAV Classification + pred_mav = classify(X_MAV(test_idx), X_MAV(train_idx), Y(train_idx)); + acc_mav = acc_mav + sum(pred_mav == Y(test_idx)) / length(pred_mav); + + % Build Confusion Matrix for MAV (just one example needed for assignment) + % Rows = True Class, Cols = Predicted Class + current_conf = confusionmat(Y(test_idx), pred_mav); + % Handle edge case if a fold misses a class + if size(current_conf,1) == 2 + conf_mav = conf_mav + current_conf; + end -% Class1 vs Class2 dataset -MAV_Data_Class1vsClass2 = [MAV_class1 MAV_class2]; -MAV_Labels_Class1vsClass2 = [ones(1,length(MAV_class1)) 2*ones(1,length(MAV_class2))]; + % VAR Classification + pred_var = classify(X_VAR(test_idx), X_VAR(train_idx), Y(train_idx)); + acc_var = acc_var + sum(pred_var == Y(test_idx)) / length(pred_var); + end + + % Average Accuracy + mean_acc_mav = (acc_mav / k) * 100; + mean_acc_var = (acc_var / k) * 100; + + % Determine Winner + if mean_acc_mav > mean_acc_var + winner = 'MAV'; + elseif mean_acc_var > mean_acc_mav + winner = 'VAR'; + else + winner = 'Tie'; + end + + fprintf('%-20s | %-11.1f%% | %-11.1f%% | %-15s\n', comp_name, mean_acc_mav, mean_acc_var, winner); + + % --- Display Confusion Matrix for MAV --- + % Only printing logic to keep output clean, answering "Observe confusion matrices" + fprintf(' Confusion Matrix (MAV) for %s:\n', comp_name); + fprintf(' True %-6s: [ %4d %4d ] (Predicted %s / %s)\n', label_names{1}, conf_mav(1,1), conf_mav(1,2), label_names{1}, label_names{2}); + fprintf(' True %-6s: [ %4d %4d ]\n\n', label_names{2}, conf_mav(2,1), conf_mav(2,2)); +end +fprintf('----------------------------------------------------------------\n'); -VAR_Data_Class1vsClass2 = [VAR_class1 VAR_class2]; -VAR_Labels_Class1vsClass2 = MAV_Labels_Class1vsClass2; - -%% -% Both feature datasets -MAVVAR_Data_Class1vsRest = [MAV_Data_Class1vsRest; VAR_Data_Class1vsRest]; -MAVVAR_Labels_Class1vsRest = MAV_Labels_Class1vsRest; - -MAVVAR_Data_Class2vsRest = [MAV_Data_Class2vsRest; VAR_Data_Class2vsRest]; -MAVVAR_Labels_Class2vsRest = MAV_Labels_Class2vsRest; - -MAVVAR_Data_Class1vsClass2 = [MAV_Data_Class1vsClass2; VAR_Data_Class1vsClass2]; -MAVVAR_Labels_Class1vsClass2 = MAV_Labels_Class1vsClass2; - -%% -% Classify all combinations (training set) -k = 10; % for k-fold cross validation -c1 = cvpartition(length(MAV_Labels_Class1vsRest),'KFold',k); -c2 = cvpartition(length(VAR_Labels_Class1vsRest),'KFold',k); -c3 = cvpartition(length(MAVVAR_Labels_Class1vsRest),'KFold',k); -c4 = cvpartition(length(MAV_Labels_Class2vsRest),'KFold',k); -c5 = cvpartition(length(VAR_Labels_Class2vsRest),'KFold',k); -c6 = cvpartition(length(MAVVAR_Labels_Class2vsRest),'KFold',k); -c7 = cvpartition(length(MAV_Labels_Class1vsClass2),'KFold',k); -c8 = cvpartition(length(VAR_Labels_Class1vsClass2),'KFold',k); -c9 = cvpartition(length(MAVVAR_Labels_Class1vsClass2),'KFold',k); - -% Repeat the following for i=1:k, and average performance metrics across all iterations -i=1; -% loop over all k-folds and avergae the performance -% for i=1:k - [TstMAVFC1Rest TstMAVErrC1Rest] = classify(MAV_Data_Class1vsRest(c1.test(i))',MAV_Data_Class1vsRest(c1.training(i))',MAV_Labels_Class1vsRest(c1.training(i))); - [TstCM_MAV_C1rest dum1 TstAcc_MAV_C1rest dum2] = confusion(MAV_Labels_Class1vsRest(c1.test(i)), TstMAVFC1Rest); - - [TstVARFC1Rest TstVARErrC1Rest] = classify(VAR_Data_Class1vsRest(c2.test(i))',VAR_Data_Class1vsRest(c2.training(i))',VAR_Labels_Class1vsRest(c2.training(i))); - [TstCM_VAR_C1rest dum1 TstAcc_VAR_C1rest dum2] = confusion(VAR_Labels_Class1vsRest(c2.test(i)), TstVARFC1Rest); - - [TstMAVVARFC1Rest TstMAVVARErrC1Rest] = classify(MAVVAR_Data_Class1vsRest(:,c3.test(i))',MAVVAR_Data_Class1vsRest(:,c3.training(i))',MAVVAR_Labels_Class1vsRest(c3.training(i))); - [TstCM_MAVVAR_C1rest dum1 TstAcc_MAVVAR_C1rest dum2] = confusion(MAVVAR_Labels_Class1vsRest(c3.test(i)), TstMAVVARFC1Rest); - - % Class2 vs Rest - [TstMAVFC2Rest TstMAVErrC2Rest] = classify(MAV_Data_Class2vsRest(c4.test(i))',MAV_Data_Class2vsRest(c4.training(i))',MAV_Labels_Class2vsRest(c4.training(i))); - [TstCM_MAV_C2rest dum1 TstAcc_MAV_C2rest dum2] = confusion(MAV_Labels_Class2vsRest(c4.test(i)), TstMAVFC2Rest); - - [TstVARFC2Rest TstVARErrC2Rest] = classify(VAR_Data_Class2vsRest(c5.test(i))',VAR_Data_Class2vsRest(c5.training(i))',VAR_Labels_Class2vsRest(c5.training(i))); - [TstCM_VAR_C2rest dum1 TstAcc_VAR_C2rest dum2] = confusion(VAR_Labels_Class2vsRest(c5.test(i)), TstVARFC2Rest); - - [TstMAVVARFC2Rest TstMAVVARErrC2Rest] = classify(MAVVAR_Data_Class2vsRest(:,c6.test(i))',MAVVAR_Data_Class2vsRest(:,c6.training(i))',MAVVAR_Labels_Class2vsRest(c6.training(i))); - [TstCM_MAVVAR_C2rest dum1 TstAcc_MAVVAR_C2rest dum2] = confusion(MAVVAR_Labels_Class2vsRest(c6.test(i)), TstMAVVARFC2Rest); - - % Class1 vs Class2 - [TstMAVFC1C2 TstMAVErrC1C2] = classify(MAV_Data_Class1vsClass2(c7.test(i))',MAV_Data_Class1vsClass2(c7.training(i))',MAV_Labels_Class1vsClass2(c7.training(i))); - [TstCM_MAV_C1C2 dum1 TstAcc_MAV_C1C2 dum2] = confusion(MAV_Labels_Class1vsClass2(c7.test(i)), TstMAVFC1C2); - - [TstVARFC1C2 TstVARErrC1C2] = classify(VAR_Data_Class1vsClass2(c8.test(i))',VAR_Data_Class1vsClass2(c8.training(i))',VAR_Labels_Class1vsClass2(c8.training(i))); - [TstCM_VAR_C1C2 dum1 TstAcc_VAR_C1C2 dum2] = confusion(VAR_Labels_Class1vsClass2(c8.test(i)), TstVARFC1C2); - - [TstMAVVARFC1C2 TstMAVVARErrC1C2] = classify(MAVVAR_Data_Class1vsClass2(:,c9.test(i))',MAVVAR_Data_Class1vsClass2(:,c9.training(i))',MAVVAR_Labels_Class1vsClass2(c9.training(i))); - [TstCM_MAVVAR_C1C2 dum1 TstAcc_MAVVAR_C1C2 dum2] = confusion(MAVVAR_Labels_Class1vsClass2(c9.test(i)), TstMAVVARFC1C2); -% end -%% +%% 4. Answer Prompts +fprintf('\n=== Automated Analysis for Part 2.3 ===\n'); +fprintf('1. Check "Pinch vs Rest" accuracy above. Is it low? (Likely yes, due to low SNR).\n'); +fprintf('2. Check "Flex vs Pinch". Can they be distinguished?\n'); +fprintf('3. Observe the Confusion Matrices: Are they balanced? \n'); +fprintf(' - If one class is predicted much more often, the classifier is biased.\n'); +fprintf('4. Feature Performance: Look at the "Best Feature" column.\n'); +fprintf(' - MAV is typically more robust for these signals.\n'); +fprintf('5. Validation Fairness (Assignment Q4):\n'); +fprintf(' - This script uses "cvpartition", which splits data RANDOMLY.\n'); +fprintf(' - Since EMG/ENG signals are time-series, random splitting causes "data leakage"\n'); +fprintf(' (training on samples immediately adjacent to test samples).\n'); +fprintf(' - Therefore, this is likely NOT a fair assessment of generalization.\n'); \ No newline at end of file diff --git a/HW_2_sp23_EE374N-385J_Neural_Eng_EMGanalysis/HW_2_sp22_EE379K-385V_Neural_Eng_EMGanalysis.pdf b/HW_2_sp23_EE374N-385J_Neural_Eng_EMGanalysis/HW_2_sp22_EE379K-385V_Neural_Eng_EMGanalysis.pdf new file mode 100644 index 0000000..0cfb07a Binary files /dev/null and b/HW_2_sp23_EE374N-385J_Neural_Eng_EMGanalysis/HW_2_sp22_EE379K-385V_Neural_Eng_EMGanalysis.pdf differ diff --git a/HW_2_sp23_EE374N-385J_Neural_Eng_EMGanalysis/ppt_HW_2_EMG_EE374N-385J_Neural_Eng.pdf b/HW_2_sp23_EE374N-385J_Neural_Eng_EMGanalysis/ppt_HW_2_EMG_EE374N-385J_Neural_Eng.pdf new file mode 100644 index 0000000..7de6109 Binary files /dev/null and b/HW_2_sp23_EE374N-385J_Neural_Eng_EMGanalysis/ppt_HW_2_EMG_EE374N-385J_Neural_Eng.pdf differ diff --git a/HW_2_sp23_EE374N-385J_Neural_Eng_EMGanalysis/subject1.mat b/HW_2_sp23_EE374N-385J_Neural_Eng_EMGanalysis/subject1.mat new file mode 100644 index 0000000..e0595ea Binary files /dev/null and b/HW_2_sp23_EE374N-385J_Neural_Eng_EMGanalysis/subject1.mat differ diff --git a/HW_2_sp23_EE374N-385J_Neural_Eng_EMGanalysis/subject2.mat b/HW_2_sp23_EE374N-385J_Neural_Eng_EMGanalysis/subject2.mat new file mode 100644 index 0000000..44436cb Binary files /dev/null and b/HW_2_sp23_EE374N-385J_Neural_Eng_EMGanalysis/subject2.mat differ