Writing Fast Matlab Code #9:Numerical Integration
Quadrature formulas are numerical approximation of integrals of the form
where the xk are called the nodes or abscissas and wk are the associated weights. Simpson’s rule is
Simpson’s rule is a quadrature formula with nodes a, a + h, b and node weights h/3,4/3h, h/3 .
Matlab has several functions for quadrature in one dimension:
| quad | adaptive Simpson | Better for low accuracy and nonsmooth integrands |
| quadl | adaptive Gauss-Lobatto | Higher accuracy with smoother integrands |
| quadgk | adaptive Gauss-Kronrod | Oscillatory integrands and high accuracy |
| quadv | adaptive Simpson | Vectorized for multi-valued integrands |
The quad- functions are robust and precise, however, they are not very efficient. They use an adaptive refinement procedure, and while this reduces the number of function calls, they gain little from vectorization and incur significant overhead.
If an application requires approximating an integral whose integrand can be efficiently vectorized, using nonadaptive quadrature may improve speed.
for n = −20:20 % Compute Fourier series coefficients c(n + 21) = quad(@(x)(exp(sin(x).ˆ6).*exp(−1i*x*n)),0,pi,1e−4); end
This code runs in 5.16 seconds. In place of quad, using Simpson’s composite rule with N = 199 nodes yields results with comparable accuracy and allows for vectorized computation. Since the integrals are all over the same interval, the nodes and weights only need to be constructed once.
N = 199; h = pi/(N−1); x = (0:h:pi).'; % Nodes w = ones(1,N); w(2:2:N−1) = 4; w(3:2:N−2) = 2; w = w*h/3; % Weights for n = −20:20 c(n + 21) = w * ( exp(sin(x).ˆ6).*exp(−1i*x*n) ); end
This version of the code runs in 0.02 seconds (200 times faster). The quadrature is performed by the dot product multiplication with w. It can be further optimized by replacing the for loop with one vector-matrix multiply:
[n,x] = meshgrid(−20:20, 0:h:pi); c = w * ( exp(sin(x).ˆ6).*exp(−1i*x.*n) );
For this example, quadv can be used on a multi-valued integrand with similar accuracy and speed,
n = −20:20;
c = quadv(@(x)exp(sin(x).ˆ6).*exp(−1i.*x.*n),0,pi,1e−4);
9.1 One-Dimensional Integration
|
|
is approximated by composite Simpson’s rule with
h = (b − a)/(N−1); x = (a:h:b).'; w = ones(1,N); w(2:2:N−1) = 4; w(3:2:N−2) = 2; w = w*h/3; I = w * f(x); % Approximately evaluate the integral
where N is an odd integer specifying the number of nodes.
A good higher-order choice is 4th-order Gauss-Lobatto [4] (as used by quadl), based on
N = max(3*round((N−1)/3),3) + 1; % Adjust N to the closest valid choice h = (b − a)/(N−1); d = (3/sqrt(5) − 1)*h/2; x = (a:h:b).'; x(2:3:N−2) = x(2:3:N−2) − d; x(3:3:N−1) = x(3:3:N−1) + d; w = ones(1,N); w(4:3:N−3) = 2; w([2:3:N−2,3:3:N−1]) = 5; w = w*h/4; I = w * f(x); % Approximately evaluate the integral
The number of nodes N must be such that (N − 1)/3 is an integer. If not, the first line adjusts N to the closest valid choice. It is usually more accurate than Simpson’s rule when f has six continuous derivatives, f ∈ C 6 (a, b).
A disadvantage of this nonadaptive approach is that the accuracy of the result is only indirectly con- trolled by the parameter N. To guarantee a desired accuracy, either use a generously large value for N or, if possible, determine the error bounds 
[6,7].
where h = (b-a)/(N-1). Note that these bounds are valid only when the integrand is sufficiently differentiable: f N −1 , must have four continuous derivatives for the Simpson’s rule error bound, and six continuous derivatives for 4th-order Gauss-Lobatto.
Composite Simpson’s rule is a fast and effective default choice. But depending on the situation, other methods may be faster and more accurate:
- If the integrand is expensive to evaluate, or if high accuracy is required and the error bounds above are too difficult to compute, use one of Matlab’s adaptive methods. Otherwise, consider composite methods.
- Use higher-order methods (like Gauss-Lobatto/quadl) for very smooth integrands and lower-order methods for less smooth integrands.
- Use the substitution u = 1/(1-x) or Gauss-Laguerre quadrature for infinite-domain integrals like f ∞.
9.2 Multidimensional Integration
An approach for evaluating double integrals of the form 
%%% Construct Simpson nodes and weights over x %%% h = (b − a)/(Nx−1); x = (a:h:b).'; wx = ones(1,Nx); wx(2:2:Nx−1) = 4; wx(3:2:Nx−2) = 2; wx = w*h/3; %%% Construct Simpson nodes and weights over y %%% h = (d − c)/(Ny−1); y = (c:h:d).'; wy = ones(1,Ny); wy(2:2:Ny−1) = 4; wy(3:2:Ny−2) = 2; wy = w*h/3; %%% Combine for two−dimensional integration %%% [x,y] = meshgrid(x,y); x = x(:); y = y(:); w = wy.'*wx; w = w(:).'; I = w * f(x,y); % Approximately evaluate the integral
Similarly for three-dimensional integrals, the weights are combined with
[x,y,z] = meshgrid(x,y,z); x = x(:); y = y(:); z = z(:); w = wy.'*wx; w = w(:)*wz; w = w(:).';
|
|
When the integration region is complicated or of high dimension, Monte Carlo integration techniques are appropriate. The disadvantage is that an N -point Monte Carlo quadrature has error on the order O(1/sqrt(n)), so many points are necessary even for moderate accuracy. Suppose that N points,x1 , x2 , . . . , xN , are uniformly randomly selected in a multidimensional region (or volume) Ω. Then
To integrate a complicated region W that is difficult to sample uniformly, find an easier region Ω that contains W and can be sampled [5]. Then
χW (x) is the indicator function of W : χW (x) = 1 when x is within region W and χW (x) = 0 otherwise. Multiplying the integrand by χW sets contributions from outside of W to zero.
%%% Uniformly randomly sample points (x,y) in Ω %%% x = 4*rand(N,1) − 2; y = 4*rand(N,1) − 2; %%% Restrict the points to region W %%% i = (cos(2*sqrt(x.ˆ2 + y.ˆ2)).*x <= y & x.ˆ2 + y.ˆ2 <= 4); x = x(i); y = y(i); %%% Approximately evaluate the integrals %%% area = 4*4; % The area of rectangle Ω M = (area/N) * length(x); Mx = (area/N) * sum(x); My = (area/N) * sum(y);
Region W sampled with N = 1500. The center of mass is ≈ (0.5, 0.7).
More generally, if W is a two-dimensional region contained in the rectangle defined by a ≤ x ≤ b and
c ≤ y ≤ d, the following code approximates 
x = a + (b−a)*rand(N,1); y = c + (d−c)*rand(N,1); i = logical(indicatorW(x,y)); x = x(i); y = y(i); area = (b−a)*(d−c); I = (area/N) * sum(f(x,y)); % Approximately evaluate the integral
where indicatorW(x,y) is the indicator function χW (x, y) for region W .
For refinements and variations of Monte Carlo integration, see for example [1].
Popularity: 1% [?]
GUI Examples #8: Explore selection determination for a buttongroup
function [] = GUI_8() % Demonstrate how to tell which button in a uibuttongroup is selected. % Similar to GUI_6 except that a uibuttongroup which enforces exclusivity % is used. % % Suggested exercise: Make the editbox change the selected radiobutton. % Be sure to check that user input is valid. % % % Author: Matt Fig % Date: 7/15/2009 S.fh = figure('units','pixels',... 'position',[300 300 250 200],... 'menubar','none',... 'name','GUI_8',... 'numbertitle','off',... 'resize','off'); S.bg = uibuttongroup('units','pix',... 'pos',[20 100 210 90]); S.rd(1) = uicontrol(S.bg,... 'style','rad',... 'unit','pix',... 'position',[20 50 70 30],... 'string','Radio 1'); S.rd(2) = uicontrol(S.bg,... 'style','rad',... 'unit','pix',... 'position',[20 10 70 30],... 'string','Radio 2'); S.rd(3) = uicontrol(S.bg,... 'style','rad',... 'unit','pix',... 'position',[120 50 70 30],... 'string','Radio 3'); S.rd(4) = uicontrol(S.bg,... 'style','rad',... 'unit','pix',... 'position',[120 10 70 30],... 'string','Radio 4'); S.ed = uicontrol('style','edit',... 'unit','pix',... 'position',[100 60 50 30],... 'string','1'); S.pb = uicontrol('style','push',... 'unit','pix',... 'position',[75 20 100 30],... 'string','Get Current Radio',... 'callback',{@pb_call,S}); function [] = pb_call(varargin) % Callback for pushbutton. S = varargin{3}; % Get the structure. % Instead of switch, we could use num2str on: % find(get(S.bg,'selectedobject')==S.rd) (or similar) % Note the use of findobj. This is because of a BUG in MATLAB, whereby if % the user selects the same button twice, the selectedobject property will % not work correctly. switch findobj(get(S.bg,'selectedobject')) case S.rd(1) set(S.ed,'string','1') % Set the editbox string. case S.rd(2) set(S.ed,'string','2') case S.rd(3) set(S.ed,'string','3') case S.rd(4) set(S.ed,'string','4') otherwise set(S.ed,'string','None!') % Very unlikely I think. end
Download full list of examples
Popularity: 1% [?]
Bioelectromagnetism Matlab Toolbox: working with EEG/ERP and MRI images
Taken from: http://eeg.sourceforge.net/
License
This toolbox is released under the GNU General Public License (GPL, see http://www.gnu.org/licenses/gpl.html). This is a copyleft license, which means you have the freedom to use, distribute and modify the code, but only on the condition that you must pass on this freedom. You can integrate this code into proprietary packages, but you must do so according to this rule. That is, some parts of your proprietary package will not have this freedom, but those parts derived from this code must retain that freedom. You must use, distribute and develop the code herein in accordance with the GPL.
EEG Features
Firstly, this is not a signal processing toolbox. Of course, once the data is loaded, there are many matlab functions available for data processing, but few of them are integrated into a GUI interface here. At present, there are no specific functions for processing raw EEG, such as filtering, averaging, etc. For examples of signal processing tools, see the matlab signal processing toolbox and the links below, especially EEGLAB.
This toolbox has been developed to facilitate quick and easy import, visualisation and measurement for ERP data. The toolbox can open and visualise ERP averaged data (Neuroscan, ascii formats), 2D/3D electrode coordinates and 3D cerebral tissue tesselations (meshes). All the features can be explored quickly and easily using the example data provided in the toolbox. The GUI interface is simple and intuitive. The following lists the features already available and some items that could be developed.
ERP Visualisation
- ERP data can be read and plotted as a time series
- Automated or GUI entry of ERP epoch/sampling etc. parameters
- Interactive, precise measurement of ERP waveform values
- Interactive ERP peak detection and plotting/measurement
- Interactive ERP topographic mapping
Data Import/Export Support
- Neuroscan EEG formats (.avg,.eeg,.cnt)
- Neuroscan electrode formats (.tri, .3dd ascii)
- EMSE electrode and mesh formats (.elp/.wfr/.reg)
- FreeSurfer mesh formats (.tri/.asc/.surf/.curv/etc)
- BrainStorm formats
- All data is stored internally in one large, convenient data structure (p), which is available from the matlab workspace.
Topographic Mapping
If the electrode position data is available or adapted from the standardized electrode positions available, the toolbox can generate topographic maps. There are various topography options, including 2D/3D surface mapping with various controls for contour mapping, scaling, and colour maps. If a scalp tesselation is available, the toolbox can load and visualise the ‘mesh’ and interpolating from the electrodes onto the mesh (only when they are already coregistered – the functions for coregistration are in early stages of development).
- standardized extended 10/20 electrode coordinates available
- example realistic geometry, with 124 channel electrode coordinates and associated scalp/skull/cortex tissue meshes from MRI volume provided
- latency selection for topographic mapping based on single values, either entered manually or interactively selected
- animation of topographic maps
- automatic or user-defined amplitude scales
- various color or bw topographic maps (linear or polynomial color scales)
- contour topographic mapping, with automatic or user-defined intervals or numbers of contours specified (rudimentary at the moment – needs refinement)
- Printing or saving graphics files (various formats)
- 3D rotation and left,right,front,back views of 3D topographic maps
The following graphic illustrates 3D scalp topography (with interpolation from 124 electrodes onto a scalp mesh). As of May 2002, the methods are integrated with the GUI interface (they are available in the mesh_laplacian.m and mesh_laplacian_interp.m functions). Many thanks to Robert Oostenveld for assistance in validating these functions.
Data Transforms/Analysis
- Identification/replacement of bad electrodes
- ERP peak detection for all electrodes
- ERP peak detection for regions of electrodes
MRI Features
There are useful functions to load and visualize MRI volumes in Analyze format (or the Freesurfer COR- format and GE Signa files). The Analyze avw* functions have been developed to carefully handle the orientation and implement a strict interpretation of the original Analyze 7.5 specification. This specification is available here in two very informative pdf documents:
If you need to, use the orient option in the avw* functions to handle different image orientations, but read the above documents and this discussion on the issue first (you will be wise in no time).
Also, when working with format conversions, consider these enlightening notes from Mark Jenkinson!
It is expected these MRI functions, together with mesh functions, will provide the opportunity to visualize mesh overlays with MRI volumes. It is also creates an avenue for conversion of MRI volumes. There are some MRI processing functions freely available for matlab, some of them are bundled into the CVS archives, but none are integrated into GUI interfaces yet.
For further MRI processing functions, see the matlab image processing toolbox, the SPM toolbox for matlab, and the FSL tools (in c/c++ with source code available).
System Requirements – Development Platform
The development of this matlab toolbox is in its infancy. It is not very clear what the system requirements are, although matlab 6+ is required. I understand from one report that the toolbox GUI does not work under matlab 5.x, but many command line functions should be OK. For most ERP plotting, the toolbox creates about 4-8Mb of data in the workspace and GUI. For more elaborate mesh plotting and interpolation, the toolbox can create up to 40Mb of workspace data (probably that much again in the GUI itself).
The toolbox has been developed on matlab 6.x on a windows platform. I have noticed some minor problems with mesh plotting and interpolation on systems without OpenGL graphics.
Download
Translator
Follow Us
Video Posts
