Archive

Archive for the ‘matlab’ Category

We contributed to MatLab (wavelet toolbox)

January 25th, 2017

We use MatLab a lot! It’s the major program for brain imaging data analysis in our lab. However, I never thought we could actually contribute to MatLab’s development.

In MatLab 2016, there is a toolbox called Wavelet Toolbox. If you read the document on wavelet coherence (link below), you will find that they used our NIRS data as an example:

https://www.mathworks.com/help/wavelet/examples/compare-time-frequency-content-in-signals-with-wavelet-coherence.html

Back in 2015/4/9, Wayne King from MathWorks contacted us, saying that they are developing the wavelet toolbox and asking if we can share some data as an example. We did. I’m glad that it’s part of the package now.

The following section are from the page above:


Find Coherent Oscillations in Brain Activity

In the previous examples, it was natural to view one time series as influencing the other. In these cases, examining the lead-lag relationship between the data is informative. In other cases, it is more natural to examine the coherence alone.

For an example, consider near-infrared spectroscopy (NIRS) data obtained in two human subjects. NIRS measures brain activity by exploiting the different absorption characteristics of oxygenated and deoxygenated hemoglobin. The data is taken from Cui, Bryant, & Reiss (2012) and was kindly provided by the authors for this example. The recording site was the superior frontal cortex for both subjects. The data is sampled at 10 Hz. In the experiment, the subjects alternatively cooperated and competed on a task. The period of the task was approximately 7.5 seconds.

load NIRSData;
figure
plot(tm,NIRSData(:,1))
hold on
plot(tm,NIRSData(:,2),'r')
legend('Subject 1','Subject 2','Location','NorthWest')
xlabel('Seconds')
title('NIRS Data')
grid on;
hold off;

Obtain the wavelet coherence as a function of time and frequency. You can use wcoherence to output the wavelet coherence, cross-spectrum, scale-to-frequency, or scale-to-period conversions, as well as the cone of influence. In this example, the helper function helperPlotCoherence packages some useful commands for plotting the outputs of wcoherence.

[wcoh,~,F,coi] = wcoherence(NIRSData(:,1),NIRSData(:,2),10,'numscales',16);
helperPlotCoherence(wcoh,tm,F,coi,'Seconds','Hz');

In the plot, you see a region of strong coherence throughout the data collection period around 1 Hz. This results from the cardiac rhythms of the two subjects. Additionally, you see regions of strong coherence around 0.13 Hz. This represents coherent oscillations in the subjects’ brains induced by the task. If it is more natural to view the wavelet coherence in terms of periods rather than frequencies, you can use the ‘dt’ option and input the sampling interval. With the ‘dt’ option, wcoherence provides scale-to-period conversions.

[wcoh,~,P,coi] = wcoherence(NIRSData(:,1),NIRSData(:,2),seconds(0.1),...
    'numscales',16);
helperPlotCoherence(wcoh,tm,seconds(P),seconds(coi),...
    'Time (secs)','Periods (Seconds)');

Again, note the coherent oscillations corresponding to the subjects’ cardiac activity occurring throughout the recordings with a period of approximately one second. The task-related activity is also apparent with a period of approximately 8 seconds. Consult Cui, Bryant, & Reiss (2012) for a more detailed wavelet analysis of this data.

Conclusions

In this example you learned how to use wavelet coherence to look for time-localized coherent oscillatory behavior in two time series. For nonstationary signals, it is often more informative if you have a measure of coherence that provides simultaneous time and frequency (period) information. The relative phase information obtained from the wavelet cross-spectrum can be informative when one time series directly affects oscillations in the other.

References

Cui, X., Bryant, D.M., and Reiss. A.L. “NIRS-Based hyperscanning reveals increased interpersonal coherence in superior frontal cortex during cooperation”, Neuroimage, 59(3), pp. 2430-2437, 2012.

Grinsted, A., Moore, J.C., and Jevrejeva, S. “Application of the cross wavelet transform and wavelet coherence to geophysical time series”, Nonlin. Processes Geophys., 11, pp. 561-566, 2004.

Maraun, D., Kurths, J. and Holschneider, M. “Nonstationary Gaussian processes in wavelet domain: Synthesis, estimation and significance testing”, Phys. Rev. E 75, pp. 016707(1)-016707(14), 2007.

Torrence, C. and Webster, P. “Interdecadal changes in the ESNO-Monsoon System,” J.Clim., 12, pp. 2679-2690, 1999.

Author: Xu Cui Categories: brain, matlab, nirs, programming Tags:

Communications between two MatLabs (2): over socket

October 17th, 2016

