function outMat = prepareEEGfwdModel(outData, tRange, wSize, stepSize, dirPad, modelParam, bkOpt, bkWin)
%
% A routine for preparing data to apply forward encoding/decoding model analysis
%
% outData: output of reOrganizeSortedEEG.m function
% tRange: minimum and maximum time from stimulus motion onset ex) [-200 400]
% wSize: time window size for averaging (in milliseconds)
% stepSize: step size of sliding time window 
% dirAngles: Directions for foward model in degree. ex) [0:30:350]
% modelParam: Parameters for encoding/decoding forward model (tuning function)
% 
% prepareEEGfwdModel.m
% 03/01/2017 Joonyeol Lee
% joonyeol@skku.edu
%

behavior = [];
behavSpd = [];
behavDir = [];
accelSpd = [];
accelDir = [];
fitQual = [];


dirIdx = [];
dirValues = [];
trIdx = [];
eegData = [];
stimDirs = [];
chNum = 64;

if exist('bkWin') & ~isempty(bkWin)
	bkgWin = bkWin;
else
	bkgWin = [-199:-50];
end

if isfield(outData, 'info')
	preStim = outData.info.preStim;
else
	preStim = 400;
end

tVec = [tRange(1)+preStim:stepSize:tRange(2)+preStim-wSize];

timeTick = [tRange(1)+round(wSize/2):stepSize:tRange(2)-round(wSize/2)];

for i = 1:length(outData.data)
	stimDirs = [stimDirs;outData.data(i).stim.direction];
end

tDir = find(stimDirs < 0);
stimDirs(tDir) = stimDirs(tDir)+360;

sDir1 = [];
sDir2 = [];
sDir1 = find(stimDirs <= 120 & stimDirs >= 0);
sDir2 = find(stimDirs >= 240);

if ~isempty(sDir1) & ~isempty(sDir2)
	stimDirs(sDir2) = stimDirs(sDir2) - 360;
end

for i = 1:length(outData.data)
	cData = outData.data(i);
	if isfield(cData, 'behavior')
		sze = length(cData.behavior.time);
		behavior = cat(1, behavior, [cData.behavior.spd cData.behavior.dir+cData.behavior.rotDeg cData.behavior.time]);
		tmpH = cData.behavior.hVel;
		tmpV = cData.behavior.vVel;

		tmpHV = ft_preproc_lowpassfilter(tmpH, 1000, 15, 2, 'but', 'twopass');
		tmpVV = ft_preproc_lowpassfilter(tmpV, 1000, 15, 2, 'but', 'twopass');

		tmpH_accel = diff(tmpHV, 1, 2).*1000;
		tmpV_accel = diff(tmpVV, 1, 2).*1000;
		
		smH = nan(size(tmpH,1), length(tVec));
		smV = nan(size(tmpV,1), length(tVec));
		smH_accel = nan(size(tmpH,1), length(tVec));
		smV_accel = nan(size(tmpH,1), length(tVec));
		for t = 1:length(tVec)
			smH(:, t) = mean(tmpH(:, tVec(t):tVec(t)+wSize-1),2);
			smV(:, t) = mean(tmpV(:, tVec(t):tVec(t)+wSize-1),2);
			
			smH_accel(:, t) = mean(tmpH_accel(:, tVec(t):tVec(t)+wSize-1),2);
			smV_accel(:, t) = mean(tmpV_accel(:, tVec(t):tVec(t)+wSize-1),2);
		end
		behavSpd = cat(1, behavSpd, sqrt(smH.^2+smV.^2));
		behavDir = cat(1, behavDir, atan2d(smV, smH));
		accelSpd = cat(1, accelSpd, sqrt(smH_accel.^2+smV_accel.^2));
		accelDir = cat(1, accelDir, atan2d(smV_accel, smH_accel));
		fitQual = cat(1, fitQual, cData.behavior.thresh);
	else
		sze = size(cData.eeg, 3);
	end
	dirValues = cat(1, dirValues, stimDirs(i)*ones(sze, 1));
	trIdx = cat(1, trIdx, transpose([1:sze]));
	eegData = cat(3, eegData, cData.eeg);
