Custom processing for ABR weighted averaging

[mcp0228]

Hello
I am trying to use custom processing for ABR weighted averaging and weighted residual noise calculations as described by "Don and Elberling 1994 Evaluating residual background noise in ABRs".
The "weighted average" that already exists in Brainstorm is not what I am looking for.

The following is where I am up to:

1. I have imported 6000 epochs/files, each epoch going from 0 to 30 ms.
2. All the epochs are already preprocessed and ready for the next steps.

The following is what I want to do.
I attach an excel file of what I mean to further demonstrate what I've described:

1. Calculation of epoch weighting (column SF) = 1 / (Variance (epoch no. Y))
2. Classic residual noise calculation (column SD) =
   standard deviation (1st data point of epoch no. 1 to last data point of epoch no. Y) / square root (epoch no. Y)
3. Application of epoch weighting to data points (columns SI onwards) = all data points of epoch no. Y * epoch no. specific weighting (calculated in 1)
4. Weighted residual noise calculation (column SE) =
   (standard deviation (1st epoch weighted data point of epoch no. 1 (calculated in 3) to last epoch weighted data point of epoch no. Y (calculated in 3)) / square root (epoch no. Y)) / average (1st epoch specific weighting (calculated in 1) to no. Y epoch specific weighting (calculated in 1))
5. Calculation of weighted averaging = average (data points in 3).
6. Plot classic residual noise and weighted residual noise on the y-axis against the x-axis epoch no.
7. Plot classic averages (already on Brainstorm) and weighted averages on the y-axis against the x-axis of time in (ms).

In terms of coding this is where I am up to:

1. The "process" that is closest to what I am trying to do is "process_average.m".
2. I copied this file into "$HOME/.brainstorm/process" and renamed it "process_weighted_average.m".
3. I've changed "sProcess.Comment = 'Average files';" into " sProcess.Comment = 'Weighted averaging';"

I am really not sure what to do from here onwards.
What and where do I insert the relevant code so that the process described upwards works?
Any pointers would be greatly appreciated.

ABR weighted averaging and RN calculation.xlsx (151.9 KB)

[Raymundo.Cassani]

Hi,

The sort answer is: the "processing" code goes into the Run function of the process script.

As it was mentioned in a previous post answer, the place to look for information on how to write your own processes is:
https://neuroimage.usc.edu/brainstorm/Tutorials/TutUserProcess

While the process_average.m gives a similar result to your goal, to handle the many options in the process, there are extra functions besides the Run function. I'd suggest you to look at the simpler process file, process_example_customavg.m, for inspiration. This script is provided in the Examples section.

Best,
Raymundo

[mcp0228]

Thank you for that answer.
Where would I find codes for things like "variance", "standard deviation", "counting no. of epochs", "average", "square root" etc.?

[Francois]

Where would I find codes for things like "variance", "standard deviation", "counting no. of epochs", "average", "square root" etc.?

On the Mathworks website (Standard deviation - MATLAB std - MathWorks France), Wikipedia, or in the existing functions of Brainstorm. Of course, you need to know how to code in Matlab in order to do this.

[mcp0228]

Thank you Francois

I am now using the "Pre-process > Run Matlab command and managed to work out that the following code does what I need it to do for weighted epoching.

EpochW = ones(2,1)./var((Data),0,2);
Data = EpochW.*(Data);

Because the process category of Run Matlab command is "Filter", the process considers each file one at a time independently and saves each file. As I have 6000 epochs, after running through this process, I now have a new set of 6000 processed epochs which is what I was expecting.

Now an issue arose when doing averaged weighted epoching.
I need to take the sum of all the newly generated 6000 processed epochs giving me a single row vector containing the sum of each column (time). I then need to divide this by the sum of the outcome of the code "EpochW = ones(2,1)./var((Data),0,2)".

What would be the possible code for this?
Can Run Matlab command handle something like this when its process category is "Filter"?
I understand that process category "Custom" would work better in this situation as all the input files are passed at once to the process.

Thoughts?

[Francois]

As recommended initially by @Raymundo.Cassani, I think the best option is for you to write a custom process that does all of what you need at once, instead of creating 6000 intermediate and useless files.
Start from the example process_example_customavg.m.

[mcp0228]

@mcp0228

EpochW = ones(2,1)./var((Data),0,2);
Data = EpochW.*(Data);

Thank you for that Francois!
I managed to get everything working and successfully compared and contrasted its output with my VERY slow method of excel sheet.
Thankfully the outputs are exactly the same which is very reassuring.

The "Avg: 1 (6000)" is the original classic average and the WAvg is the weighted average.

How do I change the code so that Brainstorm writes the event number (in this case "1") and the number of files (6000)?

Many Thanks!

