%This script is licensed under Creative Commons Attribution Non-Commercial (CC-BY-NC)
%Copyright 2019 Caltech

%Instructions for running this script:
%Put this script with Inc.Data.*.mat and IncSham.Data.*.mat files into the same folder
%Type analyzeIncTest in the MATLAB command window or click Run

%What the data files contain:
%data for inclination rotation experiment with matched shams (Table 1 of paper)
%eeg is raw continuous recording (channels x samples)
%trg is triggers with condition labels (condition, sample)
%run is information about the experiment

%What this script does:
%Computes time-freq power for each trial of each participant and does trial rejection
%Computes time-freq power for each condition of each participant
%Computes average alpha-erd for each condition of each participant
%Does two-way repeated-measures anova for interaction of inclination
%rotation (sweep, fixed) and magnetic stimulation (active, sham)

%Outputs:
%ANOVA table figure
%Bar graph figure
%Summary table and additional statistics variables saved in stat.IncTest.mat file

%Important note:
%Due to subsampling trials to balance trial numbers across conditions for
%some participants, the final results will vary slightly for each run of
%this script. The conclusions are not affected by slight changes in
%p-value.

%Note:
%The functions in this script work in MATLAB 2018a. If you are using an
%older version of MATLAB and get any compatibility errors, please update
%MATLAB.

