% extraction des data pour une waveform analysis

clear all; clc;
task_name={'DT5','FR5','FS5'};
task_Dataset={'R_DT5','R_FR5','R_FS5'};
figure;
for task=1:length(task_Dataset)
    clear R selection X idx sizeplot color Celltype SesCelltype
    load(task_Dataset{task})

    WFlength=800; %total duration of the waevform 800microseconds
    bins=25;%bins=25microsecond
    time=(bins:bins:WFlength);
    runAnalysis=1;
    
    for i=1:size(R.WF,1)
        if sum(R.WF(i,:))==0
            R.WFanalysis(i,2)=0;
            R.WFanalysis(i,3)=0;
        else
            [maxValue, MaxIndex] = max(R.WF(i,:));
            [minValue, MinIndex] = min(R.WF(i,:));
            
            PeakValleyW=abs(time(MaxIndex)-time(MinIndex));
            PeakValleyM=abs(maxValue-minValue);
            R.WFanalysis(i,2)=PeakValleyM;
            
            idx = find(diff(R.WF(i,:) >= minValue/2));
            x2 = time(idx) + (minValue/2 - R.WF(i,idx)) .* (time(idx+1) - time(idx)) ./ (R.WF(i,idx+1) - R.WF(i,idx));
            if length(x2)==1
                R.WFanalysis(i,3)=NaN;
            else
                ValleyFWHM=x2(2)-x2(1);
                R.WFanalysis(i,3)=ValleyFWHM;
            end
        end
    end
    
    selection=R.WFanalysis(:,2)~=0;
    X=R.WFanalysis(selection,:);
     
    idx(1:sum(selection),1)=1;
    idx(X(:,1)>20,1)=2; % threshold 20Hz to separate MSN from FSI
    idx(idx(:,1)==2 & X(:,3)>150,1)=0;
    idx(idx(:,1)==1 & X(:,3)>450,1)=3;
    
    idx(X(:,1)>5 & idx(:,1)==3,1)=0;
    idx(X(:,1)<20 & X(:,1)>12.5 & idx(:,1)==1,1)=0; % intermediate frequencies = not classified
    idx(X(:,3)<450 & X(:,3)>400 & idx(:,1)==1,1)=0; % intermediate frequencies = not classified
    
    nbExcludedN=sum(idx(:,1)==0);
    
    
    for Y=1:size(idx,1)
        if idx(Y,1)==1
            color(Y,:)=[255/255 0/255 0/255];%MSN
            sizeplot(Y,1)=50;
        elseif idx(Y,1)==2
            color(Y,:)=[0 255/255 0/255];%FSI
            sizeplot(Y,1)=50;
        elseif idx(Y,1)==3
            color(Y,:)=[0/255 0/255 255/255]; %TAN
            sizeplot(Y,1)=50;
        elseif idx(Y,1)==0
            color(Y,:)=[0/255 0/255 0/255];%unclassified
            sizeplot(Y,1)=5;
        end
    end
    
    subplot(3,4,[task 4+task 8+task])
    scatter(X(:,1),X(:,3),sizeplot,color)
    hold on
    title(task_name{task})
    xlabel('Firing rate') % x-axis label
    ylabel('half-width')% y-axis label
    axis([0 50,0 600])
    hold off
    
    Celltype(1:length(selection),1)=R.SESSION;
    Celltype(1:length(selection),2)=0;
    Celltype(selection,2)=idx(:,1);
    
    
    for i=1:max(Celltype(:,1))
        firstUnitID=find(Celltype(:,1)==i,1,'first');
        LastUnitID=find(Celltype(:,1)==i,1,'last');
        SesCelltype(i).Celltype=Celltype(firstUnitID:LastUnitID,2);
    end
    
    save([task_Dataset{task},'.mat'],'R','Celltype','SesCelltype')    
end

%% examples of waveforms

subplot(3,4,4)
TAN_id=find(Celltype(:,2)==3);
plot(time,R.WF(TAN_id(2),:))
axis([0 800 -0.15 0.15 ])
 
subplot(3,4,8)
MSN_id=find(Celltype(:,2)==1);
plot(time,R.WF(MSN_id(2),:))
axis([0 800 -0.4 0.4])
 
subplot(3,4,12)
FSI_id=find(Celltype(:,2)==2);
plot(time,R.WF(FSI_id(2),:))
axis([0 800 -0.2 0.1])