eeg_interp_sph_spline_c - Spherical Spline Coefficients
Useage: [C] = eeg_interp_sph_spline_c(Zi,G)
Zi - Nx1 is an EEG/ERP measurement at time t from N sensors
G - NxN Legendre function of the cosine matrix for N sensor positions
eg, COS = elec_cosines(Ei,Ei); where Ei is Nx3 [X Y Z]
G = eeg_interp_sph_spline_g(COS)
C - spherical spline coefficients for Ei (includes co = C(1))
Notes: This function calculates the spherical spline coefficients of
Perrin et al (1989). Electroenceph. & Clin.
Neurophysiology, 72: 184-187. Corrigenda (1990),
Electroenceph. & Clin. Neurophysiology, 76: 565.
(see comments in the .m file for details).
see elec_cosines, eeg_interp_sph_spline, eeg_interp_sph_spline_g0001 function [C] = eeg_interp_sph_spline_c(Zi,G) 0002 0003 % eeg_interp_sph_spline_c - Spherical Spline Coefficients 0004 % 0005 % Useage: [C] = eeg_interp_sph_spline_c(Zi,G) 0006 % 0007 % Zi - Nx1 is an EEG/ERP measurement at time t from N sensors 0008 % G - NxN Legendre function of the cosine matrix for N sensor positions 0009 % eg, COS = elec_cosines(Ei,Ei); where Ei is Nx3 [X Y Z] 0010 % G = eeg_interp_sph_spline_g(COS) 0011 % 0012 % C - spherical spline coefficients for Ei (includes co = C(1)) 0013 % 0014 % Notes: This function calculates the spherical spline coefficients of 0015 % Perrin et al (1989). Electroenceph. & Clin. 0016 % Neurophysiology, 72: 184-187. Corrigenda (1990), 0017 % Electroenceph. & Clin. Neurophysiology, 76: 565. 0018 % (see comments in the .m file for details). 0019 % 0020 % see elec_cosines, eeg_interp_sph_spline, eeg_interp_sph_spline_g 0021 0022 % $Revision: 1.2 $ $Date: 2005/07/13 06:03:50 $ 0023 0024 % Licence: GNU GPL, no implied or express warranties 0025 % History: 08/2001 Darren.Weber_at_radiology.ucsf.edu, with 0026 % mathematical advice from 0027 % Dr. Murk Bottema (Flinders University of SA) 0028 % 10/2003 Darren.Weber_at_radiology.ucsf.edu, with 0029 % mathematical advice and LegendreP function from 0030 % Dr. Tom.Ferree_at_radiology.ucsf.edu 0031 % revised, tested & initial verification complete 0032 % 0033 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0034 0035 eegversion = '$Revision: 1.2 $'; 0036 fprintf('EEG_INTERP_SPH_SPLINE_C [v %s]\n',eegversion(11:15)); tic 0037 0038 %-------------------------------------------------------------------------- 0039 % Check for correct size & orientation of Zi & G 0040 [g1,g2] = size(G); 0041 [z1,z2] = size(Zi); 0042 if g1 ~= g2, error('G matrix must be square'); end 0043 if z1 < z2, Zi = Zi'; [z1,z2] = size(Zi); end 0044 if ~and(isequal(g1,z1),isequal(z2,1)), 0045 error('...G must be NxN & Zi must be Nx1'); 0046 end 0047 nElectrodes = z1; % The number of electrodes 0048 clear g1 g2 z1 z2; 0049 0050 0051 %-------------------------------------------------------------------------- 0052 % Computations described by Perrin et al. (1989) 0053 % 0054 % Let z(i) be the potential value measured at the ith electrode whose 0055 % spherical projection is denoted Ei. The value at any point E of the 0056 % sphere, of the spherical spline U that interpolates from zi at all Ei, 0057 % is given by (Eq. 1): 0058 % 0059 % U(E) = Co + ( for i=1:nElectrodes, sum = (sum + (C(i) * g(cos(E,E(i)))) ) 0060 % 0061 % Where the column array C(i) contains the spherical spline coefficients. 0062 % These are solutions of the simultaneous equations (Eq. 2): 0063 % 0064 % a) GC + TCo = Z 0065 % b) T'C = 0 0066 % 0067 % with: 0068 % T' = ( 1, 1,..., 1), 0069 % C' = (c1,c2,...,cn), - unknown 0070 % Z' = (z1,z2,...,zn), - potentials at Ei electrodes 0071 % G = (Gij) = (g(cos(Ei,Ej))) - cosine matrix betweeen E and Ei 0072 % 0073 % That is, cos(Ei,Ej) denotes the cosine of the angle between the 0074 % electrode projections Ei and the sphere interpolation points Ej. 0075 % 0076 % The function g(x) is the sum of the series (Eq. 3), 0077 % 0078 % g(x) = (1/(4*pi) * ... 0079 % for n=1:inf, sum = sum + ( ( (2*n+1)/(n^m * (n+1)^m) ) * Pn(x) ) 0080 % 0081 % where m is a constant > 1 and Pn(x) is the nth degree Legendre 0082 % polynomial. Perrin et al. (1989) evaluated m=1:6 and recommend m=4, 0083 % for which the first 7 terms of Pn(x) are sufficient to obtain 0084 % a precision of 10^-6 for g(x) (ie, the above for loop is n=1:7). 0085 % Perrin et al. (1989) recommend tabulation of g(x) for 0086 % x = linspace(-1,1,2000) to be used as a lookup given actual 0087 % values for cos(E(i),E(j)). 0088 0089 % Initialise T and C 0090 T = ones(nElectrodes,1); 0091 C = zeros(nElectrodes,1); % so Eq. 2b, T' * C = 0 0092 0093 % G is NxN, where N = nElectrodes 0094 % C is Nx1 coefficients that will solve Eq. 2a 0095 % Z is Nx1 potential values 0096 0097 % GC is Nx1, a linear combination of the columns of G 0098 % TCo is Nx1, a constant column vector 0099 % 0100 % We can rearrange Eq. 2a so we have 0101 % GC = Z, effectively incorporating TCo into G,C and Z. 0102 % First, we consider Co at the top of C, so we get C' = (co,c1,c2,...cn) 0103 % and C becomes an (nElectrodes+1)x1 column array (ie [0;C]). (Note that 0104 % we solve for C below, so this concatenation is not actually done.) 0105 % Then we have to concatenate [0 T'] onto the top of G (ie, [ [0 T']; G]), 0106 % so that [0 T'] * C = 0 (this is the first line of Gx multiplied by C). 0107 % Lastly, we have to concatenate a column of ones into the first column 0108 % of G and place a zero into G(1,1), so that effectively we add co to the 0109 % sum of multiplying G by C. The resulting linear system is: 0110 % 0111 % [ 0 1 .. 1] [co] [ 0] 0112 % [ 1 ] [c1] [z1] 0113 % [ . ] [c2] = [z2] 0114 % [ . G ] [..] [..] 0115 % [ 1 ] [cn] [zn] 0116 % 0117 % so... 0118 0119 % add ones to first row of G, call it Gx 0120 Gx = [T'; G]; 0121 % now add ones to first column of Gx, but 0122 % with a zero at Gx(1,1) 0123 Gx = [[0;T] Gx]; 0124 0125 % Append 0 to the top of Zi (Z becomes (nElectrodes+1)x1 ) 0126 Z = [0;Zi]; 0127 0128 % Find the inverse of Gx 0129 %invGx = inv(Gx); 0130 0131 % or 0132 invGx = Gx\eye(size(Gx)); 0133 0134 % or 0135 %lambda = 0; 0136 %tolerance = 10^-7 * 10^(lambda*-1); 0137 %invGx = pinv(Gx,tolerance); 0138 0139 % Calculate inverse of Gx 0140 % Ginv = inv(GT); 0141 % EMSE NOTE: The inversion of G is sensitive to errors because 0142 % electrode locations are not really on a sphere. Therefore, 0143 % the solution should be regularized. In EMSE v4.2, the 0144 % regularization involves Singluar Value Decomposition: 0145 % G = UWV' 0146 % inv(G) approx = V W- U' 0147 % 0148 % W- = [w-(i,j)], w-(i,j) { 1/w(ij) for w(i,j) > cutoff 0149 % 0 otherwise 0150 % 0151 % cutoff is determined by a smoothing parameter: 0152 % 0153 % cutoff = k * 10 ^ -lambda 0154 % 0155 % where k = { 10^-5 for laplacian 0156 % 10^-7 for voltage 0157 % and lambda is a the smoothing parameter. 0158 % Smaller values of lambda (<0), fewer terms are 0159 % truncated and the solution is less smooth. 0160 % 0161 % [U,W,S] = svd(Gx); %Gx = U*W*S'; 0162 % invW = inv(W); 0163 % lambda = 0; 0164 % cutoff = 10^-7 * 10^(lambda*-1); 0165 % Winv(find(invW<cutoff)) = 0; 0166 % invGx = S * invW * U.'; 0167 0168 0169 0170 % now we solve for C = inv(Gx)*Z 0171 C = invGx * Z; 0172 0173 Co = C(1); 0174 Ci = C(2:end); 0175 0176 t=toc; fprintf('...done (%6.2f sec)\n\n',t); 0177 0178 return