%% EEG TIME FREQUENCY ANALYZE SCRIPT
% computing EEG time-frequency PHASE-LOCKED/NON-PHASE-LOCKED/TOTAL-POWER

%INPUT:
%   - EEGlab study set.
%OUTPUT:
%   - n - dimensional matrix; size( channels x frequencies x time-points)

%Bartlomiej Panek,  pnk.bartlomiej@gmail.com
%%
clear
close
eeglab
%%
    
    EEG_temp = EEG;  
    EEG_temp.data = double(EEG_temp.data); 
    pnts    = EEG_temp.pnts;
    trials  = EEG_temp.trials;
    nbchans = EEG_temp.nbchan;
             
%baseline    
    timevec  =  -1500:2:800; % down-sampling for many trials
    timeidx = dsearchn(EEG_temp.times',timevec');

%%  Wavelet 
    mtime = -2:1/EEG_temp.srate:2;
    htime = (length(mtime)-1)/2;
    num_frex    = 30;
    min_freq    = 1;
    max_freq    = 30;
    frex        = logspace(log10(min_freq),log10(max_freq),num_frex); 
    rcycles     = [3 8];
    ncycles     = logspace(log10(rcycles(1)),log10(rcycles(end)),num_frex);
    
% FFT parameters
    kernel = length(mtime);
    ndata  = pnts*trials;
    nconv  = kernel+ndata-1;
    
 % WAVELET FFT    
    cmwFFT   = zeros(length(frex),nconv);
    chanidx  = zeros(nbchans);
    mw_phase = zeros(nbchans,length(frex),length(timevec),trials);
    
for fi=1:num_frex
    s = ncycles(fi)/(2*pi*frex(fi));
    cmw = exp(2*1i*pi*frex(fi).*mtime) .* exp(-mtime.^2./(2*s^2)); %complex morlet wavelet creation

    temp = fft(cmw,nconv);
    cmwFFT(fi,:) = temp./max(temp); %uV normalization units  
end

   

%% MAIN LOOP

for chani=1:nbchans
    dataFFT = fft(reshape(EEG_temp.data(chani,:,:),1,[]),nconv);
    
    for fi=1:length(frex) 
        as = ifft(cmwFFT(fi,:).*dataFFT,nconv);
        as = as(htime+1:end-htime);
        as = reshape(as,pnts,trials);

        mw_phase(chani,fi,:,:) = angle(as(timeidx,:));
        
    end
 end 


%% SAVE
    % Analyze settings
    MW_settings.mwtime   =  {save_time};
    MW_settings.LAP      =  {'yes'};
    MW_settings.wavelet  =  [{rcycles}, {frex}];
    MW_settings.ntrials  =  {trials};
    %T-F matrices
    MW_results.mw_phase    = single(mw_phase);   
    
    nameTag_raw = sprintf(strcat(EEG_temp.setname,'_PHASE'),'%s');
    save (nameTag_raw,'frex','chanlocs','num_frex','MW_results','MW_settings','timevec')

%%
phase_vals(1,:,:,:) = squeeze(MW_results.mw_phase(25,:,:,:)); %PO7
phase_vals(2,:,:,:) = squeeze(MW_results.mw_phase(47,:,:,:)); %FCz