% this script creates the training data and variables required for each classifier

%% load data
ch = 5;
name = 'RawData_422b';
load([ name '.mat'])

%% offset

for i = 1:6 
    evalc(['thres' num2str(i) '=find(data' num2str(i) '(:,ch+1)==1);']);   % locate times at which each movement begins
    evalc(['data' num2str(i)  '=detrend(data' num2str(i) '(:,1:ch));']);    % remove offset from raw data
end

%% feature extraction
set_t = 500;     % analysis window size (units = samples)
buffer = [1000 1250 1500 1750 2000];

for i = 1:6     % number of classes
    if i == 6
        thres = 4750:250:9500;     % allocate arbitrary windows for resting state
    else
        evalc(['thres = thres' num2str(i) ';']);    % locate windows for each repetition
        thres_buf1 = thres + buffer(1);     % add time buffer to account for human reaction time
        thres_buf2 = thres + buffer(2);
        thres_buf3 = thres + buffer(3);
        thres_buf4 = thres + buffer(4);
        thres_buf5 = thres + buffer(5);
        thres = [thres_buf1; thres_buf2; thres_buf3; thres_buf4; thres_buf5];
    end

    evalc(['emg = data' num2str(i) '(:,1:ch)']);
    
    emg_abs = abs(emg);     %full rectification of EMG signal
    
    % mean absolute value
    MAV = zeros(length(thres),ch);
    for ai = 1:ch   % cycle through channels
        for bi = 1:length(thres)    % cycle through samples
            MAV_temp = emg_abs(thres(bi):thres(bi)+(set_t-1),ai);
            MAV(bi,ai) = mean(MAV_temp);
        end
    end
    
    clear ai bi MAV_temp
    
    % zero crossings
    zmat = zeros(set_t,ch);
    z = zeros(length(thres),ch);
    for ai = 1:ch   % cycle through channels
        for bi = 1:length(thres)    %cycle through samples
            zmat(:,ai) = emg(thres(bi):thres(bi)+(set_t-1),ai);     %temp matrix of
            for ci = 1:set_t-1  %cycle through points
                if zmat(ci,ai)<0 && zmat(ci+1,ai)>0 && abs(zmat(ci,ai)-zmat(ci+1,ai))>=0.01
                    z(bi,ai) = z(bi,ai) + 1;
                elseif zmat(ci,ai)>0 && zmat(ci+1,ai)<0 && abs(zmat(ci,ai)-zmat(ci+1,ai))>=0.01
                    z(bi,ai) = z(bi,ai) + 1;
                end
            end
        end
    end
    
    clear ai bi ci zmat
    % slope sign change
    m_mat = zeros(set_t,ch);
    m_sign = zeros(length(thres),ch);
    for ai = 1:ch   % cycle through channels
        for bi = 1:length(thres)    % cycle through samples
            m_mat(:,ai) = emg(thres(bi):thres(bi)+(set_t-1),ai);
            for ci = 1:set_t-2  % cycle through points
                if m_mat(ci+1,ai)>m_mat(ci+2,ai) && m_mat(ci+1,ai)>m_mat(ci,ai) && abs(m_mat(ci+1,ai)-m_mat(ci,ai))>=0.01
                    m_sign(bi,ai) = m_sign(bi,ai) + 1;
                elseif m_mat(ci+1,ai)<m_mat(ci+2,ai) && m_mat(ci+1,ai)<m_mat(ci,ai) && abs(m_mat(ci+1,ai)-m_mat(ci,ai))>=0.01
                    m_sign(bi,ai) = m_sign(bi,ai) + 1;
                end
            end
        end
    end
    
    clear ai bi ci m_mat
    
    % waveform length
    l_mat = zeros(set_t,ch);
    l = zeros(length(thres),ch);
    for ai = 1:ch   % cycle through channels
        for bi = 1:length(thres)    % cycle through samples
            l_mat(:,ai) = emg(thres(bi):thres(bi)+(set_t-1),ai);
            for ci = 1:set_t-1  % cycle through points
                l(bi,ai) = l(bi,ai) + abs(l_mat(ci+1,ai)-l_mat(ci,ai));
            end
        end
    end
    
    clear ai bi l_mat
    
    % ave matrix
    data_act = [MAV z m_sign l];    % concatenate feature vectors
    evalc (['datamat' num2str(i) '= data_act']);
    clear data_act
end

%% define thresholds
seg = 250;
ave = zeros(10,ch);

for i = 1:6
    evalc(['data = data' num2str(i) '.^2;']);
    if i == 6
        thres = 5000:250:9750;
    else
        evalc(['thres = thres' num2str(i) ';']);
        thres_buf1 = thres + buffer(1);     % add time buffer to account for human reaction time
        thres_buf2 = thres + buffer(2);
        thres_buf3 = thres + buffer(3);
        thres_buf4 = thres + buffer(4);
        thres_buf5 = thres + buffer(5);
        thres = [thres_buf1; thres_buf2; thres_buf3; thres_buf4; thres_buf5];
    end
    for j = 1:10
        evalc(['ave' num2str(i)...
            '(j,:) = mean(data(thres(j):thres(j)+seg-1,1:ch),1)']);     % calculate movement threshold with 250ms moving average
    end
end

thresVAL6 = mean(ave6,1);
thresSTD6 = std(ave6,0,1);

for i = 1:5;
    evalc(['avediff = mean(ave' num2str(i) ',1)-thresVAL6;']);
    thresCH=find(avediff==max(avediff));
    thresVAL=max(avediff)+thresVAL6(thresCH);
    evalc(['thresMAT' num2str(i) '=[thresCH thresVAL];']);
end

