ssa_rebound_rasters_v2.m 4.0 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123
  1. function ssa_rebound_rasters_v2(spikeTimes, alignmentTimes, unitData, timeWindow, micData, recDate, uType)
  2. %%%%
  3. % this plots a figure showing unit response to different perturbations
  4. %%%%
  5. % spike times is a simple vector of times when spike occurs (sec)
  6. % alignment Times is cell array where each entry contains a vector of start
  7. % times for all perturbations of a given type
  8. % unitData is unit #
  9. % window is time (sec) before and after 0 for each plot
  10. % micData is microphone data, used to generate spectrogram
  11. %%%%
  12. % function generates single figure window:
  13. % - A DIFFERENT COLUMN FOR EACH PERTURBATION TYPE! THIS COULD BE AS MANY AS
  14. % 19 COLUMNS! (but only 13 for first try because we messed up
  15. % perturbations)
  16. % - each column has spectrogram at top, raster in middle, psth at bottom
  17. %%%%
  18. figPath = 'D:\ssa_expmts\figures'; % UPDATE THIS
  19. fs = 40000;
  20. these_ATs = alignmentTimes;
  21. this_window = timeWindow;
  22. micRows = zeros(length(these_ATs),round((this_window(2)-this_window(1))*fs)+1);
  23. for trialIdx = 1:length(these_ATs)
  24. thisIdxRange = round((these_ATs(trialIdx)+this_window)*fs);
  25. micRows(trialIdx,:) = micData(thisIdxRange(1):thisIdxRange(2))';
  26. end
  27. avgMic = micRows(1,:);
  28. winLen = round(.005*fs); % 1ms window on spec
  29. freqs = 1:10:8001; % step to 8kHz in 10Hz bins
  30. figure(18);
  31. subplot(4,1,1)
  32. spectrogram(avgMic, winLen, floor(winLen/2), freqs, fs, 'yaxis');
  33. % this stuff to be inserted with actual plotting
  34. lim = caxis;
  35. span = diff(caxis);
  36. newspan = span/2;
  37. newlim1 = lim(2) - newspan;
  38. caxis([newlim1 lim(2)])
  39. set(gca, 'Xticklabel', []);
  40. set(gca, 'Yticklabel', []);
  41. title(['Unit ' num2str(unitData)])
  42. xlabel([])
  43. ylabel([])
  44. colorbar('off')
  45. %% now do PSTHs and ssa idx
  46. these_ATs = alignmentTimes;
  47. [~, bins, ~, ~, ~, ba] = psthAndBA(spikeTimes, these_ATs, timeWindow, 0.001); %1ms bin size
  48. [tr,b] = find(ba);
  49. [rasterX,yy] = rasterize(bins(b));
  50. rasterY = yy+reshape(repmat(tr',3,1),1,length(tr)*3); % yy is of the form [0 1 NaN 0 1 NaN...] so just need to add trial number to everything
  51. % scale the raster ticks
  52. rasterY(2:3:end) = rasterY(2:3:end)+floor(numel(these_ATs)/100);
  53. % PSTH smoothing filter
  54. % gw = gausswin(round(*6),3); % 5 ms gaussian window
  55. % smWin = gw./sum(gw);
  56. %
  57. % % smooth ba
  58. % baSm = conv2(smWin,1,ba', 'same')'./.001;
  59. psthSm = mean(ba); % mean(baSm) when smoothing
  60. figure(18);
  61. subplot(4,1,[2 3 4])
  62. plot(rasterX, rasterY, 'k', 'MarkerSize', 1);
  63. xlim(timeWindow);
  64. xlabel('Time (s)')
  65. ylim([1 length(these_ATs)+1]);
  66. % makepretty;
  67. % if pertIdx ~= 1
  68. % pTime = pert_time(pertIdx);
  69. % pLen = pert_lens(pertIdx);
  70. % hold on;
  71. % lim_of_y = ylim;
  72. % if songIdx == 1
  73. % yl = [lim_of_y(1) 75+1];
  74. % else
  75. % yl = [lim_of_y(1) 50+1];
  76. % end
  77. % v1 = [pTime yl(1); pTime yl(2); pTime+pLen yl(1); pTime+pLen yl(2)];
  78. % f1 = [1 2 4 3];
  79. %
  80. % patch('Faces', f1, 'Vertices', v1, 'FaceColor', 'k', 'FaceAlpha', 0.25);
  81. % if songIdx == 1
  82. % yl = [75 lim_of_y(2)];
  83. % else
  84. % yl = [50 lim_of_y(2)];
  85. % end
  86. %
  87. % v1 = [pTime yl(1); pTime yl(2); pTime+pLen yl(1); pTime+pLen yl(2)];
  88. % f1 = [1 2 4 3];
  89. %
  90. % patch('Faces', f1, 'Vertices', v1, 'FaceColor', 'k', 'FaceAlpha', 0.25);
  91. %
  92. set(gca,'YTick',[])
  93. h1 = refline(0,6);
  94. h1.Color = 'k';
  95. h2 = refline(0,16);
  96. h2.Color = 'k';
  97. h3 = refline(0,11);
  98. h3.Color = 'r';
  99. hold off;
  100. % end
  101. ylabel(['event number'])
  102. xlabel(['Time (mot. on, s)'])
  103. % if pertIdx == 1
  104. % % ylabel(['Avg. Firing Rate (Hz)'])
  105. % end
  106. xlim(timeWindow)