%spectragram correlation
%crosscorrelate the two spectrogram images
%and look for the peak value

clear all

global MaxSpecTable MaxSpeedTable



%scrape
filename{1} = 'hdosata31to700' ;
filename{2} = 'hdosata21to490' ;
filename{3} = 'hdosata11to800';
filename{4} = 'Hdosata4750to1550' ;
filename{5} = 'Hdosata5a490to1470' ;
filename{6} = 'Hdosata6b280to1050' ;
%thump
filename{7} = 'hdosata3820to1150' ;
filename{8} = 'hdosata2550to880' ;
filename{9} = 'hdosata11290to1540' ;
filename{10} = 'Hdosata3b1520to1800' ;
filename{11} = 'Hdosata41560to1820';
filename{12} = 'Hdosata5a1480to1750'; %noise
filename{13} = 'Hdosata5b1040to1230';
filename{14} = 'Hdosata61470to1640';
filename{15} = 'Hdosata6c770to1130';
filename{16} = 'Hdosata71510to1700';
%buzz
filename{17} = 'hdosata31400to1900' ;
filename{18} = 'hdosata21140to1940';
filename{19} = 'hdosata1b1260to1800' ;%different
filename{20} = 'Hdosata5b1490to2040';
filename{21} = 'Hdosata6c1310to1986';
nfile = length(filename)

lowcut = 4 ;

for i=1:nfile
    i
    for j=i:nfile
        load(['c:\documents and settings\bruce land\my documents\matlab\damian\waveforms6dec03\',...
                filename{i},...
                'gram.mat'],'spdhist');
        
        load(['c:\documents and settings\bruce land\my documents\matlab\damian\waveforms6dec03\',...
                filename{i},...
                'wave.mat'],'s');
        spdhist1 = log10(spdhist(:,lowcut:end));
        spdhist1(spdhist1<0) = 0 ;
        nT = size(spdhist,1);
        avgspdhist = repmat(sum(spdhist1,1)/nT,nT,1);
        spdhist1 = spdhist1 - avgspdhist ;
        diffhist1 = spdhist1-mean(mean(spdhist1));
        s1 = s-mean(s);
        
        load(['c:\documents and settings\bruce land\my documents\matlab\damian\waveforms6dec03\',...
                filename{j},...
                'gram.mat'],'spdhist');
        
        load(['c:\documents and settings\bruce land\my documents\matlab\damian\waveforms6dec03\',...
                filename{j},...
                'wave.mat'],'s');
        spdhist2 = log10(spdhist(:,lowcut:end));
        spdhist2(spdhist2<0) = 0 ;
        nT = size(spdhist,1);
        avgspdhist = repmat(sum(spdhist2,1)/nT,nT,1);
        spdhist2 = spdhist2 - avgspdhist ;
        diffhist2 = spdhist2-mean(mean(spdhist2));
        s2 = s-mean(s);
        
        %cross-corr the speed-o-grams
        len1=length(diffhist1);
        len2=length(diffhist2);
        h = size(diffhist1,2);
        if len2<len1
            diffhist2 = [diffhist2; mean(mean(diffhist2))*ones(len1-len2,h)] ;
        elseif len1<len2
            diffhist1 = [diffhist1; mean(mean(diffhist1))*ones(len2-len1,h)] ; ;
        end
        
        %the actual circular cross-corr
        hcorr=[];
        normfactor = sqrt(sum(sum(diffhist1.*diffhist1)).*sum(sum(diffhist2.*diffhist2)));
        for ii=1:length(diffhist1)
            hcorr(ii) = sum(sum(diffhist1.*diffhist2))/normfactor ;
            diffhist1 = [diffhist1(end,:); diffhist1(1:(end-1),:)];
        end
        
        %make the speed sequences the same length
        len1=length(s1);
        len2=length(s2);
        if len2<len1
            s2 = [s2; mean(s2)*ones(len1-len2,1)];
        elseif len1<len2
            s1 = [s1; mean(s1)*ones(len2-len1,1)];
        end
        
        %compute the crosscorr, normalized by the autocorr of each trace
        %         autos1 = xcov(s1,s1,'unbiased');
        %         autos2 = xcov(s2,s2,'unbiased');
        %         xcv = xcov(s1,s2,'unbiased')/sqrt(autos1(length(s1))*autos2(length(s2)));
        normfactor = sqrt(sum(s1.*s1).*sum(s2.*s2));
        for ii=1:length(s1)
            xcv(ii) = sum(s1.*s2)/normfactor ;
            s1 = [s1(end); s1(1:(end-1))];
        end
