EEG_FAMvs.IMW/hem_powerBand_pwelchHamming1sEpochs.m

%% Compute average power in band for HEM ICs during the 2x2 TaskXCond using pwelch() on 1" epochs tapered with Hanning window

% Task labels
task = {'MOM-INS', 'MOM-SEP', 'MTT-INS', 'MTT-SEP'};

% bands
alpha = [8 12];
beta = [12 25];
gamma = [35 45];
delta = [1 4];
theta = [4 8];

% band_limits = [8 12; 12 25; 35 45; 1 4; 4 8];
band_limits = [alpha; beta; gamma; delta; theta];
band_names = {'alpha', 'beta', 'gamma', 'delta', 'theta'};

% sampling frequency
Fs= 250;


% Load signals for HEM ICs (n.tasks, n.subjects, n.points)
load hem_data.mat;
% permute (n.points, n.subjects, n.tasks)
hem_data = permute( hem_data, [3, 2, 1]);

tot_sub=  size( hem_data, 2);
tot_epochs= floor( size( hem_data, 1) / Fs);

% for each task
for n_task= 1: 4
    % for each subject
    for n_sub= 1: tot_sub
        
        % data for current subject, task
        y= hem_data(:, n_sub, n_task);
        % for each 1s-epoch
        for t= 1: tot_epochs
            % 1s-epoch
            x = y(1+(t-1)*250 : t*250);

            N = length(x);
            freq = 0:Fs/length(x):Fs/2;
            % power spectral density (PSD) estimate using welch estimator
            [pxx, f] = pwelch(x, 250, 0, freq, Fs);

            psd(t,:) = pxx;
        end
        
        % compute average power in each band for each 1s-epoch
        for n_band = 1:length(band_limits)
        
            % lower bound for selected band
            b = band_limits(n_band,1);
            % higher bound for selected band
            B = band_limits(n_band,2);
            % container for frequencies of specdata in selected band
            freq_b = zeros(size(freq));
            % for all frequencies of specdata
            for j = 1:length(freq)
                % extract frequencies
                if and( freq(j)>=b, freq(j)<=B )
                    freq_b(j) = 1;
                end
            end

            % band frequencies matrix for all 1s-epoch
            freq_bR = repmat(freq_b, [size(psd, 1), 1] );
            % data power in selected band
            dp = freq_bR .* psd;

            % extract non-zero columns
            % (https://it.mathworks.com/matlabcentral/answers/40018-delete-zeros-rows-and-columns)
            dp( :, ~any(dp,1) ) = [];

            % power in band for each 1s-epoch
            t= sprintf('Average power in band %s for subject %d in task %s...', band_names{n_band}, n_sub, task{n_task}); disp(t);
            pAvg_inTime_band = mean( dp, 2);
            
            % power in band of whole VEM IC for current subject, task
            pAvg_band( n_task, n_sub, n_band) = mean( pAvg_inTime_band);
        end
        
    end
end


% Save data
% Rearrange matix (n.band, n.subject, n.taskXcond)
data = permute( pAvg_band, [3,2,1] );
save hem_powerBand_pwelch.mat data;

% figure;plot(pAvg_inTime_band');legend('alpha','beta','gamma','delta','theta')
% 
% plot(freq,10*log10(psdx))
% grid on
% title('Periodogram Using FFT')
% xlabel('Frequency (Hz)')
% ylabel('Power/Frequency (dB/Hz)')