fin

.gitignore

@@ -205,3 +205,5 @@ cython_debug/
 marimo/_static/
 marimo/_lsp/
 __marimo__/
+
+.DS_Store

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;
+flexSignal=Flex.signal;
-labels=Flex.trigger;
+pinchSignal = Pinch.signal;
-TRIG = gettrigger(labels,0.5);
+VFSignal = VF.signal;
-TRIGend = gettrigger(-labels,-0.5);
+flexLabels=Flex.trigger;
+pinchLabels=Pinch.trigger;
+VFLabels = VF.trigger;
 
-figure('units','normalized','Position',[0.1,0.1,0.7,0.4])
+%% Plot the raw signals
-plot((1:length(signal))./fs,zscore(signal));
+% Grouping variables into cell arrays to loop through them cleanly
-hold on;
+signals = {flexSignal, pinchSignal, VFSignal};
-plot((1:length(signal))./fs,zscore(labels),'y');
+labels_all = {flexLabels, pinchLabels, VFLabels};
-stem(TRIG./fs,ones(length(TRIG),1)*max(zscore(labels)),'Color','g');
+titles = {'Raw Flex Signal with Stimulation Pattern (Yellow)', ...
-stem(TRIGend./fs,ones(length(TRIG),1)*max(zscore(labels)),'Color','r');
+          'Raw Pinch Signal with Stimulation Pattern (Yellow)', ...
-grid on; grid minor;
+          'Raw VF Signal with Stimulation Pattern (Yellow)'};
-xlim([0,length(signal)./fs])
+
-xlabel('Time (s)')
+% Create one large figure for all three subplots
-ylabel('Amplitude (uV)')
+figure('units','normalized','Position',[0.1, 0.1, 0.7, 0.8])
-title('Raw VF signal with labels for stimulation periods')
+
+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])
+figure('Name', 'Raw PSD Estimates', 'units','normalized','Position',[0.1,0.1,0.5,0.5])
-[rows_act,cols_act,values_act] = find(labels>0);
+h = spectrum.welch; 
-[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';
 
-% %% 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 
+fc1 = 800; % first cutoff frequency in Hz 
-% fc2 = 100; % second 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

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');

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
+
+

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 
+% Sweep parameters
-label = ; % labels of stimulus locations
+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
+% Get trigger points once (indices in the raw signal)
-Olap = ; % overlap percentage
+Rise1 = gettrigger(label, 0.5);   % Start of stimulation
+Fall1 = gettrigger(-label, -0.5); % End of stimulation
 
-%% Extracting Features over overlapping windows
+% Create Figure
-
+figure('Name', 'Feature Extraction Sweep with MAV & VAR SNR', 'units', 'normalized', 'Position', [0, 0, 1, 1]);
-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
 
+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

HW_1_PNS_EE379K-385V_Neural_Eng/c3_classification.m

@@ -1,107 +1,178 @@
-close all;clc;
+%% c3_classification_complete.m
-%%
+% This script performs 10-fold cross-validation for Part 2.3 of the assignment.
-% Inputs: 
+% It answers:
-% --------
+% 1. Classifiability of Stimuli vs Rest
-% MAVClass1: the features of the VF case (stimulus and rest features)
+% 2. Classifiability of Stimuli vs Stimuli
-% MAVClass2: the features of the Pinch case (stimulus and rest features)
+% 3. Comparison of MAV vs VAR features
-% TriggerClass1: labels for VF features (stimulus or rest label)
+% 4. Evaluation of Confusion Matrices
-% TriggerClass2: labels for Pinch features (stimulus or rest label)
 
-% Build the datasets
+clearvars -except filteredFlex filteredPinch filteredVF flexLabels pinchLabels VFLabels fs;
-MAV_class1 = MAVClass1(find(TriggerClass1==1));
+clc;
-MAV_rest1 = MAVClass1(find(TriggerClass1==0));
 
-VAR_class1 = VARClass1(find(TriggerClass1==1));
+% Check if data is loaded
-VAR_rest1 = VARClass1(find(TriggerClass1==0));
+if ~exist('filteredFlex', 'var')
+    error('Error: Filtered signals not found. Please run c1_dataVis.m first.');
+end
 
-MAV_class2 = MAVClass2(find(TriggerClass2==1));
+%% 1. Feature Extraction (WSize=100ms, Olap=0)
-MAV_rest2 = MAVClass2(find(TriggerClass2==0));
+fprintf('1. Extracting Features (100ms Window, 0%% Overlap)...\n');
 
-VAR_class2 = VARClass2(find(TriggerClass2==1));
+WSize_sec = 0.1; 
-VAR_rest2 = VARClass2(find(TriggerClass2==0));
+Olap_pct = 0; 
+WSize = floor(WSize_sec * fs);
+nOlap = floor(Olap_pct * WSize);
+hop = WSize - nOlap;
 
-% Concantenate the rest classes
+% Organize data for looping
-MAV_rest = [MAV_rest1 MAV_rest2];
+% Index 1=VF, 2=Flex, 3=Pinch
-VAR_rest = [VAR_rest1 VAR_rest2];
+sigs = {filteredVF, filteredFlex, filteredPinch};
+lbls = {VFLabels, flexLabels, pinchLabels};
+names = {'VF', 'Flex', 'Pinch'};
 
+feats = struct(); % Structure to hold features
 
-%%
+for k = 1:3
-% Class1 vs Rest dataset
+    sig = sigs{k};
-MAV_Data_Class1vsRest = [MAV_class1 MAV_rest];
+    lab = lbls{k};
-MAV_Labels_Class1vsRest = [ones(1,length(MAV_class1)) 2*ones(1,length(MAV_rest))];
+    
+    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];
+%% 2. Define Comparisons
-VAR_Labels_Class1vsRest = MAV_Labels_Class1vsRest;
+% 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
+%% 3. Classification Loop (10-Fold CV)
-MAV_Data_Class2vsRest = [MAV_class2 MAV_rest];
+fprintf('\n2. Running 10-Fold Cross Validation...\n');
-MAV_Labels_Class2vsRest = [ones(1,length(MAV_class2)) 2*ones(1,length(MAV_rest))];
+fprintf('----------------------------------------------------------------\n');
+fprintf('%-20s | %-12s | %-12s | %-15s\n', 'Comparison', 'Acc (MAV)', 'Acc (VAR)', 'Best Feature');
+fprintf('----------------------------------------------------------------\n');
 
