Adding the stim from this window to the average

8 vues (au cours des 30 derniers jours)
Seungryul Lee
Seungryul Lee le 19 Déc 2022
Réponse apportée : Rik le 19 Déc 2022
Hi, I'm doing a project, and I am trying to keep a sum of the data in the windows that is added to in each iteration of the for loop, and then divide by the number of interations after the for loop. And it's for taking the average.
nPoints = 10000; % number of points in stimulus (stim) and response (resp)
stim = randn(1, nPoints); % stim is random noise drawn from a normal distribution
resp = stim > 2; % vector of spikes (when stim is > 2)
resp = [zeros(1,100), resp(1:end-100)]; % added 100 zeros in front of resp, to mimic latency that would exist in a real system (spikes are not instant)
subplot(2,1,1), plot(stim(1:1000))
line([1, 1000], [2, 2], 'LineStyle', '--', 'color', 'k')
ylabel('Mag, arb')
title('Stimulus')
subplot(2,1,2);
plot(resp(1:1000), 'r')
ylabel('Spikes')
title('Response')
xlabel('Time, pts')
windowSize = 300;
resp(1:windowSize)=0; % discard spikes in first window
nEvs = sum(resp); % binary vector, we can simply add to find number of events
evIdx = find(resp); % gives indexes of spikes (where a 1 was found)
% For each event
for w = 1: nEvs
% Find the indexes of the time window preceding the event
wId = evIdx(w)-windowSize : evIdx(w)-1;
% Add the stim from this window to the average
avg = avg + stim(wIdx);
end
'avg' is used in Include Code of User-Defined Supporting Functions in Report.
avg = avg./sum(resp);
figure
plot(avg);
title(['Average of', num2str(nEvs), 'windows']);
xlabel('Time, points');
ylabel('Mage, arb');
And I'm getting an error for
avg = avg + stim(wIdx);
this line.
Can I please get help for getting average, with the method that I'm trying to use?
Thank you.
%%REVISED
  5 commentaires
Rik
Rik le 19 Déc 2022
As you can see, you're getting a fairly cryptic error message which doesn't match what you posted.
I think the reason is that you don't actually define avg anywhere. If I'm wrong, can you point to where you do define it?
Seungryul Lee
Seungryul Lee le 19 Déc 2022
I think I have not defined.
If the average that I want is ave spike from all time windows, do you know the correct code for the average calculation? I'm having trouble with calculating the average.

Connectez-vous pour commenter.

Réponse acceptée

Rik
Rik le 19 Déc 2022
After you correct the variable name to wId, you end up with a vector of 300 elements. What avg should be isn't completely clear to me, but perhaps you intend it to initialize with zeros. At least the code will run without errors. You will have to check whether the code will indeed do what you expect.
I added a call to reshape to prevent dynamic expansion if the dimensions flip.
nPoints = 10000; % number of points in stimulus (stim) and response (resp)
stim = randn(1, nPoints); % stim is random noise drawn from a normal distribution
resp = stim > 2; % vector of spikes (when stim is > 2)
resp = [zeros(1,100), resp(1:end-100)]; % added 100 zeros in front of resp, to mimic latency that would exist in a real system (spikes are not instant)
subplot(2,1,1), plot(stim(1:1000))
line([1, 1000], [2, 2], 'LineStyle', '--', 'color', 'k')
ylabel('Mag, arb')
title('Stimulus')
subplot(2,1,2);
plot(resp(1:1000), 'r')
ylabel('Spikes')
title('Response')
xlabel('Time, pts')
windowSize = 300;
resp(1:windowSize)=0; % discard spikes in first window
nEvs = sum(resp); % binary vector, we can simply add to find number of events
evIdx = find(resp); % gives indexes of spikes (where a 1 was found)
% For each event
avg=zeros(1,windowSize);
for w = 1: nEvs
% Find the indexes of the time window preceding the event
wId = evIdx(w)-windowSize : evIdx(w)-1;
% Add the stim from this window to the average
avg = avg + reshape(stim(wId),1,windowSize);
end
avg = avg./sum(resp);
figure
plot(avg);
title(sprintf('Average of %d windows',nEvs));
xlabel('Time, points');
ylabel('Mage, arb');

Plus de réponses (0)

Catégories

En savoir plus sur Data Type Conversion dans Help Center et File Exchange

Tags

Produits


Version

R2022b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by