Aaron Piccirilli

Aaron Piccirilli

After the previous blog post Communications between two MatLabs (1) over file, Aaron Piccirilli in our lab suggested a more efficient way to communicate between two matlabs, i.e. over socket. Below is the source code provided by Aaron:

udpSocket = udp('127.0.0.1', 'LocalPort', 2121, 'RemotePort', 2122);
fopen(udpSocket);
udpCleaner = onCleanup(@() fclose(udpSocket));
for ii = 1:100
    fprintf(udpSocket, '%d%f', [ii ii/100]);
    disp(['Sending ' num2str(ii)]);
    pause(0.1);
end
udpSocket = udp('127.0.0.1', 'LocalPort', 2122, 'RemotePort', 2121);
fopen(udpSocket);
udpCleaner = onCleanup(@() fclose(udpSocket));
    while(1)
    if udpSocket.BytesAvailable
        ii = fscanf(udpSocket, '%d%f');
        disp(['Received' num2str(ii(1)) '-' num2str(ii(2))])
    end
end

More words from Aaron:

I also wrote a couple of quick tests to compare timing by having each method pass 10,000 random integers as quickly as they could. Using UDP is over four times faster on my work machine, and would be sufficient to keep up with sampling rates up to about 900 Hz, whereas the file-based transfer became too slow at about 200 Hz.

Obviously these rates and timings are going to be system and data-dependent, but the UDP implementation is about the same amount of code. It has some added benefits, too. First is what I mentioned before - that this allows you to communicate between different languages. Second, though, is what might be more important: buffer management. If your data source is sending data faster than you can process it, even for just a moment, the UDP method handles that gracefully with a buffer. To get the same functionality with the file method you have to write your own buffer maintenance - not too tricky, but adds another layer of complexity and probably won’t be as efficient.

I did a similar timing test passing 40 floats each time (say for 20 channels of NIRS data) instead of a single integer and the timing did not really change on my machine.

I also tested the above scripts, and they work beautifully! I definitely recommend this method over the ‘file’ method. One thing to note: when you Ctrl-C to quit the program, remember to close the socket (fclose(udpSocket)) AND clean the variables (udpSocket, udpCleaner); otherwise you will run into the “Unsuccessful open: Unrecognized Windows Sockets error: 0: Cannot bind” error.

Note from Aaron:

One note: the onCleanup function/object is designed as a callback of sorts: no matter how the function exits (normally, error, crash, Ctrl-C), when the onCleanup object is automatically then destroyed, its anonymous function should run. Thus, the UDP connection should be closed no matter how you exit the function. This won’t work for a script, though, or if you were just running the code on its own in a command window, as the onCleanup object wouldn’t be automatically destroyed. I would just exclude that line completely if you weren’t running it as a function.

Author: Xu Cui Categories: matlab, programming Tags:

Communications between two MatLabs (1) over file

October 3rd, 2016

Ref to: Communications between two MatLabs (2): over socket

