rmeddis@0
|
1 function [PH, binTimes]=UTIL_periodHistogram(A, dt, frequency)
|
rmeddis@0
|
2 % UTIL_makePeriodHistogram converts a time sequence into a period histogram
|
rmeddis@0
|
3 %*********************
|
rmeddis@0
|
4 % The period is 1/frequency.
|
rmeddis@0
|
5 %
|
rmeddis@0
|
6 % usage:
|
rmeddis@0
|
7 % PH=UTIL_periodHistogram(A, dt, frequency)
|
rmeddis@0
|
8 %
|
rmeddis@0
|
9 % Input:
|
rmeddis@0
|
10 % A is a channel x time matrix of spikes (or other stuff)
|
rmeddis@0
|
11 % frequency determines the period of the histogram
|
rmeddis@0
|
12 %
|
rmeddis@0
|
13 % Output
|
rmeddis@0
|
14 % PH is a channel by periodhistogram matrix
|
rmeddis@0
|
15 % bintimes is useful for plotting the output
|
rmeddis@0
|
16 %
|
rmeddis@0
|
17
|
rmeddis@0
|
18 periodInSeconds=1/frequency;
|
rmeddis@0
|
19 [numChannels signalNpoints]=size(A);
|
rmeddis@0
|
20
|
rmeddis@0
|
21 % retrict data array to a multiple of the period.
|
rmeddis@0
|
22 pointsPerPeriod= round(periodInSeconds/dt);
|
rmeddis@0
|
23 NcompletePeriods=floor(signalNpoints/pointsPerPeriod);
|
rmeddis@0
|
24 totalPointsUsed=NcompletePeriods*pointsPerPeriod;
|
rmeddis@0
|
25
|
rmeddis@0
|
26 % check that the period is a whole number of epochs
|
rmeddis@0
|
27 aliasing=NcompletePeriods*(periodInSeconds/dt-pointsPerPeriod);
|
rmeddis@0
|
28
|
rmeddis@0
|
29 if aliasing>.1
|
rmeddis@0
|
30 error('UTIL_periodHistogram: irregular period length')
|
rmeddis@0
|
31 end
|
rmeddis@0
|
32
|
rmeddis@0
|
33 if NcompletePeriods<1
|
rmeddis@0
|
34 error('UTIL_periodHistogram: too few datapoints')
|
rmeddis@0
|
35 end
|
rmeddis@0
|
36
|
rmeddis@0
|
37 % transpose data so that time is down a column
|
rmeddis@0
|
38 A=A(:,1:totalPointsUsed)';
|
rmeddis@0
|
39
|
rmeddis@0
|
40 % knock it into shape
|
rmeddis@0
|
41 A=reshape(A,pointsPerPeriod, NcompletePeriods, numChannels);
|
rmeddis@0
|
42 % each period is a separate column
|
rmeddis@0
|
43 % imagesc(squeeze(A)) % should have horizontal stipe
|
rmeddis@0
|
44
|
rmeddis@0
|
45 % channels are now the third dimension.
|
rmeddis@0
|
46 PH=squeeze(sum(A,2))';
|
rmeddis@0
|
47 % PH=PH/NcompletePeriods;
|
rmeddis@0
|
48
|
rmeddis@0
|
49 binTimes=dt:dt:pointsPerPeriod*dt;
|