|
|
@@ -1,98 +0,0 @@
|
|
|
-%% 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)')
|
