ChronicGui/alignSpkWavXcovMulti.m at master · jamesjun/ChronicGui · 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
function tr = alignSpkWavXcovMulti(tr, varargin)
% using neo

P = funcDefStr(funcInStr(varargin{:}), 'viChan', [], 'maxLag', 5, 'nRepeat', 10, 'nRepeatMin', 3, 'nTrim', 4);

nSpks = size(tr,3);
nSamples = size(tr,1);
nChans = size(tr,2);
viRange0 = 1:nSamples;
if isempty(P.viChan)
    P.viChan = 1:nChans;
end

switch numel(size(tr))
    case 3
        nShiftPrev = inf;
        for iRepeat=1:P.nRepeat
            nShift = 0;
            for iChan=1:nChans
                mr = tr2mr(tr,iChan);
                vr0 = median(mr,2); %template
                for iSpk=1:nSpks
                    vr = mr(:,iSpk);
                    [vrXcv, vnLags] = xcov(vr0, vr, P.maxLag, 'coef');
                    [~, iLag] = max(vrXcv);
                    if vnLags(iLag) ~= 0
                       nShift = nShift+1;
                    end
                    viRange = mod(viRange0 - vnLags(iLag) - 1, nSamples)+1;
                    tr(:,:,iSpk) = tr(viRange,:,iSpk); %waveform shifted, move other channels
                end
            end
            fprintf('nShift=%d\n', nShift);
            if nShiftPrev <= nShift && iRepeat > P.nRepeatMin
                break;
            end
            nShiftPrev = nShift;
        end
        if ~isempty(P.nTrim)
            tr = tr(P.nTrim+1:end-P.nTrim,:,:);
        end
    case 2
        nSpks = size(tr,2);
        for iRepeat=1:P.nRepeat
            vr0 = mean(tr,2); %template
            nShift = 0;
            for iSpk=1:nSpks
                vr = tr(:,iSpk);
                [vrXcv, vnLags] = xcov(vr0, vr, P.maxLag, 'coef');
                [~, iLag] = max(vrXcv);
                if vnLags(iLag) ~= 0
                   nShift = nShift+1;
                end
                viRange = mod(viRange0 - vnLags(iLag) - 1, nSamples)+1;
                tr(:,iSpk) = vr(viRange); %waveform shifted
            end
            fprintf('nShift=%d\n', nShift);
        end
        if ~isempty(P.nTrim)
            tr = tr(P.nTrim+1:end-P.nTrim,:);
        end
end %switch

end %function


function test()
mrA = repmat(sin(0:.1:2*pi)', [1,1,10]);
mrA(:,:,5) = [mrA(5:end,:,1); mrA(1:4,:,1)];
figure; imagesc(mrA(:,:,5));
end