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

clear all

global MaxSpecTable MaxSpeedTable


filename{1} = 'hpugcat_1' ;
filename{2} = 'hpugcat_2' ;
filename{3} = 'hpugcat_3' ;
filename{4} = 'hpugcat_4' ;
filename{5} = 'hpugcat_5' ;

filename{6} = 'hpuggal_1b';
filename{7} = 'hpuggal_2a';
filename{8} = 'hpuggal_2b';
filename{9} = 'hpuggal_3';
filename{10} = 'hpuggal_4a';
filename{11} = 'hpuggal_4b';
filename{12} = 'hpuggal_6';
filename{13} = 'hpuggal_7';

filename{14} = 'hpugsr_1a';
filename{15} = 'hpugsr_1b';
filename{16} = 'hpugsr_2';
filename{17} = 'hpugsr_4a';
filename{18} = 'hpugsr_4b';
filename{19} = 'hpugsr_5';

filename{20} = 'hpugata_1' ;
filename{21} = 'hpugata_2' ;
filename{22} = 'hpugata_3' ;
filename{23} = 'hpugata_4' ;


nfile = length(filename);

lowcut = 4 ;%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 = 4;
grp{1} = 1:5  ; 
grp{2} = 6:13 ;
grp{3} = 14:19 ;
grp{4} = 20:23 ;

%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