EEG potential for multiple dipoles at a single EEG sensor
in spherical brain model
eeg = vb_eeg_one_shell(P, Q, R, sigma)
INPUT
Dipole current source
P : current dipole position ( NP x 3 )
Q : current dipole moment ( NP x 3 )
EEG sensor
R : one of sensor position ( 1 x 3 ) or ( 3 x 1 )
sigma : Conductivity inside the sphere
OUTPUT
eeg : EEG Lead field matrix ( NP x 1 )
NP : # of dipole points
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL.46, 245-259, 1999
John C. Mosher, Richard M. Leahy and Paul S. Lewis
2004-12-15 Made by M. Sato
Copyright (C) 2011, ATR All Rights Reserved.
License : New BSD License(see VBMEG_LICENSE.txt)0001 function eeg = vb_eeg_one_shell(P, Q, R, sigma) 0002 % EEG potential for multiple dipoles at a single EEG sensor 0003 % in spherical brain model 0004 % eeg = vb_eeg_one_shell(P, Q, R, sigma) 0005 % INPUT 0006 % Dipole current source 0007 % P : current dipole position ( NP x 3 ) 0008 % Q : current dipole moment ( NP x 3 ) 0009 % EEG sensor 0010 % R : one of sensor position ( 1 x 3 ) or ( 3 x 1 ) 0011 % 0012 % sigma : Conductivity inside the sphere 0013 % 0014 % OUTPUT 0015 % eeg : EEG Lead field matrix ( NP x 1 ) 0016 % 0017 % NP : # of dipole points 0018 % 0019 % IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL.46, 245-259, 1999 0020 % John C. Mosher, Richard M. Leahy and Paul S. Lewis 0021 % 0022 % 2004-12-15 Made by M. Sato 0023 % 0024 % Copyright (C) 2011, ATR All Rights Reserved. 0025 % License : New BSD License(see VBMEG_LICENSE.txt) 0026 0027 % Dipole position vector 0028 P1 = P(:,1); 0029 P2 = P(:,2); 0030 P3 = P(:,3); 0031 0032 % Sensor position vector 0033 R = R(:); 0034 R1 = R(1); 0035 R2 = R(2); 0036 R3 = R(3); 0037 0038 % Difference from dipole to sensor : (NP x 3) 0039 dP = [(R1-P1), (R2-P2), (R3-P3)]; 0040 0041 % Square norm : (NP x 1) 0042 dd = sum(dP.^2, 2); 0043 d = sqrt(dd); 0044 ddd = d.*dd; 0045 pp = sum(P.^2, 2); 0046 0047 rr = sum(R.^2); 0048 r = sqrt(rr); 0049 0050 % Inner product 0051 dp = sum(dP.*P, 2); 0052 dr = dP*R; 0053 pr = P*R; 0054 0055 qr = Q*R; 0056 qp = sum(P.*Q, 2); 0057 0058 % Denominator of Sarvas Eq. 0059 f = dd .* r + d .* dr ; % Eq.(11) 0060 c1 = ( 2*dp./ddd + 1./d - 1./r ); % Eq.(25) 0061 c2 = ( 2./ddd + (d + r)./(r.*f) ); % Eq.(26) 0062 0063 % EEG potential 0064 0065 eeg = (( c1 - c2.*pr )./pp).*qp + c2.*qr; % Eq.(27) 0066 0067 eeg = eeg./(4.0*pi*sigma); 0068 0069 return 0070 0071 %% Execution speed comparison with repmat 0072 % Time 0073 % 0.23 dP = [(R1-P1), (R2-P2), (R3-P3)]; 0074 % 0075 %% This code is faster than the following repmat code for large data 0076 %% since no memory allocation is needed 0077 % ( R1, R2, R3 are scalars) 0078 % 0079 % 0.36 R = repmat(R',[size(P,1) 1]); 0080 % 0.17 dP = R - P;