[Francois]

@mcp0228

How do I change the code so that Brainstorm writes the event number (in this case "1") and the number of files (6000)?

If you want the file to be displayed with a different in the database explorer, you need to change the Comment field of the data structures.

• https://neuroimage.usc.edu/brainstorm/Tutorials/Epoching#On_the_hard_drive
• https://neuroimage.usc.edu/brainstorm/Tutorials/Scripting

[mcp0228]

Thanks Francois!
I tried this and it worked well

I sometimes get number of epochs that are slightly different to 6000.
So rather than changing the "DataMat.Comment" field every time I get a different overall number of epoch, is there a way for Brainstorm to add as comment the number of epochs already counted?
In this particular case it would be "length(sInputs)"

However changing DataMat.Comment into 'WAvg length(sInputs)' only gives the name of the eventual file to be 'WAvg length(sInputs)' .

What can I do here?
Many Thanks!

[Francois]

'WAvg length(sInputs)'

Maybe you want to use something like:
DataMat.Comment = sprintf('WAvg (%d)', length(sInputs));

[mcp0228]

Thank you so much Francois!
Everything is working beautifully - I will now get onto the calculation of residual noises.

On a different note, I'd imagine there would be other researchers around the globe wanting to use weighted averaging as described by "Don and Elberling 1994 Evaluating residual background noise in ABRs" which is what my script does.

What would be the best way to share this?

[Francois]

What would be the best way to share this?

Writing a Brainstorm process, documenting it, and posting the code here:
bst-users/processes at master · brainstorm-tools/bst-users · GitHub

Or post it here if you don't know how to use github.

Thanks!

[mcp0228]

Hello Francois

I am now working on residual noise calculation and the following is what I have to do:

1. I realised that Brainstorm first loads the files into one large load matrix consisting of [Nchannels x Ntime x Nepochs]. I have 2 recording channels and each channel consists of 493 time data points and there are 6000 trials (or epochs). Therefore the 3D array that I'm working on has the size [2 x 493 x 6000] double. Let's call this 3D array "data".

2. To calculate residual noise I have to: calculate the standard deviation of the data points (in rows) from the 1st data point of epoch no.1 to the last data point of epoch no. X.
   The 1st data point of epoch no.1 can be given by data(:,1,1)

• this means all the rows (there are two recording channels), of time point 1 at array no. 1.
  The last data point of epoch no. X can be given by data(:,end,X)
• this means all the rows (two recording channels), of last time point (given by 'end') at array no. X.

Graphically this is what I'm trying to achieve. Each colour denotes the values used to calculate that std.

3. The problem that I've run into is how to specify using coding to calculate the std from epoch1 to epoch no X all the way to epoch no. 6000.
   I imagine the code would look something like the following.
   I've tried this but it does NOT work - this is just a demonstration.
   std(data,0,2,[data(:,1,1), data(:,end,1)]);
   std(data,0,2,[data(:,1,1), data(:,end,2)]);
   std(data,0,2,[data(:,1,1), data(:,end,3)]);
   std(data,0,2,[data(:,1,1), data(:,end,4)]);
   std(data,0,2,[data(:,1,1), data(:,end,5)]);
   .
   .
   .
   std(data,0,2,[data(:,1,1), data(:,end,6000)]);

Therefore my two questions are as follows:

1. What would be the code to specify std calculation from array 1 to array X in a 3D array?
2. What would be the best way to achieve this without writing 6000 lines?

[Raymundo.Cassani]

@mcp0228

1. What would be the code to specify std calculation from array 1 to array X in a 3D array?

These are two ways you could do it:

• If your Matlab R2018b or later, it's possible to indicate 'all' to use all the elements used in the std() function: Standard deviation - MATLAB std
• Otherwise, you can use reshape() to change the data from 2D [nEpochs, nTimes] to 1D [1, nEpochs×nTimes] and then use std()

@mcp0228

2. What would be the best way to achieve this without writing 6000 lines?

If an action is repetitive, it's a machine job. In this case you can use a for loop.