-VAR_Data_Class2vsRest = [VAR_class2 VAR_rest];
+for c = 1:size(comparisons, 1)
-VAR_Labels_Class2vsRest = MAV_Labels_Class2vsRest;
+    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
+        % VAR Classification
-MAV_Data_Class1vsClass2 = [MAV_class1 MAV_class2];
+        pred_var = classify(X_VAR(test_idx), X_VAR(train_idx), Y(train_idx));
-MAV_Labels_Class1vsClass2 = [ones(1,length(MAV_class1)) 2*ones(1,length(MAV_class2))];
+        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];
+%% 4. Answer Prompts
-VAR_Labels_Class1vsClass2 = MAV_Labels_Class1vsClass2;
+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');
-% Both feature datasets
+fprintf('3. Observe the Confusion Matrices: Are they balanced? \n');
-MAVVAR_Data_Class1vsRest = [MAV_Data_Class1vsRest; VAR_Data_Class1vsRest];
+fprintf('   - If one class is predicted much more often, the classifier is biased.\n');
-MAVVAR_Labels_Class1vsRest = MAV_Labels_Class1vsRest;
+fprintf('4. Feature Performance: Look at the "Best Feature" column.\n');
-
+fprintf('   - MAV is typically more robust for these signals.\n');
-MAVVAR_Data_Class2vsRest = [MAV_Data_Class2vsRest; VAR_Data_Class2vsRest];
+fprintf('5. Validation Fairness (Assignment Q4):\n');
-MAVVAR_Labels_Class2vsRest = MAV_Labels_Class2vsRest;
+fprintf('   - This script uses "cvpartition", which splits data RANDOMLY.\n');
-
+fprintf('   - Since EMG/ENG signals are time-series, random splitting causes "data leakage"\n');
-MAVVAR_Data_Class1vsClass2 = [MAV_Data_Class1vsClass2; VAR_Data_Class1vsClass2];
+fprintf('     (training on samples immediately adjacent to test samples).\n');
-MAVVAR_Labels_Class1vsClass2 = MAV_Labels_Class1vsClass2;
+fprintf('   - Therefore, this is likely NOT a fair assessment of generalization.\n');
-
-%%
-% 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
-%%