Category ArchiveTips
Tips Martin Pyka on 03 May 2012
Changing the starting position of the EM for DCM
When you have already preassumptions with regard to the coupling strength of the neural parameters or the hemodynamic parameters you can start the estimation of a DCM with a certain starting position.
All you have to do is, create a variable DCM.M.P and assign new values to it, e.g:
load dcm_file
DCM.M.P = DCM.M.pE;
% alter DCM.M.P
save dcm_file DCM
The EM-algorithm will begin with the parameters specified in DCM.M.P.
Tips Martin Pyka on 01 May 2012
Understanding the framework of DCM
DCM is implemented in a framework that can be used to develop and test other models using the Expectation Maximization algorithm and Bayesian techniques. To understand how to use other models within this framework, I can recommend to look at spm_nlsi.m. After the return-command, at the end of the script, there is a nice short example, how to setup everything in order to test your own model.
And by the way: the code for the dynamic causal model can be found in spm_fx_dcm.m (the neural signal is in y(:,1), some hemodynamic parameters are stored in y(:,2:5)) and spm_gx_dcm.m (the hemodynamic forward model).
Tips Martin Pyka on 23 Mar 2012
PPI
An interesting read about Psychophysiological Interactions and some issues related to deconvolution and orthogonality of the regressors can be found here at the FMRIB Centre.
Tips Martin Pyka on 16 Sep 2011
Screen
Gerrit from the blog Sustainable Research posted a nice article about the screen-command for the linux console that allows you to run scripts in the background without an open console window. This can be also useful, when you want to run SPM-scripts without a user interface and remote on a server. You just start the script within a screen, exit the screen and close the remote connection. And with Putty for Android or iPhone you can supervise your scripts from everywhere :).
Tips Martin Pyka on 31 Aug 2011
Display and save brain images
I suppose that this is not the most efficient way to do this, but it works. With the following commands you can easily write a script that displays any kind of nifti-image (raw files, mask, beta-maps) across all subjects and saves it as a JPG image. This can be handy for quality assurance:
spm_image('init','niftifile.nii') % display a nifti-file
spm_orthviews('Reposition', [22 0 -21]) % jump to a certain coordinate
fs = get(0,’Children’);
res = getframe(fs(1));
imwrite(res.cdata, ‘niftifile.jpg’); % save the window as JPG
This should work for SPM5 and SPM8. Any suggestions how to improve the “save as jpg” procedure?
Tips Martin Pyka on 26 Aug 2011
Execute SPM scripts without popups
I guess I am not the only one who uses matlab scripts to analyze large amounts of fMRI data. To get rid of the popups of SPM during the processing, you can start Matlab either with
matlab -nodesktop
to solely start matlab on the command line, or you start it with
matlab -noFigureWindows
which starts Matlab as normal but prevents the depiction of any additional windows.
This post was inspired by a question that I received by email. Thanks!
Tips Martin Pyka on 25 Aug 2011
Manual high-pass filtering
The high-pass filter can have a great effect on the subsequent parameter estimation. To explore, how a certain setting (default is 128) influences the time course of the regressor you can use the spm_filter function. First, we build a HRF-regressor:
reg = [zeros(10,1); ones(10,1); zeros(30,1); ones(10,1); zeros(30,1); ones(10,1); zeros(30,1)];
bf = spm_get_bf %TR is 2.2 in our case
U.u = reg;
U.name = {'reg'};
% create the regressor convolved with the HRF
convreg = spm_Volterra(U, bf.bf);
In your SPM.mat you can find this kind of regressors under SPM.xX.X.
spm_filter is a function that either creates the low frequencies to be removed or removes those frequencies from the time course. The input variable K defines all necessary parameters:
K.RT = 2.2;
K.row = 1:130; % 130 corresponds to the length of convreg
K.HParam = 128; % cut-off period in seconds
Here, I defined K.row as it also happens internally. But I guess, ones(1,130) should also work, as it seems that K.row just serves as a mask. Now, we compute the low frequencies that should be removed from convreg.
nK = spm_filter(K);
nK.X0 contains the frequencies that can be removed. A second call of spm_filter removes these frequencies.
Y = spm_filter(nK, convreg);
In SPM.mat these regressors are stored in SPM.xX.nKX and are depicted in the design matrix. As you can see, the high-pass filter cut-off period can greatly affect the regressors that are used for the parameter estimation:
K.HParam = 256;
figure; plot(spm_filter(spm_filter(K), convreg));
K.HParam = 64;
figure; plot(spm_filter(spm_filter(K), convreg));
Tips Martin Pyka on 23 Aug 2011
spm_dcm_create
If you experiment with DCM, it might be interesting to know that you can build synthetic input data for DCM with the spm_dcm_create-function. Just type spm_dcm_create on the console, press enter and follow the steps. Other interesting units for DCM are: spm_dcm_estimate (the estimation procedure) and spm_dcm_compare (SPM8) for comparing two or more models.
Tips Martin Pyka on 19 Aug 2011
Manual convolution with the HRF
With the help of two SPM-functions you can manually convolve time series in Matlab. This might help to understand how the convolution with the HRF (or any other function) affects the time series. Here is an example:
% create a time series
reg = [zeros(10,1); ones(10,1); zeros(10,1); ones(10,1); zeros(10,1)];
plot(reg)
% create a basis function, in this case the hrf
bf = spm_get_bf
% with a tr of 2.2 and HRF as basis function you receive
dt: 2.2000
name: 'hrf'
length: 33
order: 1
bf: [15x1 double]
% these are necessary definitions for the spm_Volterra-function
U.u = reg;
U.name = {’reg’}; % U needs a name, but the string has no meaning
% convolve reg with the hrf
convreg = spm_Volterra(U, bf.bf);
plot(convreg)
Tips Martin Pyka on 31 May 2010
Convert matrix coordinates to image coordinates and vice versa
To convert matrix coordinates for a given Nifti/Analyze-image into the image coordinates and vice versa, you need to load the transformation matrix of the image into the workspace:
vol = spm_vol('image.img')
vol.mat % displays the transformation matrix
% trasform from matrix/voxel coordinates to image coordinates (the fourth element has to be 1)
vol.mat * [x y z 1]‘
% transform image coordinates (e.g. mni space) to matrix coordinates
inv(vol.mat) * [x y z 1]‘