%         disp([filename{i},'//',filename{j},...
%                 '  maxspeccorr=',num2str(max(hcorr),'%5.2f'),...
%                 '  maxspeedcorr=',num2str(max(xcv),'%5.2f')])
        
        MaxSpecTable(i,j)=max(hcorr);
        MaxSpeedTable(i,j)=max(xcv);
        MaxSpecTable(j,i)=MaxSpecTable(i,j);
        MaxSpeedTable(j,i)=MaxSpeedTable(i,j);
    end
    figure(3)
    clf
    colormap(gray(30))
    imagesc(spdhist1')
    axis image
    set(gca,'ydir','normal')
    drawnow
    %pause
end

figure(1)
clf
colormap(gray(50))
imagesc((1-MaxSpecTable).^2)
title('Max Correlation Spectral Table')

figure(2)
clf
colormap(gray(50))
imagesc((1-MaxSpeedTable).^2)
title('Max Correlation Speed Table')


%pick a distance weighting exponent
k = -2 ;

%define response group membership as a CELL array
%( note curly brackets{}! )
%so that different groups can have different numbers of traces
ngrp = 3;
grp{1} = 1:6  ; 
grp{2} = 7:17 ;
grp{3} = 17:21 ;


%make array of avg dist of each train to each group 
%as an (nNeuron,ngrp) array
%excluding the current train from the average
for i=1:nfile
    for j = 1:ngrp
        traceIndex = grp{j}(find(grp{j}~=i)) ;  %excludes the current trace
        avgDist(i,j) = ( mean((1-MaxSpeedTable(i,traceIndex)).^k) ).^(1/k) ;
    end
end

%find minimum in each row (for each train)
%make N matrix (see Victor)
%with a row per stimulus group and a column per response group
N=zeros(ngrp);
for i=1:nfile
    [minv,minResponseGrp] = min(avgDist(i,:));
    %find out which stim group the trace is in
    for j=1:ngrp
        if (sum(grp{j}==i)==1)
            StimGrp = j;
        end
    end
    %form the confusion matrix
    N(StimGrp,minResponseGrp) = N(StimGrp,minResponseGrp) + 1;
end


%get the number of bits available--
%max info is log2(ngrp)
Ntotal = sum(sum(N)); 

for i=1:ngrp
    for j=1:ngrp
        h(i,j) = N(i,j) * ...
            ( log2(N(i,j)+.0001) - log2(sum(N(:,j))+.0001) - ...
            log2(sum(N(i,:))+.0001) + log2(Ntotal) );     
    end
end

Hspeed = sum(sum(h))/Ntotal  
    

%make array of avg dist of each train to each group 
%as an (nNeuron,ngrp) array
%excluding the current train from the average
for i=1:nfile
    for j = 1:ngrp
        traceIndex = grp{j}(find(grp{j}~=i)) ;  %excludes the current trace
        avgDist(i,j) = ( mean((1-MaxSpecTable(i,traceIndex)).^k) ).^(1/k) ;
    end
end

%find minimum in each row (for each train)
%make N matrix (see Victor)
%with a row per stimulus group and a column per response group
N=zeros(ngrp);
for i=1:nfile
    [minv,minResponseGrp] = min(avgDist(i,:));
    %find out which stim group the trace is in
    for j=1:ngrp
        if (sum(grp{j}==i)==1)
            StimGrp = j;
        end
    end
    %form the confusion matrix
    N(StimGrp,minResponseGrp) = N(StimGrp,minResponseGrp) + 1;
end


%get the number of bits available--
%max info is log2(ngrp)
Ntotal = sum(sum(N)); 

for i=1:ngrp
    for j=1:ngrp
        h(i,j) = N(i,j) * ...
            ( log2(N(i,j)+.0001) - log2(sum(N(:,j))+.0001) - ...
            log2(sum(N(i,:))+.0001) + log2(Ntotal) );     
    end
end

Hspec = sum(sum(h))/Ntotal