load('data.mat');

%

size(data_eeg) % [ sample x channel ]

%

clab_eeg

%

figure
plot(data_eeg(1e5 + (1:5:50000),:)) % 10 second of data, all channels

%

figure
plot(data_emg(1e5 + (1:5:50000),1)) % 10 second of data, all channels

%

figure
plot(detrend(data_eeg(1e5 + (1:5:50000),:), 0)) % de-mean

%

% Create epochs
EEG_EPOCH_LIMITS = [-1.000 -0.005]; %s
EMG_EPOCH_LIMITS = [-0.500  0.500]; %s
EMG_BASELINE_WIN = [-0.100 -0.005]; %s
epochs_eeg.trial = zeros(diff(EEG_EPOCH_LIMITS*fs)+1, length(ix), size(data_eeg, 2));
epochs_eeg.dimord = 'time_rpt_chan';
epochs_eeg.time = (EEG_EPOCH_LIMITS(1)*fs:EEG_EPOCH_LIMITS(2)*fs)/fs;
epochs_eeg.clab = clab_eeg;
epochs_emg = [];
epochs_emg.trial = zeros(diff(EMG_EPOCH_LIMITS*fs)+1, length(ix), size(data_emg, 2));
epochs_emg.dimord = 'time_rpt_chan';
epochs_emg.time = (EMG_EPOCH_LIMITS(1)*fs:EMG_EPOCH_LIMITS(2)*fs)/fs;
epochs_emg.clab = clab_emg;
for e = 1:length(ix) % [time, trial, channel]
    epochs_eeg.trial(:,e,:) = data_eeg(ix(e)+(EEG_EPOCH_LIMITS(1)*fs:EEG_EPOCH_LIMITS(2)*fs), :);
    epochs_emg.trial(:,e,:) = data_emg(ix(e)+(EMG_EPOCH_LIMITS(1)*fs:EMG_EPOCH_LIMITS(2)*fs), :);
end
epochs_eeg.trial = detrend(epochs_eeg.trial, 0); % de-mean
epochs_emg.trial = epochs_emg.trial - mean(epochs_emg.trial(epochs_emg.time > EMG_BASELINE_WIN(1) & epochs_emg.time < EMG_BASELINE_WIN(2),:,:), 1); % baseline correction

%

epochs_eeg

%

% downsample eeg data to 1 kHz

%

figure
plot(squeeze(mean(epochs_eeg.trial, 2)))

%

figure
plot(squeeze(mean(epochs_eeg.trial(:,:,ismember(epochs_eeg.clab, {'C3', 'CP1', 'CP5', 'FC1', 'FC5'})), 2)))

%

% Spectral Estimation using PSD multi-taper method
fs = 1/mean(diff(epochs_eeg.time));
ix = find(ismember(epochs_eeg.clab, {'C3', 'CP1', 'CP5', 'FC1', 'FC5'}))';
figure; ax = axes; hold on;
for c = ix
    [Pxx, F] = pmtm(epochs_eeg.trial(:,:,c), 4/2, 2^nextpow2(1.5 * size(epochs_eeg.trial, 1)), fs);
    plot(F, pow2db(mean(Pxx, 2)))
end
xlabel('Frequency (Hz)'),
ylabel('Power (dB)')
xlim([6 120])
ax.XScale = 'log';
ax.XTick = [5 10 15 30 50 100];
legend(ax, epochs_eeg.clab(ix))
box on; grid on

%

% Apply spatial filter montage
SPF_CHANNELS = {'C3', 'CP1', 'CP5', 'FC1', 'FC5'};
SPF_WEIGHTS = [1 -0.25 -0.25 -0.25 -0.25];
epochs_spf = [];
epochs_spf.trial = reshape(epochs_eeg.trial(:,:,ismember(epochs_eeg.clab, SPF_CHANNELS)), [], numel(SPF_CHANNELS)) * SPF_WEIGHTS';
epochs_spf.trial = reshape(epochs_spf.trial, size(epochs_eeg.trial, 1), size(epochs_eeg.trial, 2), []);
epochs_spf.trial = detrend(epochs_spf.trial, 0); % de-mean (probably not necessary)
epochs_spf.dimord = 'time_rpt_chan';
epochs_spf.time = epochs_eeg.time;
epochs_spf.clab = {'C3 Hjorth'};

%

epochs_spf

%

figure
plot(epochs_spf.time, epochs_spf.trial(:,1:5,1))

%

% Spectral Estimation using PSD multi-taper method
fs = 1/mean(diff(epochs_eeg.time));
figure; ax = axes;
[Pxx, F] = pmtm(epochs_eeg.trial(:,:,1), 4/2, 2^nextpow2(1.5 * size(epochs_eeg.trial, 1)), fs);
plot(F, pow2db(mean(Pxx, 2)))
xlabel('Frequency (Hz)'),
ylabel('Power (dB)')
xlim([6 120])
ax.XScale = 'log';
ax.XTick = [5 10 15 30 50 100];
legend(ax, epochs_spf.clab)
box on; grid on

%

% Determine MEP size
MEP_LIMITS = [0.020 0.040];

figure
plot(epochs_emg.time, squeeze(mean(epochs_emg.trial(:,:,1), 2))), hold on
plot(MEP_LIMITS, [0 0], 'k--')
xlim([0.005 0.05]);
title('average MEP')

epochs_emg.trialinfo = squeeze(range(epochs_emg.trial(epochs_emg.time >= MEP_LIMITS(1) & epochs_emg.time <= MEP_LIMITS(2),:,1)));

%

figure
plot(epochs_emg.trialinfo, '.')

%

y = log10(epochs_emg.trialinfo);

figure
plot(y, '.')
hold on
line([0 size(y, 2)], [1 1] .* median(y), 'Color', 'red', 'LineWidth', 1.5)
title('Response amplitude over time')

%

figure
plot(mean(epochs_spf.trial(:,y > median(y)), 2))
hold on
plot(mean(epochs_spf.trial(:,y < median(y)), 2))