%DISPLAY spikes over filtered signal
for i=53:53
%    subplot(3,3,i-52);
    plot(1:length(filtered(:,i)),filtered(:,i),'b',1:length(filtered(:,i)),spikes(:,i)/200,'r');
%    xlabel('Time (ms)');
 %   ylabel('Voltage / Spike');
 %   title('Comparison of Filtered Signal and Extracted Spikes');
     end
         %subplot

%CONSTANTS
maxtime= 8;  %2e-7 *20000;        
elec=53;
fac=2;start=1;

%MANUAL IDENTITY of remap02pin02 - trial 53
%iden = [1 1 2 2 2 1 1 1 1 1 1 1 2 2 1 1 2 2 1 2 2 1 1 2 1 1 1 1 2 2 2 1 1 2 2 2 1 1 1 2 1 1];
iden = [1 1 2 2 2 1 1 1 1 1 1 1 2 2 1 1 2 2 1 2 2 1 1 2 1 1 1 1 2 2 2 1 1 1 2 2 1 1 1 2 1 1];

%plot(start:length(filtered(:,elec))/fac,filtered(start:end/fac,elec),'b',start:length(filtered(:,elec))/fac,spikes(start:end/fac,elec)/200,'r');


%PICTURE SPIKES
temp=[];
j=1;
i=11;
for l=11:length(spikes(:,elec))
    if spikes(i,elec)==1
        temp=filtered(i-maxtime:i+maxtime,elec);
        subplot(6,7,j);
       H= plot([-maxtime:maxtime],temp);
       title(iden(j));
        axis off;
        j=j+1;
        i=i+9;
    end
    if j==43
        break;
    end
    i=i+1;
end

%plot([1:maxtime*2+1],temp)


%----------------------------------
%CLASSIFY SPIKES
maxtime= 10;  %2e-7 *20000;        
elec=53;
fac=2;start=1;

%plot(start:length(filtered(:,elec))/fac,filtered(start:end/fac,elec),'b',start:length(filtered(:,elec))/fac,spikes(start:end/fac,elec)/200,'r');
k=0;
i=11;
while 1
    if i >= length(spikes(:,elec))
        break;
    end
    if spikes(i,elec)==1
        k=k+1;
        i=i+9;
    end
    i=i+1;
end
numspikes=k

    %picture them
temp=[];
k=1;
features=zeros(numspikes,5); %maxheight, peak-peakheight
i=11;
while 1
if i >= length(spikes(:,elec))
    break;
end
    if spikes(i,elec)==1
        max=0;max2=0;maxp=0;min=0;tmax=0;tmin=0;
        temp=filtered(i-maxtime:i+maxtime,elec);
        %maxheight
        for j=1:maxtime*2+1
            if (temp(j))^2 > max
                max2=max;
                max= temp(j);
            elseif (temp(j))^2 > max2
                max2 = temp(j);
            end
        end
        features(k,1)=max;
        features(k,2)=max-max2;
        %peak-peak height
        for j=1:maxtime*2+1
            if temp(j) > maxp
                maxp = temp(j);
                tmax=j;
            elseif temp(j) < min
                min = temp(j);
                tmin=j;
            end
        end
        features(k,3)= maxp+min;
        features(k,4)= tmax-tmin;
        counta=1;countb=1;swidth=0;
        for j=2:maxtime*2+1
            if temp(j)*temp(j-1) < 0
                counta= countb;
                countb = j;
                if countb-counta > swidth
                    swidth = countb-counta;
                end
            end
        end
        features(k,5)=swidth;
        k=k+1;
        i=i+9;
    end
    i=i+1;
end

%------------------
%PLOT SPIKE BY CLASSIFICATION
%--------------------
figure(2);
ssize=(features(:,5)-2)*20;
scolor=[1-(features(:,1))*180, zeros(length(features(:,1)),2)];
scolor(40,1,1) =1;
H=scatter3(features(:,4),features(:,3),features(:,2),ssize, scolor,'filled');
axis vis3d;
pause
for i=1:7
   rotate(H,[2,0,4],i);
    pause(1);
end

%-----------------
%PLOT SPIKE BY CLASSIFICATION AND BY INSPECTION TYPE
%----------------
figure(3);
features_1=[];features_2=[];
for i=1:length(features(:,1))
    if iden(i)==1
        features_1=[features_1;features(i,:)];
    else
        features_2=[features_2;features(i,:)];
    end
end

figure;
ssize_1=(features_1(:,5)-2)*20; ssize_2=(features_2(:,5)-2)*20;
scolor_1=[1-(features_1(:,1))*180, zeros(length(features_1(:,1)),2)]; scolor_2=[1-(features_2(:,1))*180, zeros(length(features_2(:,1)),2)];

G=scatter3(features_1(:,4),features_1(:,3),features_1(:,2),ssize_1, scolor_1,'x');
hold
H=scatter3(features_2(:,4),features_2(:,3),features_2(:,2),ssize_2, scolor_2,'o');
axis vis3d;
for i=1:100
    pause
   rotate(H,[1,0,0],i);
   rotate(G,[1,0,0],i);
end

scatter(features_1(:,4),features_1(:,3),'x');
hold
scatter(features_2(:,4),features_2(:,3),'o');

%EXTRAS
k=1;
for i=1:4
    for j=i+1:4
        subplot(2,3,k);
        scatter(features(:,i),features(:,j));
        k=k+1;
    end
end
%1,2 1,3 1,4 2,3 2,4 3,4

k=1;
for i=1:2
    for j=i+1:3
        for n=j+1:4
            subplot(2,2,k);
            scatter3(features(:,i),features(:,j), features(:,n));
            k=k+1;
        end
    end
end
%1,2,3 1,2,4 1,3,4 2,3,4

subplot(1,3,1);
scatter(features(:,1),features(:,3),'filled');
subplot(1,3,2);
scatter(features(:,1),features(:,4),'filled');
subplot(1,3,3);
scatter(features(:,3),features(:,4),'filled');

figure(2);
H=scatter3(features(:,1),features(:,3),features(:,5),'filled');
axis vis3d;
for i=1:18
   rotate(H,[0,0,1],i);
    pause(.5);
end
for i=1:180
   rotate(H,[0,0,1],i);
    pause(.5);
end  

k=1;
corrc=zeros(4);
for i=1:4
    for j=1:4
        temp=corrcoef(features(:,i),features(:,j));
        corrc(i,j)=temp(1,2);
    end
end