Ulty Multy - EEG Data Analysis and Visualization
================================================

Welcome to the tutorial of Ulty Multy!

This project utilizes the open data from OpenNeuro to analyze and visualize EEG data.

For detailed information, please visit the [Github Repo](https://github.com/mixuanpan/Ulty_Multy).

Foreword
--------

Electroencephalogram (EEG) data analysis plays a crucial role in understanding brain activity, enabling advancements in neuroscience, biomedical engineering, and brain-computer interfaces. This user manual serves as a comprehensive guide to Ulty Multy, a MATLAB and EEGLAB-based system designed to streamline EEG data processing, feature extraction, and visualization.

Developed to improve the efficiency and accuracy of EEG analysis, this tool integrates automated preprocessing workflows, advanced signal processing techniques, and intuitive graphical displays. Whether you are a researcher, student, or practitioner working with EEG data, this manual provides detailed instructions on system functionalities, from data loading to final interpretation.

By following this guide, users will gain a clear understanding of how to leverage the tool for extracting meaningful insights from EEG recordings. We hope this manual enhances your experience and facilitates a deeper exploration of brainwave patterns and cognitive analysis.

Main Script
-----------

This is the file to run to start the program

% asks for user consent to perform data analysis 
consent = ip_welcome(); 

% When Consent is given 
if consent
    % load subject data
    par = readtable("data/participants.tsv", "FileType","text", "Delimiter", "\t"); 
    
    % display user info 
    subject_num = ip_subject_info(par); 
    
    % filter data 
    [data, low, upper, noise_num] = ip_data_filter(subject_num); 
    
    % perform time frequency analysis and visualize the user-selected 
    % channel(s)
    [power, freqs1] = ip_time_frequency_analysis(data); 
    
    % view the relative location of each channel in a 3D interactive graph
    ip_brain_plot(data, subject_num); 
    
    % compare the data of the target subject to the entire subject pool in 
    % each channel
    ip_final_display(subject_num, low, upper, noise_num, freqs1, power); 
    fprintf("\nAnalysis Successful.\n\n");

% When consent is not given 
else 
    fprintf("\nPermission Denied!\n");
    fprintf("\nAnalysis Unsuccessful.\n\n");
end

Welcome Function
----------------

This function begins with a welcome message explaining the data collection process, providing users with background information to help interpret the analysis results. The program then asks for consent to analyze the user’s local data. If consent is denied, the program ends without performing any analysis. While the open data used for this project is stored locally, requesting consent is essential when handling real user data in future applications.

function consent = ip_welcome()
    % INITIALIZATION   
    fileName = "README"; % the name of the file to be displayed 
    dir = "data"; % file directory 
    filePath = fullfile(pwd, dir, fileName); % path to the file 
    
    % check if the file exists 
    if isfile(filePath) 
        content = fileread(filePath); % read the file as a char array 
    else
        content = '';
        fprintf("\nFile not found!\n");
    end
    
    % Printed Output 
    fprintf("Welcome to the EEG Data Analysis Tool.\n")
    fprintf("This tool is designed to perform different analysis on the given data as follows: \n\n")
    disp(content); % Data Collection Process 

    % Ask for user's consent 
    consent = input("Do you allow this tool to perform data analysis (y/n)? ", "s");
    if strcmp(lower(consent), 'y')
        fprintf("\nThank you for your permission!\n")
        fprintf("A series of analysis will be performed shortly.\n");
        consent = true;
    else 
        consent = false;
    end
end

### Welcome Message Displayed:
Welcome to the EEG Data Analysis Tool.
This tool is designed to perform different analysis on the given data as follows:

Data collection took place at the NeuroCognition Laboratory (NCL) in San Diego, California under the supervision of Dr. Phillip Holcomb.
This project followed the San Diego State University’s IRB guidelines.

Participants sat in a comfortable chair in a darkened sound attenuated room throughout the experiment.
They were given a gamepad for button pressing.
They were instructed to watch the LCD video monitor that was at a viewing distance of 150cm.

Participants were presented with 300 prime-target pairs.
All targets were four-letter English words.
Of the 300 critical trials, 100 had English word primes, 100 had ASL sign primes, and 100 had fingerspelled word primes.
Half of the primes in each condition were related to the targets.
Related English word primes were identity primes to the English word, related fingerspelled word primes were also identity primes, and related ASL primes were ASL translations of the English word targets.
The other half of the primes were unrelated to the targets.

Participants were instructed to focus on the purple fixation cross that appeared on the screen for 800ms.
This fixation cross then turned white for 500ms.
Then, one of three prime conditions was presented: an English word, an ASL sign, or a fingerspelled word.
English prime words were presented for 300ms.
Signed (M = 565ms) and fingerspelled (M = 1173ms) video primes had variable durations.
All target stimuli were 4-letter English words presented for 500ms.
Related primes were either identity or translations.

Press any of the 4 buttons on the right of the gamepad whenever you see an animal.
It doesn’t matter if the animal is presented as a sign, a word, or fingerspelled.
Press for ANY animal. You can blink whenever you see purple. A purple + means you have time for a quick blink.
A purple (--) means you can blink as much as you want.

Subject Number & Info Function
------------------------------

This function prompts the user to select a target subject (between 1 and 27) for analysis. If the input is invalid or out of range, the program will continue asking until a valid input is provided. The function then reads a .tsv file to display the selected subject’s demographics and data collection settings. This helps the user identify potential influencing factors when comparing data across subjects.

function subject_num = ip_subject_info(par)
    
    % INITIALIZATION
    invalid_num = [14 15 25]; % no subjects found 
    len = length(invalid_num); % number of invalid numbers 
    checker = true; % while loop control variable 
    bool = true; % loop control variable 
    
    % subject_num error check 
    while checker 
        % ask for subject number from user input
        subject_num = input("Which subject are you looking for? Enter a number between 1 and 27 (inclusive) --> "); 
        
        % When the subject number is in the given range 
        if (subject_num >= 1 && subject_num <= 27) 
            for i = 1:len % length of invalid_num 
                if subject_num == invalid_num(i)
                    fprintf("\nSubject not found!\n\n");
                    bool = false; 
                end
            end
            if bool
                checker = false; % stop the while loop 
            end
            bool = true;
        
        % When the subject number is out of range 
        else 
            fprintf("\nError! Enter a number between 1 to 27 (inclusive)!\n\n");
        end 
    end 
    
    % target subject ID 
    id = "sub-" + int2str(subject_num); 
    
    % print out subject info 
    fprintf("\n<strong>Subject Info:</strong>\n");

    % match: similar to the "find" function but comparing in tables 
    disp(par(matches(par.participant_id, id), :)); 

    % print out subject settings 
    fprintf("\n<strong>Subject Settings:</strong>\n");
    fileName = fullfile("data", sprintf("sub-%d", subject_num), "eeg", sprintf("sub-%d_task-SemanticCategorization_eeg.json", subject_num));
    jsonData = jsondecode(fileread(fileName));
    disp(jsonData); 
end 

### Subject Info Displayed

Data Filtering Function
-----------------------

This function includes bandpass filtering and notch filtering, both of which help reduce noise and improve signal clarity in the EEG data.

Bandpass filtering isolates signals within a specific frequency range, removing frequencies outside of this range to focus on meaningful brain activity. The user is prompted to input lower and upper bounds between 0.1 Hz and 40 Hz — the typical range for EEG analysis, as it captures key brainwave patterns such as delta, theta, alpha, beta, and gamma waves. If the user inputs values outside this range or sets the upper bound lower than the lower bound, the program triggers an error and asks for valid inputs. Once the program receives valid values, it uses the pop\_eegfiltnew() function to apply the bandpass filter, preserving only the signals within the chosen frequency band. For instance, setting the range to 1–30 Hz would eliminate low-frequency drifts and high-frequency noise, making brainwave patterns more discernible.

Notch filtering, on the other hand, targets narrow-band interference, such as powerline noise, which typically occurs at 50 Hz (in most countries) or 60 Hz (in North America). The user can select either 50 Hz or 60 Hz for the notch filter, and the program will apply it to remove the constant-frequency interference that often appears as spikes in EEG data. Removing powerline noise helps uncover subtle patterns in brain activity that might otherwise be masked by electrical interference. By applying both bandpass and notch filtering, the program ensures that the processed EEG data is cleaner and more suitable for accurate time- frequency analysis.

function [data, low, upper, noise_num] = ip_data_filter(subject_num)

    fprintf("\n<strong>EEG Data Filtering Begins</strong>\n"); 
    currentPath = "/Users/mixuan/Desktop/Ulty_Multy"; % current directory 
    
    % the file for reading EEG data
    fileName = "sub-" + int2str(subject_num) + "_task-SemanticCategorization_eeg.set"; 

    % the path for the EEG data file 
    filePath = fullfile(currentPath, "data", "sub-" + int2str(subject_num), "eeg"); 
    
    cd(filePath); % change to file directory 

    % load eeg data list from the data of user-chosen subject
    % use char() to avoid '.set' being truncated 
    EEG = pop_loadset(char(fileName)); 
        
    cd(currentPath); % change to working directory 
    checker = true; % while loop control variable 
    checker2 = true; % another loop control variable 
    low = 0; % initialize the lower bound
    upper = 0; % initialize the upper bound 
    noise_num = 0; % initialize the powerline to be removed 

    % bandpass filter bounds from user input 
    while checker 
        indicator = input("Do you want to bandpass filter (y/n)? ", 's');
        
        % When the user wants a bandpass filter 
        if strcmp(lower(indicator), "y")

            % ask for lower bound input 
            while checker2
                low = input("Enter the lower bound in Hz: ");
                if low < 0.1
                    fprintf("Entered lower bound is too small!\n")
                    fprintf("Please enter a positive number greater or equal to 0.1 Hz!\n")
                else
                    checker2 = false;
                end
            end

            % ask for upper bound input 
            checker2 = true;
            while checker2
                upper = input("Enter the upper bound in Hz: ");
                if upper > 40
                    fprintf("Entered upper bound is too high!\n")
                    fprintf("Please enter a positive number smaller than 40 Hz!\n")
                elseif upper < low
                    fprintf("Entered upper bound is too low!\n")
                    fprintf("Upper bound must be higher than lower bound!\n")
                else
                    checker2 = false;
                end
            end

            % Bandpass filter in the desired Hz range 
            EEG = pop_eegfiltnew(EEG, low, upper); 

            checker = false;
        % When the user doesn't want a bandpass filter 
        elseif strcmp(lower(indicator), "n")
            checker = false;
        end
    end
    
    % noise removal 
    checker = true;
    while checker 
        indicator = input("Do you want to remove a specific powerline (y/n)? ", 's');

        % When the user wants a noise removal 
        if strcmp(lower(indicator), 'y')
            noise_num = menu("Enter the notch filter to be filtered out", "50 Hz (common in Europe and most countries)", "60 Hz (common in North America)"); 
            
            % convert 1, 2 (index) to 50, 60 (notch filter)
            noise_num = 40 + 10 * noise_num; 

            % notch filtering 
            EEG = pop_eegfiltnew(EEG, noise_num - 1, noise_num + 1, [], 1); 
            checker = false;
        
        % When the user doesn't want a noise removal     
        elseif strcmp(lower(indicator), 'n')
            checker = false;
        
        % When the user falsely indicates if they want a noise removal or
        % not, the program skips the noise removal processing 
        else
            fprintf("Please enter either 'y' or 'n'\n"); 
        end
    end
    
    % function output 
    data = EEG; 
    fprintf("\nEEG Data Filtering Successful\n");
end 

### Notch Filtering Options
### Data Filtering Full Display

Time Frequency Analysis Function
--------------------------------

After filtering, the program automatically conducts a time-frequency analysis using FFT with Hanning window tapering. It adjusts the frequency range to match FFT output (1 Hz to 60 Hz) and generates 200 time points. This process is repeated 32 times — once for each channel — producing 30,934 estimated frequencies.

The function generates a heat map to visualize how the power of different frequency components changes over time for each EEG channel. The user can select any of the 32 channels for display, and the program will create a heat map using the imagesc() function. In the heat map, the x-axis represents time, the y-axis represents frequency, and the color intensity indicates the magnitude of power at each time-frequency point. Warmer colors, such as red and yellow, reflect higher power, while cooler colors, like blue and green, indicate lower power. In this example, the user chooses channel 21, 7, 1, and 32 to visualize.

function [power, freqs] = ip_time_frequency_analysis(input_data)

    % inform the user of the analysis type 
    fprintf("\n<strong>Time Frequency Analysis</strong>\n"); 

    EEG = input_data;
    channels_num = size(EEG.data, 1); % number of channels 
    
    % cell arrays contain elements with different types 
    power = cell(1, channels_num); % a cell array to store power 
        
    freqs = []; % frequency initialization 
    checker = true; % loop control variable 

% Wavelet Decomposition for each (typically 32) channel 
    for ch = 1:channels_num
        [power{ch}, freqs] = timefreq(EEG.data(ch, :), EEG.srate, 'freqs', [1, 60]);
    end

% Power Spectrum for individual channels 
    while checker 
        checker2 = true; % loop control variable 
        while checker2 
            channel = input("\n\nWhich channel would you like to visualize? ");
            if ~isnumeric(channel)
                fprintf("\nPlease enter a numeric value!\n");
            elseif isnumeric(channel) & (channel < 1 | channel > channels_num)
                fprintf("\nPlease enter a number between 1 and %i!\n", channels_num);
            else
                checker2 = false;
            end
        end

        % Plotting 
        figure;
        imagesc([], freqs, abs(power{channel})); % plot power magnitude 
        xlabel("Time (ms)");
        ylabel("Frequency (Hz)");
        title(["Time-Frequency Representation: Channel ", num2str(channel)]);
        colorbar;
        
        % wait until the figure is closed by the user
        waitfor(gcf); 

        % Plot another channel 
        checker2 = input("Do you want to visualize another channel (y/n)? ", "s");
        if strcmp(lower(checker2), "n")
            checker = false;
        elseif ~strcmp(lower(checker2), "y")
            fprintf("\nInput Error!!\n\n");
            checker = false;
        end
    end
    fprintf("\nTime Frequency Analysis and Visualization Successful.\n");
end 

### Analysis for Each of 32 Channels
### User's Choices of Channel Visualization
### Time Frequency Analysis Visualization Examples

3D Brain Plot Function
----------------------

This function creates a 3D interactive brain plot where the user can rotate, zoom, and explore the brain's surface. Highlighted channels are shown in red, making it easy to track their position. The interactive nature allows the user to inspect different angles and better understand the spatial relationship between channels. This feature enhances data interpretation by providing a clear, dynamic view of brain activity distribution.

function ip_brain_plot(data, subject_num)
    
    % INITIALIZATION
    EEG = data; % EEG data from user input 
    chanlocs = EEG.chanlocs; % extract channel locations 
    checker = true; % loop control variable 
    
    % extract the channels for the given subject number from the file
    channels = readtable(sprintf("data/sub-%d/eeg/sub-%d_task-SemanticCategorization_channels.tsv", subject_num, subject_num), "FileType", "text", ... 
        "Delimiter", "\t");
    
    % extract the name column of the channels to show in the final 3D graph
    labels = channels.name; 

    fprintf("\n<strong>Brain Plotting Begins</strong>\n");

    % ask the user if they want to highlight a specific channel 
    while checker 
        indicator = input("Do you want to highlight a specific channel (y/n)? ", "s");
        if strcmp(lower(indicator), "n")
            checker = false; 
        elseif strcmp(lower(indicator), "y")
            checker = false;
            highlight = menu("Choose the channel name that you want it to be highlighted.", labels);
        else
            fprintf("Please enter 'y' or 'n'!\n");
        end
    end

    figure;
    hold on;

    % iterate through each EEG channel as a point 
    for i = 1:length(chanlocs)
        if ~isempty(labels{i}) % Check if the label exists

            % locates each channel to the correct coordination
            plot3(chanlocs(i).X, chanlocs(i).Y, chanlocs(i).Z); 
            
            % determine highlight 
            if exist("highlight", "var") & i == highlight
                % variable "highlight" doesn't exist if the user didn't
                % specify it 
                color = "red";
            else 
                % the color of the highlighted channel on the 3D brain map, 
                % default black 
                color = "black"; 
            end 

            % add each label name to the data points 
            text(chanlocs(i).X, chanlocs(i).Y, chanlocs(i).Z, ...
                labels{i}, 'FontSize', 10, 'FontWeight', 'bold', 'HorizontalAlignment', 'left', 'Color', color);
        end
    end
  
    % professional naming of the plot 
    xlabel("X");
    ylabel("Y");
    zlabel("Z");
    title("3D Representation of EEG Channels");
    
    % adjust view angle
    view(3);
    % since brain is not a sphere, we need the relative scale of X, Y, & Z
    % axis equal; -> This makes the 3 axes equal 
        
    grid on;

    % wait until the graph is closed to proceed to the final display 
    waitfor(gcf); 
end

### Channel Highlight Options
### Final Brain Map Display (Interactive)
[Video Demonstration](https://github.com/mixuanpan/Ulty_Multy/blob/main/Media/01/05_Brain_Plot_03_3D_Map_Display_Video.mov)

Final Display Function
----------------------

This function compares the selected subject’s data with other subjects’ data. The selected subject’s data is displayed in transparent red, while the compared subject’s data appears in gray. Overlapping areas are shown in darker red, highlighting similarities between the two datasets. This comparison is performed channel by channel, helping to identify consistent patterns and variations across individuals. To improve efficiency and avoid excessive loading time, the program saves the comparison plots directly to a folder instead of displaying them.

function ip_final_display(subject_num, low, upper, noise_num, freqs, power_spec)
    % compare the target to all of the subjects 
    % Data Packaging 

    % Initialization 
    fprintf("\n<strong>Final Analysis Results:</strong>\n\n");
    
    % There are a total of 27 subjects with 3 missing data 
    subject_count = 27; 
    invalid_num = [14 15 25]; % no subjects found 
    current_subject = 1; % subjects to iterate through, starting from 1
    currentPath = "/Users/mixuan/Desktop/Ulty_Multy"; % data directory 
    
    % results directory 
    resultsPath = "/Users/mixuan/Desktop/Ulty_Multy/Results"; 
    channels_num = 32; % number of channels in the given data 
    
    % ip_brain_plot.m
    channels = readtable(sprintf("data/sub-%d/eeg/sub-%d_task-SemanticCategorization_channels.tsv", subject_num, subject_num), "FileType", "text", ... 
        "Delimiter", "\t");
    labels = channels.name;
    % names of the 32 channels 

    % iterate through each subject to analyze and get the data 
    while current_subject <= subject_count
        % ip_subject_info.m
        % check if the current subject is the target subject 
        current_subject = current_subject + (current_subject == subject_num) * (current_subject < subject_count);
        
        % check if the subject number is invalid 
        if current_subject >= invalid_num(1) && current_subject <= invalid_num(end)
            for i = 1:length(invalid_num)
                current_subject = current_subject + (current_subject == invalid_num(i));
            end
        end
        
        % ip_data_filter 
        % filter each subject data with the bandwith given to the target
        % subject 
        % the file for reading EEG data
        fileName = "sub-" + int2str(current_subject) + "_task-SemanticCategorization_eeg.set"; 

        % the path for the EEG data file 
        filePath = fullfile(currentPath, "data", "sub-" + int2str(current_subject), "eeg"); 
        cd(filePath); % change to file directory 

        % load eeg data list from the current subject 
        EEG = pop_loadset(char(fileName)); 
        cd(currentPath) % back to script directory 

        % bandpass filter 
        if low ~= 0 && upper ~= 0
            % both "low ~= 0" and "upper ~= 0" are numeric scales ->
            % use && instead of & for performance 
            EEG = pop_eegfiltnew(EEG, low, upper); % Bandpass filter in the desired Hz range 
        end

        % noise removal 
        if noise_num ~= 0
            EEG = pop_eegfiltnew(EEG, noise_num - 1, noise_num + 1, [], 1); % notch filter 
        end
        
        % ip_time_frequency_analysis.m 
        fig = figure("Visible", "off"); % Set the figure to be invisible 

        % % 4 rows, 8 columns for 32 subplots
        layout = tiledlayout(4, 8, 'Padding', 'compact', 'TileSpacing', 'compact'); 
        for ch = 1:channels_num
            
            % Wavelet Decomposition for each channel 
            % power data for this specific channel 
            % only need the power data, neglecting the frequency 
            power_data = timefreq(EEG.data(ch, :), EEG.srate, 'freqs', [1, 60]); 
            
            nexttile;
            
            % Compared subject data in gray 
            plot(freqs, abs(power_data), 'Color', [0.7, 0.7, 0.7, 0.03]);
            hold on; % stay on the same figure 

            % Target subject data in red
            plot(freqs, abs(power_spec{ch}), 'Color', [1, 0, 0, 0.006], "LineWidth", 1.5); 

            % Channel name for each subplot 
            title(sprintf("%s", labels{ch}), "FontWeight", "bold", 'Units', 'normalized', 'Position', [1, 1, 1]); 
        end 

        % Graph formatting for all subplots 
        xlabel(layout, "Frequency (Hz)", "FontWeight", "bold"); % y label 
        ylabel(layout, "Power (VA)", "FontWeight", "bold"); % x label 

        % title 
        title(layout, sprintf("Subject %i - Gray vs. Target Subject - Red (%i)", current_subject, subject_num)); 
        
        % save the graph file 
        saveas(fig, fullfile(resultsPath, sprintf("output_%i.png", current_subject)));
        close(fig); % close the current figure 

        fprintf("\n<strong>Final display progress: subject %i</strong>\n\n", current_subject);
        current_subject = current_subject + 1; % change to the next subject 
    end
end

### Plot Type Determination
To determine the best way to plot power (complex data) vs. frequency, I tested two approaches: plotting the magnitude of power on a real plane (left) and the original complex data on a complex plane (right). The left graph clearly shows variations in power magnitude over a frequency range, while the right graph is obscure and not useful for analysis. Therefore, I chose to plot the magnitude of power (amplitude) in the final display.

### Example Display Graphs For Subject 27
Subject 1 vs. Subject 27

Subject 3 vs. Subject 27

Subject 22 vs. Subject 27

Subject 24 vs. Subject 27

----------------------
[Published with MATLAB® R2024b](https://www.mathworks.com/products/matlab/)