This commit is contained in:
Jan Grewe 2014-10-30 17:46:29 +01:00
parent 40c847e694
commit 730e819abe
3 changed files with 21 additions and 16 deletions

View File

@ -1,11 +1,11 @@
function [sta, std_sta, valid_spikes]= sta(stimulus, spike_times, count, sampling_rate) function [st_avg, std_sta, valid_spikes]= sta(stimulus, spike_times, count, sampling_rate)
snippets = zeros(numel(spike_times), 2*count); snippets = zeros(numel(spike_times), 2*count);
valid_spikes = 1; valid_spikes = 1;
for i = 1:numel(spike_times) for i = 1:numel(spike_times)
t = spike_times(i); t = spike_times(i);
index = round(t*sampling_rate); index = round(t*sampling_rate);
if index < count || (index + count) > length(stimulus) if index <= count || (index + count) > length(stimulus)
continue continue
end end
snippets(valid_spikes,:) = stimulus(index-count:index+count-1); snippets(valid_spikes,:) = stimulus(index-count:index+count-1);
@ -14,5 +14,5 @@ end
snippets(end-(end-valid_spikes):end,:) = []; snippets(end-(end-valid_spikes):end,:) = [];
sta = mean(snippets, 1); st_avg = mean(snippets, 1);
std_sta = std(snippets,[],1); std_sta = std(snippets,[],1);

View File

@ -12,24 +12,24 @@ for i = 1:length(spike_times)
all_times = cat(1, all_times, spike_times{i}); all_times = cat(1, all_times, spike_times{i});
end end
[sta, sta_sd, num] = sta(stimulus_strong(:,2), all_times, 1000, sample_rate); [st_average, sta_sd, num] = sta(stimulus(:,2), all_times, 1000, sample_rate);
fig = figure(); fig = figure();
set(fig, 'PaperUnits', 'centimeters'); set(fig, 'PaperUnits', 'centimeters');
set(fig, 'PaperSize', [11.7 9.0]); set(fig, 'PaperSize', [11.7 9.0]);
set(fig, 'PaperPosition',[0.0 0.0 11.7 9.0]); set(fig, 'PaperPosition',[0.0 0.0 11.7 9.0]);
set(fig,'Color', 'white') set(fig,'Color', 'white')
plot(((1:length(sta))-1000)/sample_rate, sta) plot(((1:length(st_average))-1000)/sample_rate, st_average)
xlabel('time [s]') xlabel('time [s]')
ylabel('stimulus') ylabel('stimulus')
%% reverse reconstruction %% reverse reconstruction
% make binary representation of the spike times % make binary representation of the spike times
binary_spikes = zeros(size(stimulus_strong, 1), length(spike_times)); binary_spikes = zeros(size(stimulus, 1), length(spike_times));
estimated_stims = zeros(size(binary_spikes)); estimated_stims = zeros(size(binary_spikes));
for i = 1:length(spike_times) for i = 1:length(spike_times)
binary_spikes(round(spike_times{i}*sample_rate), i) = 1; binary_spikes(round(spike_times{i}*sample_rate), i) = 1;
estimated_stims(:,i) = conv(binary_spikes(:,i), sta, 'same'); estimated_stims(:,i) = conv(binary_spikes(:,i), st_average, 'same');
end end
fig = figure(); fig = figure();
@ -37,9 +37,9 @@ set(fig, 'PaperUnits', 'centimeters');
set(fig, 'PaperSize', [11.7 9.0]); set(fig, 'PaperSize', [11.7 9.0]);
set(fig, 'PaperPosition',[0.0 0.0 11.7 9.0]); set(fig, 'PaperPosition',[0.0 0.0 11.7 9.0]);
set(fig,'Color', 'white') set(fig,'Color', 'white')
plot(stimulus_strong(:,1), stimulus_strong(:,2), 'displayname','original') plot(stimulus(:,1), stimulus(:,2), 'displayname','original')
hold on hold on
plot(stimulus_strong(:,1), mean(estimated_stims,2), 'r', 'displayname', 'reconstruction') plot(stimulus(:,1), mean(estimated_stims,2), 'r', 'displayname', 'reconstruction')
xlabel('time [s]') xlabel('time [s]')
ylabel('stimulus') ylabel('stimulus')
legend show legend show
@ -49,14 +49,19 @@ legend show
% we need to downsample the data otherwise the covariance matrixs gets too % we need to downsample the data otherwise the covariance matrixs gets too
% large 20Khz to 1kHz % large 20Khz to 1kHz
downsampled_binary = zeros(size(stimulus_strong, 1)/20, length(spike_times)); % downsampled_binary = zeros(size(stimulus, 1)/20, length(spike_times));
downsampled_stim = zeros(size(downsampled_binary,1),1); downsampled_stim = zeros(size(stimulus,1)/20,1);
for i = 1:length(spike_times) % for i = 1:length(spike_times)
binary_spikes(round(spike_times{i}*1000), i) = 1; % indices = round(spike_times{i}.*1000);
end % indices(indices < 1) = [];
% downsampled_binary(indices, i) = 1;
% end
for i = 1:length(downsampled_stim) for i = 1:length(downsampled_stim)
start_index = (i-1) * 1000 + 1; start_index = (i-1) * 20 + 1;
downsampled_stim(i) = mean(stimulus_strong()*20)) downsampled_stim(i) = mean(stimulus(start_index:start_index+19,2));
end end
[st_average, ~, ~] = sta(downsampled_stim, all_times, 50, 1000);