end

NoT = length(trIdx);
NoC = size(eegData,1);

tmpInfo = [trIdx dirValues];
[sortedInfo sortIdx] = sortrows(tmpInfo, 2);

dirCond = unique(sortedInfo(:,2));
tmpDir = sort(stimDirs);
switch dirPad
case 0
	dirAngles = tmpDir;
case 1
	dirAngles = [tmpDir(1)-30;tmpDir;tmpDir(end)+30];
case 2
	dirAngles = [tmpDir(1)-2*30;tmpDir(1)-1*30;tmpDir;tmpDir(end)+1*30;tmpDir(end)+2*30];
case 3
	dirAngles = [tmpDir(1)-3*30;tmpDir(1)-2*30;tmpDir(1)-1*30;tmpDir;tmpDir(end)+1*30;tmpDir(end)+2*30;;tmpDir(end)+3*30];
end

numAngles = length(dirAngles);

tunModel = modelParam{1};

p = modelParam{2};

C = nan(numAngles, length(dirCond));

for i = 1:length(dirCond)
	currSze = length(find(sortedInfo(:,2) == dirCond(i)));
	dirIdx = cat(1, dirIdx, i*ones(currSze, 1));
	[C(:,i)]=fwdModelTuningCurve(dirAngles, dirCond(i), tunModel, p);
end

stimInfo = [sortIdx sortedInfo dirIdx];
if isfield(cData, 'behavior')
	behavior = behavior(sortIdx,:);
	behavSpd = behavSpd(sortIdx,:);
	behavDir = behavDir(sortIdx,:);
	accelSpd = accelSpd(sortIdx,:);
	accelDir = accelDir(sortIdx,:);
	fitQual = fitQual(sortIdx);
end
eegData = eegData(:,:,sortIdx);

clear dirIdx dirValues trIdx sortIdx;

Cmat = nan(numAngles, NoT);

for i = 1:NoT
	Cmat(:,i) = C(:, stimInfo(i, 4));
end

zEEG = eegData;
chV = eegData(:);
chMean = mean(chV(:));
chStd = std(chV(:));
zEEG = (eegData - chMean)./chStd;
tmpEEG = zEEG;

bWin = bkgWin+preStim;

if exist('bkOpt') & bkOpt == 1 
	ntmpEEG = tmpEEG;
	for i = 1:NoC
		for j = 1:NoT
			ntmpEEG(i,:,j) = tmpEEG(i,:,j) - mean(squeeze(tmpEEG(i,bWin,j)));
		end
	end
	tmpEEG = ntmpEEG;
end
	
binnedEEG = nan(NoC, length(tVec), NoT);

for i = 1:NoC
	if hopt == 2
		for k = 1:length(tVec)
			for j = 1:NoT
				binnedEEG(i,k,j) = circ_mean(squeeze(tmpEEG(i, tVec(k):tVec(k)+wSize-1, j)), [], 2);
			end
		end
	else
		for k = 1:length(tVec)
			for j = 1:NoT
				binnedEEG(i,k,j) = mean(squeeze(tmpEEG(i, tVec(k):tVec(k)+wSize-1, j)));
			end
		end
	end
end

zEEG = binnedEEG;

zEEG = zEEG.*1e8;

rng('shuffle', 'twister'); 

randIdx = transpose(randperm(NoT));

outMat.binEEG = zEEG(:,:,randIdx);
outMat.eegTimeTick = timeTick;
outMat.stim = stimInfo(randIdx,:);
if ~isempty(behavior)
	outMat.behavior = behavior(randIdx,:);
	outMat.behavSpd = behavSpd(randIdx,:);
	outMat.behavDir = behavDir(randIdx,:);
	outMat.accelSpd = accelSpd(randIdx,:);
	outMat.accelDir = accelDir(randIdx,:);
	outMat.fitQual = fitQual(randIdx);
end
outMat.fwdModel = Cmat(:, randIdx);
outMat.rIdx = randIdx;
outMat.angles = dirAngles;