clear all;
for n=1:2
    %get data files
    if n==1
        datafile=dir('Inc.Data.*.mat');
    elseif n==2
        datafile=dir('IncSham.Data.*.mat');
    end
    %set up NaN to be filled with average alpha-erd for each condition of each participant
    %29 participants x 2 conditions (sweep, fixed)
    nsub=length(datafile);
    val=nan(nsub,2);
    %loop through participants
    for isub=1:nsub
        eeg=[]; trg=[]; run=[];
        trl=[]; erd=[]; perm={};
        eegtrl=[]; pwrtrl=[];
        eegavg=[]; pwravg=[];
        %load participant data
        filename=datafile(isub).name;
        disp(filename);
        load(filename);
        %set up information about the experiment
        expt.srate=512;
        expt.ntrl=run.ntrl;
        expt.ttrl=run.ttrl;
        expt.nsamp=expt.ttrl*expt.srate;
        expt.tsamp=[-0.75 expt.ttrl-1+0.25];
        expt.tbase=[-0.5 -0.25];
        expt.fsamp=[1 100];
        %set up time-freq window constraints for computing alpha-erd
        time=[0 expt.ttrl-1];   %post-stimulus interval
        freq=[6 14];            %extended alpha band
        win=[0.125 2];          %window size is 250 ms and 5 Hz bandwidth
        %segment eeg into trials
        psamp=[floor(expt.tsamp(1)*run.srate) floor(expt.tsamp(2)*run.srate)-1];
        samp=repmat(cell2mat(trg(:,2)),[1 length(psamp)])+repmat(psamp,[run.ntrl 1]);
        pbase=[floor(expt.tbase(1)*run.srate) floor(expt.tbase(2)*run.srate)-1];
        base=repmat(cell2mat(trg(:,2)),[1 length(pbase)])+repmat(pbase,[run.ntrl 1]);
        for itrl=1:run.ntrl
            datasamp=eeg(:,samp(itrl,1):samp(itrl,2));
            database=eeg(:,base(itrl,1):base(itrl,2));
            trl(:,:,itrl)=datasamp-repmat(mean(database,2),[1 run.ttrl*run.srate]);
        end
        if run.srate>expt.srate
            trl=trl(:,1:run.srate/expt.srate:end,:);
        end
        srate=expt.srate;
        nsamp=expt.nsamp;
        %select trials of each condition
        trlpos=find(run.stim==1);   %down
        trlneg=find(run.stim==-1);  %up
        trloff=find(run.stim==0);   %fixed
        %subtract erp to compute non-phase-locked power later
        eegtrl=trl;
        eegtrl(:,:,trlpos)=trl(:,:,trlpos)-mean(trl(:,:,trlpos),3);
        eegtrl(:,:,trlneg)=trl(:,:,trlneg)-mean(trl(:,:,trlneg),3);
        eegtrl(:,:,trloff)=trl(:,:,trloff)-mean(trl(:,:,trloff),3);
        %compute power in alpha band for each trial
        fdata=[expt.fsamp(1):expt.fsamp(2)];
        fstd=linspace(3,9,length(fdata))./(2*pi*fdata);
        twave=[-1:1/srate:1];
        nwave=length(twave);
        hwave=(nwave-1)/2;
        ndata=nsamp*expt.ntrl;
        nconv=ndata+nwave-1;
        ffteeg=fft(reshape(eegtrl(1,:,:),1,ndata),nconv);
        for ifreq=freq(1):freq(2)
            wave=exp(1i*2*pi*fdata(ifreq).*twave).*exp((-twave.^2)./(2*fstd(ifreq)^2));
            fftwave=fft(wave,nconv);
            eegwave=ifft(ffteeg.*fftwave);
            eegwave=eegwave(hwave+1:end-hwave);
            eegwave=reshape(eegwave,nsamp,expt.ntrl);
            pwrtrl(:,ifreq-freq(1)+1,:)=eegwave.*conj(eegwave);
        end
        pwrtrl=log10(squeeze(mean(pwrtrl,2)));
        %collect trials with too high or low alpha power in pre- and post-stimulus intervals
        disp([filename ' trial rejection']);
        etrl=[];
        thr=1.5;
        pbase=[floor(diff([expt.tsamp(1) expt.tbase(1)])*srate)+1:floor(diff([expt.tsamp(1) expt.tbase(2)])*srate)];
        test=mean(pwrtrl(pbase,:),1);
        q1=quantile(test,0.25);
        q3=quantile(test,0.75);
        iqr=q3-q1;
        lo=q1-thr*iqr;
        hi=q3+thr*iqr;
        etrl=[etrl find(test<lo | test>hi)];
        pdata=[floor(diff([expt.tsamp(1) 0])*srate)+1:floor(diff([expt.tsamp(1) 1])*srate)];
        test=mean(pwrtrl(pdata,:),1);
        q1=quantile(test,0.25);
        q3=quantile(test,0.75);
        iqr=q3-q1;
        lo=q1-thr*iqr;
        hi=q3+thr*iqr;
        etrl=[etrl find(test<lo | test>hi)];
        if run.ttrl==3
            pdata=[floor(diff([expt.tsamp(1) 1])*srate)+1:floor(diff([expt.tsamp(1) 2])*srate)];
            test=mean(pwrtrl(pdata,:),1);
            q1=quantile(test,0.25);
            q3=quantile(test,0.75);
            iqr=q3-q1;
            lo=q1-thr*iqr;
            hi=q3+thr*iqr;
            etrl=[etrl find(test<lo | test>hi)];
        end
        etrl=unique(etrl)
        if isempty(etrl)
            etrl=[];
        end
        %downsample for speed over iterations (recommended to avoid computer problems)
        %down=1 is the full sampling rate of 512 Hz used for recording eeg
        %down=4 is 128 Hz
        down=4;
        srate=expt.srate/down;
        nsamp=expt.nsamp/down;
        trl=trl(:,1:down:end,:);
        %compute power for each condition
        fdata=[expt.fsamp(1):expt.fsamp(2)];
        fstd=linspace(3,9,length(fdata))./(2*pi*fdata);
        twave=[-1:1/srate:1];
        nwave=length(twave);
        hwave=(nwave-1)/2;
        ndata=nsamp*expt.ntrl;
        nconv=ndata+nwave-1;
        %subsample trials to balance trial numbers across conditions
        if length(trloff)<floor(expt.ntrl/2)
            nperm=100;
        else
            nperm=1;
        end
        %loop through subsamples
        for iperm=1:nperm
            trlpos=find(run.stim==1);
            trlneg=find(run.stim==-1);
            trloff=find(run.stim==0);
            %because sweep trials are overrepresented in some participants,
            %subsample the same number of sweep trials as fixed trials
            if length(trloff)<floor(expt.ntrl/2)
                trlpos=sort(randsample(trlpos,length(trlpos)/2));
                trlneg=sort(randsample(trlneg,length(trlneg)/2));
            end
            trlon=sort([trlpos; trlneg]);
            trlall=sort([trlon; trloff]);
            %exclude rejected trials from further calculations
            trlpos=trlpos(~ismember(trlpos,etrl));
            trlneg=trlneg(~ismember(trlneg,etrl));
            trloff=trloff(~ismember(trloff,etrl));
            trlon=trlon(~ismember(trlon,etrl));
            trlall=trlall(~ismember(trlall,etrl));
            %subtract erp to compute non-phase-locked power later
            eegavg=trl;
            eegavg(:,:,trlpos)=trl(:,:,trlpos)-mean(trl(:,:,trlpos),3);
            eegavg(:,:,trlneg)=trl(:,:,trlneg)-mean(trl(:,:,trlneg),3);
            eegavg(:,:,trloff)=trl(:,:,trloff)-mean(trl(:,:,trloff),3);
            %select midline frontal channel Fz (channel 38 in Biosemi 64 layout)
            chan=38;
            disp([filename ' chan ' num2str(chan) ' perm ' num2str(iperm)]);
            ffteeg=fft(reshape(eegavg(chan,:,:),1,ndata),nconv);
            for ifreq=freq(1):freq(2)
                wave=exp(1i*2*pi*fdata(ifreq).*twave).*exp((-twave.^2)./(2*fstd(ifreq)^2));
                fftwave=fft(wave,nconv);
                eegwave=ifft(ffteeg.*fftwave);
                eegwave=eegwave(hwave+1:end-hwave);
                eegwave=reshape(eegwave,nsamp,expt.ntrl);
                %compute power for each condition (sweep, fixed, all trials)
                pwravg(:,ifreq-freq(1)+1,1)=mean(eegwave(:,trlon).*conj(eegwave(:,trlon)),2);
                pwravg(:,ifreq-freq(1)+1,2)=mean(eegwave(:,trloff).*conj(eegwave(:,trloff)),2);
                pwravg(:,ifreq-freq(1)+1,3)=mean(eegwave(:,trlall).*conj(eegwave(:,trlall)),2);
            end
            %save condition averages for subsample
            perm(iperm,1)={pwravg};
        end
        %select time-freq window
        pdata=[floor(diff([expt.tsamp(1) time(1)])*srate)+1:floor(diff([expt.tsamp(1) time(2)])*srate)];
        pbase=[floor(diff([expt.tsamp(1) expt.tbase(1)])*srate)+1:floor(diff([expt.tsamp(1) expt.tbase(2)])*srate)];
        twin=floor(win(1)*srate);
        tmin=twin+1;
        tmax=length(pdata)-twin+1;
        fwin=win(2);
        fmin=fwin+1;
        fmax=size(pwravg,2)-fwin;
        for iperm=1:nperm
            pwr=perm{iperm};
            pwr0=pwr(:,:,3);    %all trials average
            pwr1=pwr(:,:,1);    %sweep trials average
            pwr2=pwr(:,:,2);    %fixed trials average
            pwr0=10*log10(pwr0./repmat(mean(pwr0(pbase,:),1),[size(pwr0,1) 1]));
            pwr1=10*log10(pwr1./repmat(mean(pwr1(pbase,:),1),[size(pwr1,1) 1]));
            pwr2=10*log10(pwr2./repmat(mean(pwr2(pbase,:),1),[size(pwr2,1) 1]));
            pwr0=squeeze(pwr0(pdata,:));
            pwr1=squeeze(pwr1(pdata,:));
            pwr2=squeeze(pwr2(pdata,:));
            %select most negative time-freq point in all trials average
            %center time-freq window over this point
            [tpnt,fpnt]=ind2sub(size(pwr0),find(pwr0==min(min(pwr0))));
            %select time window
            if tpnt<tmin
                tpnt=tmin;
            elseif tpnt>tmax
                tpnt=tmax;
            end
            trng=tpnt+twin*[-1 1]+[0 -1];
            %select frequency window
            if fpnt<fmin
                fpnt=fmin;
            elseif fpnt>fmax
                fpnt=fmax;
            end
            frng=fpnt+fwin*[-1 1];
            %compute average alpha-erd in time-freq window for each condition
            erd(iperm,1)=mean(mean(pwr1(trng(1):trng(end),frng(1):frng(end))));
            erd(iperm,2)=mean(mean(pwr2(trng(1):trng(end),frng(1):frng(end))));
        end
        %get alpha-erd average for each condition of each participant
        val(isub,:)=mean(erd,1);
    end
    tot(n,1)={val};
