%% SSVEP Data Analysis Source Code 
% This source code is an attachment of JoVE Manuscript 62082
% Fong, Cohen, Mahony, et al. (2021)
% The authors recommend the use of MATLAB R2015b and higher
% However source code is applicable to other software platforms

%% Variation in Code for Protocol Sections 8 and 13
% Refer to comments for variations between EDF and CSV files
% Section 8 refers to analysis of EDF files
% Section 13 refers to analysis of CSV files


%% Set Variables for Analysis 

    clc;
    close all;
    clear all;

% Cut off frequency set to 40 Hz
    cutoff = 40; 
    
% Adjust Sampling Frequency (fs) to match EEG device used
% Authors of JoVE 62082 set fs to be 250 Hz for both devices
    fs = 250;


%% Input File

% Pop up function, select specific file for analysis
[file,path] = uigetfile({'*.CSV;*.EDF;*.JSON','Data Files';'*.EDF','*.EDF';'*.CSV','*.CSV';'*.JSON','*.JSON'},'Select file to open');
workingDir = pwd;
 
 if contains(file,'.csv') || contains(file,'.CSV')
        cd(path);
        fileType = 'csv';
        data = transpose(csvread(file));
      
 elseif contains(file,'.edf') || contains(file,'.EDF')
        fileType = 'edf';
        [header, data] = readEdf(file, path);
        file = strrep(file,'.edf','.csv');
        fs = header.samples(1)/header.duration;
        data = transpose(data);
 end

        cd(workingDir);


%% Filter Data
% Butterworth Bandpass filter used
% Corner freqieunces of 5 and 40 Hz set
    [b,a] = butter(3, [5/(fs/2) 40/(fs/2)], 'bandpass');

    [~,length_data] = size(data);
    
    % For section 13, use commented section below to process 3rd channel
    % replace below with; filteredData = zeros(3,length_data);
    filteredData = zeros(2,length_data);

    for channels = 1:2
        filteredData(channels,:) = filtfilt(b,a,data(channels, :));
    end
    
%% Fast Fourier Transformation (FFT)

% Padding for Optimization
    p = nextpow2(length_data)+2;                                      
    n = 2^p;                                                             

    power = zeros(channels,n);                                          

    % For section 13, use commented section below to process 3rd channel
    % replace below with; channel = 1:3
    for channel = 1:2
        power(channel,:) = abs(fft(filteredData(channel,:),n))/length(data);
    end

%% Normalise Data (Averages)

    averageNorm(1) = mean(power(1,:));
    averageNorm(2) = mean(power(2,:));
    
    % For section 13, use commented section below to process 3rd channel
    % averageNorm(3) = mean(power(3,:));
    
    average_total =  mean(averageNorm);

    power(1,:) = power(1,:) / averageNorm(1) * average_total;            
    power(2,:) = power(2,:) / averageNorm(2) * average_total;
    
    % For section 13, use commented section below to process 3rd channel
    % power(3,:) = power(3,:) / averageNorm(3) * average_total;
    
    % Multiply by 1000000 to correct NUROCHEK data to uV
    if fs == 250
        power = power * 1000000;
    end
    

    % For section 13, use commented section below to process 3rd channel
    % replace below with; power_comb = power(1,:) + power(2,:) + power(3,:);
    power_comb = power(1,:) + power(2,:);
    freq = 0 : fs/n : fs - fs/n;

    x = freq(1:int64(cutoff/fs*n));
    y1 = power(1,1:int64(cutoff/fs*n));
    y2 = power(2,1:int64(cutoff/fs*n));
    
    % For section 13, use commented section below to process 3rd channel
    % y3 = power(3,1:int64(cutoff/fs*n));
    
    y = power_comb(1:int64(cutoff/fs*n));

    cd(path);

    average_axis(1:int64(cutoff/fs*n)) = mean(power_comb(int64(5/fs*n):int64(40/fs*n)));


%% Plotting

    figure('Name', 'Combined Electrodes');
    xlabel('Frequency (Hz)');
    ylabel('Magnitude (uV)');

    plot(x, y, 'r');
    hold on;
    plot(x, average_axis, 'k--');
    legend('Combined', 'Average');
    legend boxoff 
    

 %% Sign to Noise Calculation    
    
    average_overall = mean(power_comb(int64(5/fs*n):int64(40/fs*n)));

    [peakMag,index] = max(power_comb(1:int64(40/fs*n)));

    peakFreq = freq(index);
    stdevBands = std(power_comb(1,int64(5/fs*n):int64(40/fs*n)));
    peakRating = (peakMag-average_overall)/stdevBands ;
    snr = peakMag/average_overall;
 
    
 %% End of Source Code
 
    
    
    
    
    