% clear the memory
clear;
% close all open files and windows
close all;
% clear the command console
clc;

% loop through and load the participants
for i = 1:10

    % specify a file name
    fileName = ['Subject Wavelet ' num2str(i) '.mat'];

    % load the file
    load(fileName);

    % put the data into a variable called wavData
    wavData(:,:,:,:,i) = WAV.data;

end

% now with wavelets to avoid edge artifacts we are going to trim the first
% 100 ms and the last 100 ms
wavData(:,:,1:50,:,:) = [];
wavData(:,:,end-49:end,:,:) = [];

% take the grand average
grandWAV = mean(wavData,5);

% plot the Wavelets for both conditions side by size for channel 52
condition1WaveletData = squeeze(grandWAV(52,:,:,1));
condition2WaveletData = squeeze(grandWAV(52,:,:,2));
% now the plotting is more complex as we are working with an image with 3D
% data
subplot(1,2,1);
imagesc(condition1WaveletData);
title('Channel Pz: Condition One');
xlabel('Samples');
ylabel('Frequency (Hz)');
set(gca,'YDir','normal');
subplot(1,2,2);
imagesc(condition2WaveletData);
title('Channel Pz: Condition Two');
xlabel('Samples');
set(gca,'YDir','normal');
ylabel('Frequency (Hz)');

% compute the difference scores
differenceData = squeeze(wavData(52,:,:,1,:) - wavData(52,:,:,2,:));

% and generate a statistical map, you can just alpha using the Bonferroni
% rules.
alpha = 0.05;
statMap = doWAVStats(differenceData,alpha);

% and lets plot the result
figure;
imagesc(statMap);
title('Channel Pz: Difference');
xlabel('Samples');
ylabel('Frequency (Hz)');
set(gca,'YDir','normal');