Filter löschen
Filter löschen

Adding the stim from this window to the average

16 Ansichten (letzte 30 Tage)
Seungryul Lee
Seungryul Lee am 19 Dez. 2022
Beantwortet: Rik am 19 Dez. 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 Kommentare
Rik
Rik am 19 Dez. 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 am 19 Dez. 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.

Melden Sie sich an, um zu kommentieren.

Akzeptierte Antwort

Rik
Rik am 19 Dez. 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');

Weitere Antworten (0)

Kategorien

Mehr zu Electrophysiology finden Sie in Help Center und File Exchange

Tags

Produkte


Version

R2022b

Community Treasure Hunt

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

Start Hunting!

Translated by