end
%two-way repeated-measures anova
val1=tot{1};
val2=tot{2};
nsub1=size(val1,1);
nsub2=size(val2,1);
ncond1=size(val1,2);
ncond2=size(val2,2);
delta=[val1(:,1); val1(:,2); val2(:,1); val2(:,2)];
g1=[repmat({'sweep'},[nsub1 1]); repmat({'fixed'},[nsub1 1]); repmat({'sweep'},[nsub2 1]); repmat({'fixed'},[nsub2 1])];
g2=[repmat({'stim'},[nsub1*ncond1 1]); repmat({'sham'},[nsub2*ncond2 1])];
g3=[repmat([1:nsub1]',[ncond1 1]); repmat([1:nsub2]',[ncond2 1])];
%output anova table
[~,tbl,stat]=anovan(delta,{g1 g2 g3},'random',3,'model','interaction','varnames',{'inc' 'stim' 'sub'});
figure('visible','off');
%tukey's test
comp=multcompare(stat,'dimension',[1 2]);
avg=[mean(val1,1) mean(val2,1)];
fval=tbl{5,6};      %interaction f value
pval=tbl{5,7};     %interaction p value
sscond1=tbl{2,2};
sscond2=tbl{3,2};
ssint=tbl{5,2};
sserr=tbl{8,2};
eff=[sscond1/(sscond1+sserr) sscond2/(sscond2+sserr) ssint/(ssint+sserr)];
col1={'ANOVA #1'; 'Interaction of Inclination Rotation x Magnetic Stimulation'};
col2={'F Value of Interaction'; fval};
col3={'P Value of Interaction'; pval};
col4={'Effect Size of Interaction'; eff(3)};
summary=table2array(table(col1,col2,col3,col4));
%save statistics variables
save('stat.IncTest.mat','avg','val','fval','pval','eff','tbl','comp','summary');
%bar graph
fontsize=24;
set(0,'defaultaxesfontsize',fontsize);
set(0,'defaultaxesfontname','arial');
avg1=mean(val1,1);
avg2=mean(val2,1);
sem1=std(val1,0,1)/sqrt(nsub1);
sem2=std(val2,0,1)/sqrt(nsub2);
fig=figure('color',[1 1 1],'position',[50 50 800 600]);
x1=[0.5:0.75:1.5];
x2=[2.5:0.75:3.25];
x=[x1 x2];
y=[avg1 avg2];
for ibar=1:length(y)
    hold on;
    b=bar(x(ibar),y(ibar),'barwidth',0.5);
    if sum(ismember(x([1 3]),x(ibar)))>0
        color=[0 0.4 1];
    elseif sum(ismember(x([2 4]),x(ibar)))>0
        color=0.5*[1 1 1];
    end
    set(b,'facecolor',color,'edgecolor',color);
