quickSort/doubleSpike.m at master · BellWilson2016/quickSort · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
% Uses a pairwise template matching to try to resolve complex spike
% waveforms. Will match to single spike if that's best-fit.
function data = doubleSpike(data, timePoint)

    spikeCategories = ['A','B','C','D','E','F','G'];

    % Assemble a list of possible pairs
    nClusters = length(data.spikeAvg);
    clusterList = 1:nClusters;
    pairList = [0,0];
    for baseSpike = 1:nClusters
        pairList(end+1,:) = [baseSpike,0];
        for addSpike = (baseSpike + 1):nClusters
            pairList(end+1,:) = [baseSpike,addSpike];
        end
    end

    % Find the snippet to work on.
    centerSample = round(data.sampleRate*timePoint);
    spikeWidth = length(data.spikeAvg{1});
    stSamp = centerSample - round(1.5*spikeWidth);
    enSamp = centerSample + round(1.5*spikeWidth);
    snippet = data.dVdT(stSamp:enSamp);

    template = [];
    spikeCodes = [];
    for pairN = 1:size(pairList,1)
        ix1 = pairList(pairN,1);
        ix2 = pairList(pairN,2);
        if (ix1 > 0)
            spike1 = data.spikeAvg{ix1};
        else
            spike1 = zeros(length(data.spikeAvg{1}),1);
        end
        if (ix2 > 0)
            spike2 = data.spikeAvg{ix2};
        else
            spike2 = zeros(length(data.spikeAvg{1}),1);
        end

        spikeHalfWidth = floor(spikeWidth/2);
        halfNlags = spikeHalfWidth;
        Nlags = 2*halfNlags + 1;
        templateLength = Nlags + spikeWidth + 1;
        centerTemplate = floor(templateLength/2) + 1;

        if (ix2 > 0)
            aSamples = zeros(Nlags,1) + centerTemplate;
            bSamples = [-halfNlags:halfNlags]' + centerTemplate;

            aRaster = zeros(templateLength, Nlags);
            bRaster = zeros(templateLength, Nlags);
            for n=1:length(aSamples)
                aRaster(aSamples(n),n) = 1;
                bRaster(bSamples(n),n) = 1;
            end

            % Keep track of what spikes are when
            someCodes = [ones(Nlags,1)*ix1,ones(Nlags,1)*ix2, aSamples, bSamples];
            spikeCodes = cat(1,spikeCodes,someCodes);
        else
            aSamples = centerTemplate;
            bSamples = centerTemplate;
            aRaster = zeros(templateLength, 1);
            bRaster = zeros(templateLength, 1);
            aRaster(aSamples,1) = 1;
            bRaster(aSamples,1) = 1;

            someCodes = [ones(1,1)*ix1,ones(1,1)*ix2, aSamples, bSamples];
            spikeCodes = cat(1,spikeCodes,someCodes);
        end

         C1 = convn(aRaster, spike1, 'same');
         C2 = convn(bRaster, spike2, 'same');
         template = cat(2,template,C1 + C2);
    end

%     figure;
%     subplot(2,1,1);
%     image(aRaster' + bRaster','CDataMapping','scaled');
%     subplot(2,1,2);
%     image(template','CDataMapping','scaled')


    for startSample = 1:(length(snippet) - size(template,1))
        subSnip = snippet(startSample:(startSample+size(template,1)-1));

        snipMatrix = subSnip*ones(1,size(template,2));
        diffMatrix = snipMatrix - template;
        matchFactor(startSample,:) = std(diffMatrix,1)./std(snipMatrix,1);
    end

    [col, row] = find(matchFactor == min(matchFactor(:)));
    col = col(1); row = row(1); % Only take the first match.

%     close all;
% figure;
% subplot(2,1,1);
% image(-matchFactor','CDataMapping','scaled'); hold on;
% scatter(col,row,'bo');
% subplot(2,1,2);
% plot(stSamp:enSamp,data.dVdT(stSamp:enSamp)); hold on;
% plot(stSamp + col - 1 + (1:templateLength),template(:,row),'m');
% scatter(spikeCodes(row,3) + col + stSamp,data.dVdT(spikeCodes(row,3) + col + stSamp),'ro');
% scatter(spikeCodes(row,4) + col + stSamp,data.dVdT(spikeCodes(row,4) + col + stSamp),'go');hold off;
% axis tight;
% pause;

% Put the new spikes back in the structure
if spikeCodes(row,1) > 0
    data.spikeClusters(end+1) = spikeCodes(row,1);
    data.spikeSamples(end+1)  = spikeCodes(row,3) + col + stSamp - 2;
    data.spikeEmbedding(end+1,:) = [NaN,NaN];
    if spikeCodes(row,2) > 0
        data.spikeClusters(end+1) = spikeCodes(row,2);
        data.spikeSamples(end+1)  = spikeCodes(row,4) + col + stSamp - 2;
        data.spikeEmbedding(end+1,:) = [NaN,NaN];
        disp(['pair-matched ',spikeCategories(spikeCodes(row,1)),' / ',spikeCategories(spikeCodes(row,2))]);
    else
        disp(['pair-matched ',spikeCategories(spikeCodes(row,1)),' / null']);
    end
else
    disp(['pair-matched null / null']);
end


% addSpikes(1).clustN = spikeCodes(row,1);
% addSpikes(2).clustN = spikeCodes(row,2);
% addSpikes(1).sample = spikeCodes(row,3) + col + centerSample;
% addSpikes(2).sample = spikeCodes(row,4) + col + centerSample;