rmeddis@0: function [PH, binTimes]=UTIL_periodHistogram(A, dt, frequency)
rmeddis@0: % UTIL_makePeriodHistogram converts a time sequence into a period histogram
rmeddis@0: %*********************
rmeddis@0: % The period is 1/frequency.
rmeddis@0: %
rmeddis@0: % usage:
rmeddis@0: %	PH=UTIL_periodHistogram(A, dt, frequency)
rmeddis@0: %
rmeddis@0: % Input:
rmeddis@0: % A is a channel x time matrix of spikes (or other stuff)
rmeddis@0: % frequency determines the period of the histogram
rmeddis@0: %
rmeddis@0: % Output
rmeddis@0: %  PH is a channel by periodhistogram matrix
rmeddis@0: % bintimes is useful for plotting the output
rmeddis@0: %
rmeddis@0: 
rmeddis@0: periodInSeconds=1/frequency;
rmeddis@0: [numChannels signalNpoints]=size(A);
rmeddis@0: 
rmeddis@0: % retrict data array to a multiple of the period.
rmeddis@0: pointsPerPeriod= round(periodInSeconds/dt);
rmeddis@0: NcompletePeriods=floor(signalNpoints/pointsPerPeriod);
rmeddis@0: totalPointsUsed=NcompletePeriods*pointsPerPeriod;
rmeddis@0: 
rmeddis@0: % check that the period is a whole number of epochs
rmeddis@0: aliasing=NcompletePeriods*(periodInSeconds/dt-pointsPerPeriod);
rmeddis@0: 
rmeddis@0: if aliasing>.1
rmeddis@0:     error('UTIL_periodHistogram: irregular period length')
rmeddis@0: end
rmeddis@0: 
rmeddis@0: if NcompletePeriods<1
rmeddis@0:     error('UTIL_periodHistogram: too few datapoints')
rmeddis@0: end
rmeddis@0: 
rmeddis@0: % transpose data so that time is down a column
rmeddis@0: A=A(:,1:totalPointsUsed)';
rmeddis@0: 
rmeddis@0: % knock it into shape
rmeddis@0: A=reshape(A,pointsPerPeriod, NcompletePeriods, numChannels);
rmeddis@0: % each period is a separate column
rmeddis@0: % imagesc(squeeze(A)) % should have horizontal stipe
rmeddis@0: 
rmeddis@0: % channels are now the third dimension.
rmeddis@0: PH=squeeze(sum(A,2))';
rmeddis@0: % PH=PH/NcompletePeriods;
rmeddis@0: 
rmeddis@0: binTimes=dt:dt:pointsPerPeriod*dt;