end
axis([x([1 end])+0.5*[-1 1] 2*[-1 0.25]]);
box on;
ylabel(['Alpha-ERD (dB)']);
xtick=x;
xticklabel=repmat({'SWEEP' 'FIXED'},[1 2]);
ytick=-2:0.5:0;
yticklabel=num2cell(ytick);
set(gca,'xtick',xtick,'xticklabel',xticklabel,'ytick',ytick,'yticklabel',yticklabel,'fontsize',fontsize);
hold on;
text(mean(x1),ytick(end)+0.25,'IncActive','fontsize',fontsize,'fontweight','bold','horizontalalignment','center');
text(mean(x2),ytick(end)+0.25,'IncSham','fontsize',fontsize,'fontweight','bold','horizontalalignment','center');
text(2.5,min(ytick)+0.2,['N=' num2str(nsub1)],'fontsize',fontsize,'fontweight','bold','horizontalalignment','center');
err=[sem1 sem2];
cap=0.05;
for i=1:length(err)
    hold on;
    plot(x(i)*ones([1 2]),y(i)+err(i)*[-1 1],'color',0.1*[1 1 1],'linewidth',1);
    hold on;
    plot(x(i)+cap*[-1 1],y(i)+err(i)*[-1 -1],'color',0.1*[1 1 1],'linewidth',1);
    hold on;
    plot(x(i)+cap*[-1 1],y(i)+err(i)*[1 1],'color',0.1*[1 1 1],'linewidth',1);
end
text(2.5,min(ytick)+0.4,['p=' num2str(pval(1)) ' (Interaction)'],'color','m','fontsize',fontsize,'fontweight','bold','horizontalalignment','center');
clear all;
load('stat.IncTest.mat');