For example the operation you describe in the image can be coded as:
(for sake of the example only one channel is considered, thus data is [nEpochs, nTimes]

[nEpochs, nTimes] = size(data);
residualErr = zeros(1, nEpochs);
for iEpoch = 1 : nEpochs
    residualErr(iEpoch) = std(data(1:iEpoch, 1:nTimes), 0, 'all');
end

Regarding Matlab programming, you might want to check this book:
MATLAB for Brain and Cognitive Scientists by Mike X. Cohen.

[mcp0228]

Thank you Raymundo!

I've tried using your method of "all" but it's not working as I thought it would.
Therefore I've attached a clearer diagram of what I need my code to do.

Based on the below diagram I think the code should look something like this.
I've already tried this code and it does not work - I write this to hopefully demonstrate to you what I have in mind.

% for data in 3D = 6 (No. of channels) X 439 (No. of time points) X 6000 (No. of epochs)
Nepochs = 6000
for i = 1 : Nepochs
data(i) = std(data(:, :, 1:i), 3); % dimension value 3 added for std calculation over pages
end

Hopefully, the eventual code gives me data in 2D = 6 (No. of channels) X 6000 (No. of epochs) with each element being the std output as shown by the diagram below.

Hope this makes sense. What are your thoughts?

And additionally, I need my x-axis to show epoch number instead of time.
Therefore I've changed the definition of "Time" from

Time = DataMat.Time; (which gives "Time" in seconds)
to
Time = linspace (1,(Nepochs),(Nepochs)); (which gives epoch number as x-axis depending on the number of epochs).

However, the units still say "Time (s)" - how do I change this through coding?

[Francois]

Hope this makes sense. What are your thoughts?

Wishing you luck with the coding. It will be a good training to improve your Matlab programming skills

However, the units still say "Time (s)" - how do I change this through coding?

If you want to change the label from a figure programmatically after it has been created, you just need to execute xlabel('your label'); or xlabel(hAxes, 'your label');
However, there is no solution for saving a "data" or "matrix" file where the second dimension is not time. Therefore the figure will always show "Time (s)" when creating it.

When using view_image_reg, as illustrated in view_matrix(..., 'image'), you can specify the label of the two axes.

[mcp0228]

Thanks for that Francois!
Would you have any suggestions of where I might start to do the coding described by the image?
I'm quite stuck and I think I've exhausted all my ideas.
I'm willing to try anything that you can think of.
I've even asked the Matlab community the same question and currently waiting for answers!

[mcp0228]

So further to the point above I realised that I needed a vector dimension to specify that the calculation must be across the matrix then over the pages. That is to say...

RN = zeros(Nchannels, Nepochs); % Nchannels = 6 recording channels and Nepochs = 6000 epochs
for i = 1 : Nepochs
RN(i) = std(data(:, :, 1:i), 0, [2 3]); % data = 6 X 493 X 6000 3D matrix.
end

Here the std(A, 0, [2 3]) was used and [2 3] was specified as the calculation dimension is across the matrix first then over the pages.

However, I keep getting the error "Unable to perform assignment because the left and right sides have a different number of elements."

What am I doing wrong?

I thought by making "RN = zeros(Nchannels, Nepochs);" which is a 6 X 6000 matrix, the "for" loop would then automatically insert the output data (6 X 1 matrix) 6000 times eventually completing the 6 X 6000 matrix??

[Raymundo.Cassani]

@mcp0228

Would you have any suggestions of where I might start to do the coding described by the image?

This question is now of the scope of the forum, which is topics related to the use of Brainstorm rather than Matlab scripting. There are thousands of free resources online to learn the Matlab basics online.

Good luck with your Matlab learning path!

@mcp0228

What am I doing wrong?

Lastly, the error message is very clear:
"Unable to perform assignment because the left and right sides have a different number of elements."

and plenty information comes when googling it.

In short, the result of the right size is 6 elements, but the left side is a 1-element array. Thus the solution is to assign the 6-element result on the right to a 6-element array in the left:

RN(:, i) = std(data(:, :, 1:i), 0, [2 3]); % data = 6 X 493 X 6000 3D matrix.

[mcp0228]

Thank you very much for all your help Francois and Raymundo!
Everything works very well and the code is doing what I'm intending it to do (image of code attached).

I just have two more quick questions.
Q1. Currently, I've written the code so that it calculates classic (CRN) and weighted residual noises (WRN) separately each generating 6 X 6000 double matrix (relevant code attached below).
I want Brainstorm to display both classic and weighted noises at the same time.

Therefore I've concatenated the two matrices into one - AllMat = [CRN; WRN];
And as expected, this did generate a 12 X 6000 double matrix according to "View file contents" (view attached image below).
However, once I open the time series it does not give me the output for both types of data.
I still only see output for a single 6 X 6000 data, not both as I thought it would.
I thought concatenating both data sets in this form should work because they have the same sized matrices.

What can I do?

Q2. When I review the file names of the output from the code it names my files as 'data_custom_'.
I realised that this is because of the following code.

OutputFiles{1} = bst_process('GetNewFilename', fileparts(sInputs(1).FileName), 'data_custom_');

How do I change this code so that it reflects what I want the file to be called?
It's getting a bit confusing when all my custom files are being called 'data_custom_'

I already tried renaming the code into;
OutputFiles{1} = bst_process('GetNewFilename', fileparts(sInputs(1).FileName), 'WAvg');