%calculate artifacts based on various manual and automatic detection
%methods using the Budapest-Munich sample

clear all
% filepath='E:\eeg\Nagy IQ-EEG adatbázis\'; %Budapest computer, HDD, not really good path
% filepath='F:\Nagy IQ-EEG adatbázis\'; %Fehérvár computer, external HDD
% filepath='F:\Nagy IQ-EEG adatbázis\Csak Mensa tagság ismert\';
filepath='F:\Gyerekalvás\gyerekalvás EEG\';
% filepath='F:\Nap study\Control condition\';
% filepath='D:\Sleep database\';
% filepath='F:\Nagy IQ-EEG adatbázis\Csak orsóelemzéshez (rossz felvétel vagy 16 csatorna)\';
addpath('C:\Users\Ujma Péter\OneDrive\MATLAB scripts\')

filenames=dir(fullfile(filepath,'*.edf'));
for k=1:length(filenames)
   edfs{k}=filenames(k).name;
end
edfs=edfs';
bins=0:0.1:48; %bins of interest for PSD analysis
windowlength=4; %length of PSD window in seconds
hyplength=20; %setting hypno epoch duration 
target_electrodes={'C3' 'C4' 'Fp1' 'Fp2' 'F3' 'F4' 'P3' 'P4' 'O1' 'O2'}; %all recordings should have these and maybe they are enough for this analysis
params=[1 1.5 2 2.5 3 4]; %hjorth parameter cutoffs specified as thresholds
params2=[1 1.5 2 2.5 3]; %cutoffs for two-threshold detections

%assigning names
artnames={};
for name=1:length(params)
    artnames{name}=['arts_hjorth_sd' num2str(params(name))];
    name=name+1;
end
artnames2={};
for name=1:length(params2)
    artnames2{name}=['arts_hjorth_sd' num2str(params2(name)) '_' num2str(params2(name))];
    name=name+1;
end
artnames3=['Visual rejection' artnames artnames2 'Buckelmuller' 'No artifact rejection'];

for o=1:length(edfs)
    disp(['Now working with ' edfs{o} ', file ' num2str(o) ' out of ' num2str(length(edfs))])
    %getting EEG
    EEG=pop_biosig([filepath edfs{o}],'importevent','off','importannot','off'); %reading EEG data
    electrodes={EEG.chanlocs.labels};
    clear hypnodata
    clear hypnoepochs
    if strcmp(filepath,'D:\Sleep database\') %different name format for additional files
    hypnogram=textread(['D:\Sleep database\hyps\hyps\' strrep(edfs{o},'.edf','.hyp.txt')]);
    artifacts=textread(['D:\Sleep database\arts\' strrep(edfs{o},'.edf','.art.txt')]);
    else
    fid=fopen([filepath strrep(edfs{o},'.edf','.hyp')]);
    hypnoepochs=fread(fid);
    fclose(fid);
    x=1:7;
    k=1;
    for i=311:size(hypnoepochs)
    if sum(hypnoepochs(i)==x)>0
        hypnodata(k)=hypnoepochs(i);
        k=k+1;
    end
    end
    hypnodata(length(hypnodata)+1:length(EEG.data)/EEG.srate/20+1)=1; %if the hyp file is too short, the rest is marked as junk
    hypnogram=hypnodata;
    %getting artifacts - MANUAL METHOD
    fid=fopen([filepath strrep(edfs{o},'.edf','.art')]);
    artifacts=fread(fid);
    artifacts(1:184)=[]; % the first 184 bits are the header
    artifacts(length(artifacts)+1:length(EEG.data)/EEG.srate/4+1)=2; %if the art file is too short, the rest is marked as artefactual
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%HJÖRTH METHOD%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    %getting NREM/REM epoch labels
    starttime=1;
    endtime=round(starttime+windowlength*EEG.srate);
    nremwindow=1;
    remwindow=1;
    nremtime=[];
    remtime=[];
    clear artifacts_manual_nrem
    clear artifacts_manual_rem
    while endtime<size(EEG.data,2)-hyplength*EEG.srate %trimming off the last hypnogram epoch to be sure
        if sum(hypnogram(ceil(starttime/EEG.srate/hyplength))==[4 5 6])>0 && sum(hypnogram(ceil(endtime/EEG.srate/hyplength))==[4 5 6])>0 
        nremtime(nremwindow)=starttime;
        %assigning manual artifacts to epochs
        if artifacts(ceil(starttime/EEG.srate/4))==1 && artifacts(ceil(endtime/EEG.srate/4))==1
        artifacts_manual_nrem(nremwindow)=1;
        else
        artifacts_manual_nrem(nremwindow)=0;    
        end
        nremwindow=nremwindow+1;
        elseif sum(hypnogram(ceil(starttime/EEG.srate/hyplength))==2)>0 && sum(hypnogram(ceil(endtime/EEG.srate/hyplength))==2)>0 
        remtime(remwindow)=starttime;
        %assigning manual artifacts to epochs
        if artifacts(ceil(starttime/EEG.srate/4))==1 && artifacts(ceil(endtime/EEG.srate/4))==1
        artifacts_manual_rem(remwindow)=1;
        else
        artifacts_manual_rem(remwindow)=0;    
        end
        remwindow=remwindow+1;
        end
        starttime=round(starttime+(windowlength/2)*EEG.srate)+1;
        endtime=starttime+windowlength*EEG.srate;
    end    
 
    
    
 clear spectrum
 clear hjorth_rem
 clear hjorth_nrem
 %getting spectra and Hjorth parameters from all epochs
 for ch=1:length(target_electrodes) %looping through channels
     channel=find(contains(electrodes,target_electrodes{ch}));
     channel=channel(1); %multiple matches are at the end of the order
     for nremepoch=1:nremwindow-1
        signal=EEG.data(channel,nremtime(nremepoch):round(nremtime(nremepoch)+windowlength*EEG.srate));
        [hjorth_nrem.a(ch,nremepoch),hjorth_nrem.m(ch,nremepoch),hjorth_nrem.c(ch,nremepoch)]=hjorth(signal); %Hjorth parameters
        spectrum.nrem(ch,nremepoch,:)=periodogram(signal,hamming(length(signal)),bins,EEG.srate);
     end
     for remepoch=1:remwindow-1
        signal=EEG.data(channel,remtime(remepoch):round(remtime(remepoch)+windowlength*EEG.srate));
        [hjorth_rem.a(ch,remepoch),hjorth_rem.m(ch,remepoch),hjorth_rem.c(ch,remepoch)]=hjorth(signal); %Hjorth parameters
        spectrum.rem(ch,remepoch,:)=periodogram(signal,hamming(length(signal)),bins,EEG.srate);
     end    
 end
 
 clear artifacts_singleparam_nrem
 clear artifacts_singleparam_rem
 %getting binary artifacts based on Hjörth parameters
 for ch=1:length(target_electrodes)
 for th=1:length(params) %single-threshold method
     hthresh=params(th);
     artifacts_singleparam_nrem{ch}(th,:)=(abs((hjorth_nrem.a(ch,:)-mean(hjorth_nrem.a(ch,:)))/std(hjorth_nrem.a(ch,:)))<hthresh & abs((hjorth_nrem.m(ch,:)-mean(hjorth_nrem.m(ch,:)))/std(hjorth_nrem.m(ch,:)))<hthresh & abs((hjorth_nrem.c(ch,:)-mean(hjorth_nrem.c(ch,:)))/std(hjorth_nrem.c(ch,:)))<hthresh);
     artifacts_singleparam_rem{ch}(th,:)=(abs((hjorth_rem.a(ch,:)-mean(hjorth_rem.a(ch,:)))/std(hjorth_rem.a(ch,:)))<hthresh & abs((hjorth_rem.m(ch,:)-mean(hjorth_rem.m(ch,:)))/std(hjorth_rem.m(ch,:)))<hthresh & abs((hjorth_rem.c(ch,:)-mean(hjorth_rem.c(ch,:)))/std(hjorth_rem.c(ch,:)))<hthresh);
 end
 end
 
 clear artifacts_dualparam_nrem
 clear artifacts_dualparam_rem
 for ch=1:length(target_electrodes)
 for th=1:length(params2) %single-threshold method
     hthresh=params2(th);
     hthresh2=params2(th); %there is an option here to give a different value to the two thresholds (currently the same)
     tempparam=(abs((hjorth_nrem.a(ch,:)-mean(hjorth_nrem.a(ch,:)))/std(hjorth_nrem.a(ch,:)))<hthresh & abs((hjorth_nrem.m(ch,:)-mean(hjorth_nrem.m(ch,:)))/std(hjorth_nrem.m(ch,:)))<hthresh & abs((hjorth_nrem.c(ch,:)-mean(hjorth_nrem.c(ch,:)))/std(hjorth_nrem.c(ch,:)))<hthresh);
     artifacts_dualparam_nrem{ch}(th,:)=zeros(length(tempparam),1);
     artifacts_dualparam_nrem{ch}(th,tempparam)=(abs((hjorth_nrem.a(ch,tempparam)-mean(hjorth_nrem.a(ch,tempparam)))/std(hjorth_nrem.a(ch,tempparam)))<hthresh2 & abs((hjorth_nrem.m(ch,tempparam)-mean(hjorth_nrem.m(ch,tempparam)))/std(hjorth_nrem.m(ch,tempparam)))<hthresh2 & abs((hjorth_nrem.c(ch,tempparam)-mean(hjorth_nrem.c(ch,tempparam)))/std(hjorth_nrem.c(ch,tempparam)))<hthresh2);
     tempparam=(abs((hjorth_rem.a(ch,:)-mean(hjorth_rem.a(ch,:)))/std(hjorth_rem.a(ch,:)))<hthresh & abs((hjorth_rem.m(ch,:)-mean(hjorth_rem.m(ch,:)))/std(hjorth_rem.m(ch,:)))<hthresh & abs((hjorth_rem.c(ch,:)-mean(hjorth_rem.c(ch,:)))/std(hjorth_rem.c(ch,:)))<hthresh);
     artifacts_dualparam_rem{ch}(th,:)=zeros(length(tempparam),1);
     artifacts_dualparam_rem{ch}(th,tempparam)=(abs((hjorth_rem.a(ch,tempparam)-mean(hjorth_rem.a(ch,tempparam)))/std(hjorth_rem.a(ch,tempparam)))<hthresh2 & abs((hjorth_rem.m(ch,tempparam)-mean(hjorth_rem.m(ch,tempparam)))/std(hjorth_rem.m(ch,tempparam)))<hthresh2 & abs((hjorth_rem.c(ch,tempparam)-mean(hjorth_rem.c(ch,tempparam)))/std(hjorth_rem.c(ch,tempparam)))<hthresh2);
 end
 end
 
 
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%BUECKELMULLER METHOD%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
    %getting NREM/REM epoch labels
    starttime=1;
    endtime=round(starttime+hyplength*EEG.srate);
    nremwindow_buckelmuller=1;
    remwindow_buckelmuller=1;
    nremtime_buckelmuller=[];
    remtime_buckelmuller=[];
    while endtime<size(EEG.data,2)-hyplength*EEG.srate %trimming off the last hypnogram epoch to be sure
        if sum(hypnogram(ceil(starttime/EEG.srate/hyplength))==[4 5 6])>0 && sum(hypnogram(ceil(endtime/EEG.srate/hyplength))==[4 5 6])>0 
        nremtime_buckelmuller(nremwindow_buckelmuller)=starttime;
        %assigning manual artifacts to epochs
        nremwindow_buckelmuller=nremwindow_buckelmuller+1;
        elseif sum(hypnogram(ceil(starttime/EEG.srate/hyplength))==2)>0 && sum(hypnogram(ceil(endtime/EEG.srate/hyplength))==2)>0 
        remtime_buckelmuller(remwindow_buckelmuller)=starttime;
        remwindow_buckelmuller=remwindow_buckelmuller+1;
        end
        starttime=round(starttime+(hyplength/2)*EEG.srate)+1;
        endtime=starttime+hyplength*EEG.srate;
    end    
    
 clear spectrum_buckelmuller 
 %getting spectra from all epochs
 for ch=1:length(target_electrodes) %looping through channels
     channel=find(contains(electrodes,target_electrodes{ch}));
     channel=channel(1); %multiple matches are at the end of the order
     for nremepoch=1:nremwindow_buckelmuller-1
        signal=EEG.data(channel,nremtime_buckelmuller(nremepoch):round(nremtime_buckelmuller(nremepoch)+hyplength*EEG.srate));
        spectrum_buckelmuller.nrem(ch,nremepoch,:)=periodogram(signal,hamming(length(signal)),bins,EEG.srate);
     end
     for remepoch=1:remwindow_buckelmuller-1
        signal=EEG.data(channel,remtime_buckelmuller(remepoch):round(remtime_buckelmuller(remepoch)+hyplength*EEG.srate));
        spectrum_buckelmuller.rem(ch,remepoch,:)=periodogram(signal,hamming(length(signal)),bins,EEG.srate);
     end 
 end
 
 slowfreq=find(bins>=0.5 & bins<=4.5);
 fastfreq=find(bins>=20 & bins<=40);
 clear slowpower_nrem
 clear fastpower_nrem
 clear slowpower_rem
 clear fastpower_rem
 clear slow_movmean_nrem
 clear fast_movmean_nrem
 clear slow_movmean_rem
 clear fast_movmean_rem
 clear buckelmuller_nrem_slow
 clear buckelmuller_nrem_fast
 clear buckelmuller_rem_slow
 clear buckelmuller_rem_fast
 for ch=1:length(target_electrodes) %looping through channels
 for ep=1:size(spectrum_buckelmuller.nrem,2)
     slowpower_nrem(ch,ep)=sum(spectrum_buckelmuller.nrem(ch,ep,slowfreq));
     fastpower_nrem(ch,ep)=sum(spectrum_buckelmuller.nrem(ch,ep,fastfreq));
 end
 for ep=1:size(spectrum_buckelmuller.rem,2)
     slowpower_rem(ch,ep)=sum(spectrum_buckelmuller.rem(ch,ep,slowfreq));
     fastpower_rem(ch,ep)=sum(spectrum_buckelmuller.rem(ch,ep,fastfreq));
 end
    slow_movmean_nrem(ch,:)=movmean(slowpower_nrem(ch,:),15);
    fast_movmean_nrem(ch,:)=movmean(fastpower_nrem(ch,:),15);
    slow_movmean_rem(ch,:)=movmean(slowpower_rem(ch,:),15);
    fast_movmean_rem(ch,:)=movmean(fastpower_rem(ch,:),15);
    buckelmuller_nrem_slow{ch}=slowpower_nrem(ch,:)./slow_movmean_nrem(ch,:)<2.5;
    buckelmuller_nrem_fast{ch}=fastpower_nrem(ch,:)./fast_movmean_nrem(ch,:)<2;
    buckelmuller_rem_slow{ch}=slowpower_rem(ch,:)./slow_movmean_rem(ch,:)<2.5;
    buckelmuller_rem_fast{ch}=fastpower_rem(ch,:)./fast_movmean_rem(ch,:)<2;
    allnrem=buckelmuller_nrem_slow{ch} & buckelmuller_nrem_fast{ch};
    allrem=buckelmuller_rem_slow{ch} & buckelmuller_rem_fast{ch};
    artifacts_buckelmuller_nrem{ch}=repelem(allnrem,hyplength/windowlength); %we have a problem by design if the hypno window length (the basis for Buckelmuller rejections) is not divisible by the art window length used by other methods
    artifacts_buckelmuller_rem{ch}=repelem(allrem,hyplength/windowlength); 
 end
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%POOLING ALL METHODS%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 clear meanspectrum
 %calculating mean PSDs with the various artifacts
 for ch=1:length(target_electrodes)
     type=1; %pooling manual, single-threshold and dual-threshold spectra
     meanspectrum.nrem{ch}(type,:)=squeeze(nanmean(spectrum.nrem(ch,find(artifacts_manual_nrem),:),2));
     meanspectrum.rem{ch}(type,:)=squeeze(nanmean(spectrum.rem(ch,find(artifacts_manual_rem),:),2));
     type=type+1;
     for th=1:length(params) %single-threshold method
     meanspectrum.nrem{ch}(type,:)=nanmean(spectrum.nrem(ch,find(squeeze(artifacts_singleparam_nrem{ch}(th,:))),:),2);
     meanspectrum.rem{ch}(type,:)=nanmean(spectrum.rem(ch,find(squeeze(artifacts_singleparam_rem{ch}(th,:))),:),2);
     type=type+1;    
     end
     for th=1:length(params2) %dual-threshold method
     meanspectrum.nrem{ch}(type,:)=nanmean(spectrum.nrem(ch,find(squeeze(artifacts_dualparam_nrem{ch}(th,:))),:),2);
     meanspectrum.rem{ch}(type,:)=nanmean(spectrum.rem(ch,find(squeeze(artifacts_dualparam_rem{ch}(th,:))),:),2);
     type=type+1;    
     end
     %Buckelmuller method
     meanspectrum.nrem{ch}(type,:)=nanmean(spectrum.nrem(ch,find(squeeze(artifacts_buckelmuller_nrem{ch})),:),2);
     meanspectrum.rem{ch}(type,:)=nanmean(spectrum.rem(ch,find(squeeze(artifacts_buckelmuller_rem{ch})),:),2);
     type=type+1; 
     %no artifact rejection
     meanspectrum.nrem{ch}(type,:)=nanmean(spectrum.nrem(ch,:,:),2);
     meanspectrum.rem{ch}(type,:)=nanmean(spectrum.rem(ch,:,:),2);
 end
 
 clear spectrumchange
 %calculating the change in average PSD if each artifact is not excluded
 %but included in the analyses
 for ch=1:length(target_electrodes)
     type=1; %only visual rejections are included now, can be changed later
     idx=find(artifacts_manual_nrem==0); %find NREM epochs which are not artifacts
     for ep=1:length(idx)
     artifacts_randommanual_nrem=artifacts_manual_nrem;
     artifacts_randommanual_nrem(idx(ep))=1; %setting a random artifact to non-artifact
     spectrumchange.nrem{ch,type}(ep,:)=squeeze(nanmean(spectrum.nrem(ch,find(artifacts_randommanual_nrem),:),2))';%-meanspectrum.nrem{ch}(type,:);
     end
     idx=find(artifacts_manual_rem==0); %find NREM epochs which are not artifacts
     for ep=1:length(idx)
     artifacts_randommanual_rem=artifacts_manual_rem;
     artifacts_randommanual_rem(idx(ep))=1; %setting a random artifact to non-artifact
     spectrumchange.rem{ch,type}(ep,:)=squeeze(nanmean(spectrum.rem(ch,find(artifacts_randommanual_rem),:),2));%'-meanspectrum.rem{ch}(type,:);
     end
 end
 
 save([edfs{o} '.mat'],'bins','artnames3','meanspectrum','artifacts_manual_nrem','artifacts_manual_rem','artifacts_singleparam_nrem','artifacts_singleparam_rem','artifacts_dualparam_nrem','artifacts_dualparam_rem','artifacts_buckelmuller_nrem', 'artifacts_buckelmuller_rem', 'target_electrodes')
 save([edfs{o} '_spectrumchange.mat'],'spectrumchange')
 save([edfs{o} '_spectrum+hjorth.mat'],'spectrum', 'hjorth_nrem','hjorth_rem')
end


