view userProgramsTim/poolIPI_across_channels.m @ 38:c2204b18f4a2 tip

End nov big change

function [pooledIPIdata,ctr] = poolIPI_across_channels(IPIhisttime,IPIhistweight)

%function that pools IPIdata across all channels in order to plot the data
%as in Secker-Walker JASA 1990 Fig. 6

%plots IPI-Histograms in a matrix display as in Secker-Walker JASA 1990
%Fig. 4
%
%   Tim Jürgens, February 2011
%
% input:  IPIhisttime:     matrix containing the interval times found in the
%                          analysis
%                          first dimension: frequency channel
%                          second dimension: time step (hop position)
%                          third dimension: interval dimension (max 3)
%         IPIhistweight:   matrix containing the interval weights found in the
%                          analysis, same dimensions as IPIhisttime
% output: pooledIPIdata:   matrix containing the pooled and weighted IPI
%                          histograms as a function of time (hop position)
%                          first dimension: time step (hop position)                           
%                          second dimension: IPI intervals
%         ctr:             center values of classes of IPI intervals
% ATTENTION: depending on the amplitude of the signals the value
% 'verticalshift' might have to be adjusted in order to see structures in
% the plot

verticalshift =   3000; %0.008;%vertical shift of single time series in z-coordinate units
ctr = [0:0.02:5].*1e-3; %classes (center-values) of 20 microsec width

%preallocation of variables
pooledIPIdata = zeros(size(IPIhisttime,2),length(ctr));
smoothed_pooledIPI = zeros(size(IPIhisttime,2),length(ctr)-5);

for iCounter = 1:size(IPIhisttime,2) %one for time spacing
    %cannot use matlabs hist function because weighting must be applied
   
    tmpIPIhisttime = squeeze(IPIhisttime(:,iCounter,:));
    tmpIPIhistweight = squeeze(IPIhistweight(:,iCounter,:));
    tmpIPIhisttime = tmpIPIhisttime(:);
    tmpIPIhistweight = tmpIPIhistweight(:);
    for jCounter = 1:length(tmpIPIhisttime)
        %look which class
        [tmp1,classindex] = min(abs(tmpIPIhisttime(jCounter)-ctr));
        pooledIPIdata(iCounter,classindex) = pooledIPIdata(iCounter,classindex)+tmpIPIhistweight(jCounter);
    end
    
end

%smooth data using a 5-point hamming window
hamm_window = hamming(5);
for iCounter = 1:size(pooledIPIdata,1)
    for jCounter = 3:length(ctr)-2 %start with 3 and end with 2 samples
        %less the length of ctr in order not to get in conflict with the length of
        %the hamm_window
        smoothed_pooledIPI(iCounter,jCounter-2) = ...
            pooledIPIdata(iCounter,(jCounter-2):(jCounter+2))*hamm_window./sum(hamm_window);
    end
end
smoothed_ctr = ctr(3:end-2);


figure;

Tickvector = [];
TickLabels = [];
for iCounter = 1:size(IPIhisttime,2)
    hold on;
    verticalposition = verticalshift*(iCounter-1);
    %multiply by 1000 to set abscissa to ms units
    plot(1000.*smoothed_ctr, ...
        smoothed_pooledIPI(iCounter,:)+verticalposition, ...
        'k','LineWidth',2);  
    if mod(iCounter,5) == 1 %set best frequency as a label every 10 channels
        Tickvector = [Tickvector verticalposition];
        TickLabels = [TickLabels; (iCounter-1)*3]; %time spacing is every 3ms
    end
end
set(gca,'yTick',Tickvector,'yTickLabel',num2str(TickLabels,'%4.0f'));
ylabel('Stimulus Time (ms)');
xlabel('Interval (ms)');
xlim([min(1000*ctr) max(1000*ctr)]);
ylim([-verticalshift verticalposition+3*verticalshift]);
box on;