%Function to Accompany Chapter 6 of Dickter & Kieffaber (2013) EEG Methods

%DK_PSD

%Function to accompany Chapter 6 of Dickter & Kieffaber (2013) EEG Methods

%for Social and Personality Psychology

%Usage:

%> [EEG ALLEEG CURRENTSET com]=DK_PSD(EEG,ALLEEG,CURRENTSET,varargin)

%

% Inputs:

% 1) EEG - EEGlab data structure

%

% Optional inputs:

% 'baseline' - [start stop] Start and End latencies (ms) to use as baseline.

% 'datawindow' - [start stop] Start and End latencies (ms) to use analyze

%

%Output:

% EEG - Abominized EEG dataset converted to the power spectrum

% (uV^2) of the input data set. The form of the new data is Channels X

% Frequencies X Trials. CAUTION: EEG.times and EEG.srate fields modified

% accordingly in order to make available most visualization tools.

%

% ------

% AUTHOR:

% Paul Kieffaber

% Department of Psychology &

% Program in Neuroscience

% College of William & Mary

%

% ------

%

% THIS PROGRAM IS FREE SOFTWARE INTENDED FOR EDUCATIONAL PURPOSES ONLY.

%

% This program is distributed in the hope that it will be useful,

% but WITHOUT ANY WARRANTY; without even the implied warranty of

% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

function [EEG ALLEEG CURRENTSET com]=DK_PSD(EEG,ALLEEG,CURRENTSET,varargin)

datawindow=[EEG.xmin*1000 EEG.xmax*1000]; %default = whole epoch

%Parse optional inputs

for i=1:2:length(varargin)

arg=lower(varargin{i});

val=varargin{i+1};

switch arg

case 'baseline'

if length(val)==2

baseline=val;

else

disp(['baseline argument must be two-element vector [start end]']);

end

case 'datawindow'

if length(val)==2

datawindow=val;

else

disp(['specification of DataWin requires two-element vector [start end]']);

end

end

end

datastart=find(abs(EEG.times-datawindow(1))==min(abs(EEG.times-datawindow(1))));

datastop=find(abs(EEG.times-datawindow(2))==min(abs(EEG.times-datawindow(2))));

%if baseline specified, then use same nfft for both

if exist('baseline')==1

baselinestart=find(abs(EEG.times-baseline(1))==min(abs(EEG.times-baseline(1))));

baselinestop=find(abs(EEG.times-baseline(2))==min(abs(EEG.times-baseline(2))));

nfft=min([2^nextpow2(baselinestop-baselinestart) 2^nextpow2(datastop-datastart)]);

end

%Compute signal PSD

Pxxdata=[];

for i=1:size(EEG.data,3)

for j=1:size(EEG.data,1)

if exist('nfft')==1

[Pxxdata(j,:,i),F] = pwelch(EEG.data(j,datastart:datastop,i),[],[],nfft,EEG.srate);

else

[Pxxdata(j,:,i),F] = pwelch(EEG.data(j,datastart:datastop,i),[],[],[],EEG.srate);

end

end

end

%Compute baseline PSD

if exist('baseline')==1

Pxxbaseline=[];

for i=1:size(EEG.data,3)

for j=1:size(EEG.data,1)

[Pxxbaseline(j,:,i),F] = pwelch(EEG.data(j,baselinestart:baselinestop,i),[],[],nfft,EEG.srate);

end

end

Pxxdata=10*log10(Pxxdata./Pxxbaseline);

end

%Abominize EEG structure in order to make visualization tools available

disp('Frequencies will be stored in ''EEG.times''');

EEG.data=Pxxdata;

EEG.pnts=length(F);

EEG.times=F;

EEG.setname=strcat('PX_',EEG.setname);

EEG.xmin=0;

EEG.xmax=max(F)/1000;

EEG.srate=(length(F)-1)/(max(F)/1000);

[EEG]=eeg_checkset(EEG);

[ALLEEG EEG CURRENTSET] = pop_newset( ALLEEG, EEG, CURRENTSET,'setname',strcat('PX',EEG.setname),'overwrite','off');

eeglab redraw

com=['[EEG LASTCOM]=pop_cpl_powerspectrum(EEG);'];

disp('Spectral decomposition complete....')

disp('WARNING: EEG structure has been modified. Use caution.');