thresMAT = [thresMAT1;thresMAT2;thresMAT3;thresMAT4;thresMAT5];

%% LDA 2 category boundaries train data
traindata = [];
tdata = [];

traindata = [datamat1; datamat2; datamat3; datamat4; datamat5; datamat6];   % compile feature vectors
traindatalabel(1:size(datamat1,1)) = 1;     % create training data class labels
traindatalabel(end+1:end+size(datamat2,1)) = 2;
traindatalabel(end+1:end+size(datamat3,1)) = 3;
traindatalabel(end+1:end+size(datamat4,1)) = 4;
traindatalabel(end+1:end+size(datamat5,1)) = 5;
traindatalabel(end+1:end+size(datamat6,1)) = 6;

for i = 1:6     % create two category training data labels
    traindatalabel1 = traindatalabel;
    traindatalabel1(traindatalabel~=i)=2;
    traindatalabel1(traindatalabel==i)=1;    
    [ldaClass,err,P,logp,coeff] = classify(traindata,traindata,traindatalabel1,'linear');
    co(i,:) = coeff(1,2);   % 2 category boundary coefficients
end

clear ldaClass err P logp

All = zeros(6,6);
for i = 1:6
    evalc(['datamat = datamat' num2str(i)]);
    for kk = 1:6
        for j = 1:size(traindata,1)
            All(kk,j) = (traindata(j,:)*co(kk).linear)+co(kk).const;
        end
    end
    %evalc(['aveCo(' num2str(i) ') = mean(f);']);    % compute average dot product value
end

tdata.All = All;
tdata.thresMAT = thresMAT;
tdata.co = co;

num = 50;

m1 = tdata.All(:,1:num);
m2 = tdata.All(:,num+1:2*num);
m3 = tdata.All(:,2*num+1:3*num);
m4 = tdata.All(:,3*num+1:4*num);
m5 = tdata.All(:,4*num+1:5*num);
m6 = tdata.All(:,5*num+1:end);

m1not = [m2(1,:) m3(1,:) m4(1,:) m5(1,:) m6(1,:)];
m2not = [m1(2,:) m3(2,:) m4(2,:) m5(2,:) m6(2,:)];
m3not = [m2(3,:) m1(3,:) m4(3,:) m5(3,:) m6(3,:)];
m4not = [m2(4,:) m3(4,:) m1(4,:) m5(4,:) m6(4,:)];
m5not = [m2(5,:) m3(5,:) m4(5,:) m1(5,:) m6(5,:)];
m6not = [m2(6,:) m3(6,:) m4(6,:) m1(6,:) m5(6,:)];

pd1 = fitdist(m1(1,:)','Normal');
pd1not = fitdist(m1not','Normal');
pd2 = fitdist(m2(2,:)','Normal');
pd2not = fitdist(m2not','Normal');
pd3 = fitdist(m3(3,:)','Normal');
pd3not = fitdist(m3not','Normal');
pd4 = fitdist(m4(4,:)','Normal');
pd4not = fitdist(m4not','Normal');
pd5 = fitdist(m5(5,:)','Normal');
pd5not = fitdist(m5not','Normal');
pd6 = fitdist(m6(6,:)','Normal');
pd6not = fitdist(m6not','Normal');

x = -1000:1000;

for i = 1:6
    evalc(['y' num2str(i) '= pdf(pd' num2str(i) ',x)']);
    evalc(['y' num2str(i) 'not = pdf(pd' num2str(i) 'not,x)']);
    evalc(['y' num2str(i) 'int = cumtrapz(y' num2str(i) ')' ]);
    evalc(['y' num2str(i) 'notint = cumtrapz(flip(y' num2str(i) 'not))']);
    evalc(['y' num2str(i) 'mark = length(find(y' num2str(i) 'int<.01))']);
    evalc(['y' num2str(i) 'notmark = length(find(y' num2str(i) 'notint<.001))']);
    evalc(['ymark(' num2str(i) ') = y' num2str(i) 'mark']);
    evalc(['ynotmark(' num2str(i) ') = y' num2str(i) 'notmark']);
    if x(ymark(i)) < x(end - ynotmark(i))
        ymark(i) = length(x)-ynotmark(i);
    end
end

for i = 1:6
    xmark(i) = x(ymark(i));
end

tdata.xmark = xmark;

save(['..\Full Program\Training Data\Train2Cat' name(find(name=='_'):end) '.mat'],'tdata')

%% LDA boundaries train data
clear co aveCo tdata
[ldaClass,err,P,logp,coeff] = classify(traindata,traindata,traindatalabel,'linear');
co = coeff;     % 6 category boundary coefficients

clear ldaClass err P logp

f = zeros(1,6);
sumf = zeros(1,size(datamat,1));
for i = 1:6
    evalc(['datamat = datamat' num2str(i)]);
    for j = 1:size(datamat,1)
        for k = 1:6
            if i ~= k
                f(k) = (datamat(j,:)*co(i,k).linear)+co(i,k).const;
            end
        end
        sumf(j) = sum(f(f>0));
    end
    evalc(['aveCo(' num2str(i) ') = mean(sumf);']);     % compute average dot product value
end

tdata.aveCo = aveCo;
tdata.thresMAT = thresMAT;
tdata.co = co;

save(['..\Full Program\Training Data\TrainBounds' name(find(name=='_'):end) '.mat'],'tdata')

%% full LDA train data

clear tdata
tdata.traindata = traindata;
tdata.traindatalabel = traindatalabel;
%tdata.stddata = stddata;
tdata.thresMAT = thresMAT;
% 

save(['..\Full Program\Training Data\TrainLDA' name(find(name=='_'):end) '.mat'],'tdata')

clear i j k 