It’s common that two MatLab programs needs to communicate. For instance, one program is collecting the brain imaging data but not display them, and the other program is to display the data. (Another case is at http://www.alivelearn.net/?p=1265) Sometimes it is not practical to merge the two program together (e.g. to keep the code clean). In this case we can run two MatLabs simultaneously. One keeps saving the data to a file, and the other keep reading the file.

Here I played with such a setup, and find they communicate well with small delay (small enough for hemodynamic responses). Check out the video below:

writeSomething.m

for ii=1:100
    save('data','ii');
    disp(['write ' num2str(ii)])
    pause(1)
end
readSomething.m

last_ii = 0;
while(1)
    try
        load data
        if(ii ~= last_ii)
            disp(['get data. i=' num2str(ii)])
        end
        last_ii = ii;
    end
    pause(0.1)
end

Caveat: writing/reading to/from disc is slow. So if your experiment requires real time communication without any delay (say <1ms), this method may not work. Also, the amount of data to write/read each time should be very small, and the frequency of write should be small too. The file needs to locate in your local hard drive instead of a network drive.

———- Comments ———–
Paul Mazaika from Stanford:
Cool piece of code! There may be a way to do this with one umbrella Matlab program that calls both components as subroutines. The potential advantage is that one program will keep memory in cache, not at disk, which can support rapidly updating information. For high speeds, it may be better to only occasionally update the graphical display, which otherwise may be a processing bottleneck.
-Paul

Aaron Piccirilli from Stanford:
There is, sort’ve! I think Xu’s little nugget is probably best choice for many applications, but if speed is an especially big concern then there are a couple of options that I’ve come across that will maintain some sort of shared memory.

Perhaps the easiest is to use sockets to communicate data, via UDP or TCP/IP, just like you use over the internet, but locally. You write some data to a socket in one program, and read it from that same socket in another program. This keeps all of your data in memory as opposed to writing it to disk, but there is definitely some overhead for housekeeping and to move the data from one program’s memory into the operating system’s memory then back into the other program’s memory. An added bonus here: you can communicate between different languages. If you have a logging function written in Python and a visualization program in MATLAB, they can pretty easily communicate with each other via sockets.

MATLAB doesn’t have explicit parallel computing built-in like many other languages, sadly, but we all have access here to the Parallel Computing Toolbox, which is another option for some more heavy-duty parallel processing where you have a problem you can easily distribute to multiple workers.

Finally, true shared memory might be more trouble than it’s worth for most applications, as you then have to deal with potential race conditions of accessing the same resource at the same time.

Aaron

More on this topic: Please continue to read Communications between two MatLabs (2): over socket

Author: Xu Cui Categories: brain, matlab, nirs, programming Tags:

A mistake in my False discovery rate (FDR) correction script

August 8th, 2016

I have posted an FDR script at http://www.alivelearn.net/?p=1840. I noticed that there is a small bug. In rare cases, this bug will cause the most significant voxel to be classified as ‘non-significant’ while other voxels are ’significant’.

Consider the following example:

p = [0.8147 0.9058 0.0030 0.9134 0.6324 0.0029 0.2785 0.5469 0.9575 0.9649 0.1576 0.9706 0.9572 0.4854 0.8003 0.1419 0.4218 0.9157];

The previous script will classify p(3) as significant but p(6) as non-significant.

Here is the updated version of the script:

function y = fdr0(p, q)
% y = fdr0(p, q)
%
% to calculate whether a pvalue survive FDR corrected q
%
% p: an array of p values. (e.g. p values for each channel)
% q: desired FDR threshold (typically 0.05)
% y: an array of the same size with p with only two possible values. 0
% means this position (channel) does not survive the threshold, 1 mean it
% survives
%
% Ref:
% Genovese et al. (2002). Thresholding statistical maps in functional
% neuroimaging using the false discovery rate. Neuroimage, 15:722-786.
%
% Example:
%   y = fdr0(rand(10,1),0.5);
%
% Xu Cui
% 2016/3/14
%

pvalue = p;
y = 0 * p;

[sortedpvalue, sortedposition] = sort(pvalue);
v = length(sortedposition);
for ii=1:v
    if q*ii/v >= sortedpvalue(ii)
        y(sortedposition(1:ii)) = 1;
    end
end

return;
Author: Xu Cui Categories: brain, matlab, nirs, programming Tags:

MatLab tip, save data in compressed mode

June 29th, 2016

I have a mat file which is 100+M in size. It would be fine normally but since I was trying to upload it to github, it was rejected due to its big size. Fortunately there is a way to reduce the size: save it again in the compressed mode.

After I load the file, I save the variable again using the following command:

save filename variablename -v7

The parameter -v7 is the key. How does the compression perform? Originally the file (render_ch2bet) size is 103M, with -v7 parameter the file is 40.9M (render_ch2bet_compressed). By comparison, if we compress the original file using a third party program, we got ~40M (render_ch2bet.zip). If we use the -v7 option, not only we get a much smaller file and save a lot of space, we can also load it directly in MatLab.

Compressed MAT file

Author: Xu Cui Categories: matlab Tags:

Some tips to use wavelet toolbox

June 13th, 2016

Wavelet toolbox is a useful tool to study hyperscanning data. Many recent publications on NIRS hyperscanning use wavelet coherence to quantify the relationship between two interacting brains (e.g. Baker et al 2016, Nozawa et al 2016). You can see more information about wavelet coherence at http://www.alivelearn.net/?p=1169

wavelet

wavelet

Here are some tips to use the toolbox:

1. It often takes a long time to run Monte Carlo simulation. You may use ‘mcc’=0 to disable it.

figure;wtc(signal(:,jj),signal(:,jj+22),'mcc',0);

2. If you need to get the values of the result (instead of the graphic), you may specify the return value

[Rsq,period,scale,coi,sig95] = wtc(signal1,signal2,'mcc',0); %Rsq is a complex matrix

3. If you are only interested in a certain band, you can specify the MaxScale (i.e. ms) parameter. More information at http://www.alivelearn.net/?p=1518

[Rsq,period,scale,coi,sig95] = wtc(signal1,signal2,'mcc',0, 'ms', 128);

4. If you are interested in finding the “phase” information (visualized by the arrows), you may use xwt function. The returned value is a complex matrix and you can calculate the phase.

[Rsq,period,scale,coi,sig95] = xwt(signal(:,jj),signal(:,jj+22));
figure;imagesc(angle(Rsq));

5. To visualize the power of a single signal, you may use wt, which I personally feel much better than FFT.

figure;wt(signal(:,jj));

6. To change the density of the arrows, you may specify the ArrowDensity parameter

figure;wtc(signal(:,jj),signal(:,jj+22),'mcc',0,'ArrowDensity',[30 30]);

Do you have any tips? Please let me know.

Author: Xu Cui Categories: brain, matlab, nirs, programming Tags:

Projection of a NIRS channel on a brain surface

April 11th, 2016

When analyzing the data in our concurrent NIRS-fMRI study, we are particularly interested in how the NIRS signals were correlated to the fMRI signal. To answer the question we need to create an ROI (region of interest) in brain which is directly underneath the NIRS channel (which is on the skull). So a projection of the NIRS channel on the brain surface is necessary.

It might be easy if the brain were a perfect sphere, or at least doable if it is smooth. But brain surface is anything but smooth. What we did is:

  1. Create a brain surface mask. This is easy using BET (http://fsl.fmrib.ox.ac.uk/fsl/fslwiki/BET)
  2. Loop over all voxels on the surface, calculate the distance between the NIRS channel (point A) and the voxel
  3. Find the voxel which is closest to the NIRS channel (It is point B in the figure below)

Projection of NIRS channel on brain surface

Projection of NIRS channel on brain surface

Author: Xu Cui Categories: brain, matlab, nirs, programming Tags:

T value of a single subject

March 14th, 2016

I usually report a group-level T-test image in the final publications or presentations. But in the early stage of a project when only one or a few subjects have been scanned, I often need to report an activation map for an individual subject after general linear model analysis (GLM). I could show a beta image or activation image, but a T-image is sometimes desired. The question is how do I calculate the T-value of a contrast for a single subject.

I will use the following example to show how it is calculated. Assume we have two conditions, beep and flash, and the brain signal is saved in y. Based on the timing of beep and flash, we have an independent variable called x with two columns. The 1st column is for beep, the 2nd for flash.

x=rand(10,2);
y=x(:,1)*2 + x(:,2)*3 + 0.2 + rand(size(x(:,1)))/5;
[a,b,c] = glmfit(x,y);

Now you can easily find the beta value, the T value, and the p-value of each individual condition:

beta = a(2:end)'
T = c.t(2:end)'
p = c.p(2:end)'

What about the contrast between beep and flash? Well, it’s a bit more complicated. First, we need to find the covariance matrix of the betas.

covb = c.covb(2:end, 2:end);

Then we define our contrast vector, in this case it is simply contrast = [1 -1]. The variance of contrast (beep-flash) is

V = contrast  * covb * contrast';

The T-value of the contrast is

T = (beta * contrast') / sqrt(V) * sqrt(length(x));

Note: you might be tempted to use the two-sample T formula to calculate the T-value (i.e. finding the variance of beta1, variance of beta2, and then calculate). It is not the correct way. We have to consider the case where the two conditions are correlated. Using the covariance matrix above is the right way.

Author: Xu Cui Categories: brain, matlab, nirs Tags:

Finger tapping task MatLab script

February 29th, 2016

finger tapping

finger tapping

Finger tapping is probably the most used task in brain imaging studies (fMRI and NIRS). The task is simple and elicits robust brain signal in the motor cortex. We always use it to test new devices or develop a new method.

If you want to control your NIRS device (say Hitachi ETG 4000) to start/end, please refer http://www.alivelearn.net/?p=1260

If you want to see how blood flow increases in brain motor cortex during finger tapping, check out
http://www.alivelearn.net/?p=1686

To download, please fill the form below:






Author: Xu Cui Categories: brain, matlab, nirs, programming Tags:

xjView 8.13 released. Allow to change the minimum of colorbar range

November 9th, 2015

xjView 8.14 just released. Download link: http://www.alivelearn.net/xjview

xjView 8.13 is released with the following new features:

  1. Allow to change the minimum value of the color bar range. This will enable you to do
    1. Create a symmetric color bar from negative to positive. For example, from -8 to 8
    2. Use cold color to show “negative” contrast. In this case, you simply set the max as 0, and min as a positive number such as 5
  2. Allow to resize the xjView window

A short tutorial video: https://youtu.be/jixQvi0ccHc
Download link: http://www.alivelearn.net/xjview

Author: Xu Cui Categories: matlab, programming Tags: