vb_spherical_harmo_coef

 calculate spherical harmonic coefficient
 [usage]
   A = vb_spherical_harmo_coef( Xhead, Qhead, N, Rmax, Xmeg, Qmeg)
 [input]
   Xhead : head surface (NPx3)
   Qhead : normal vector of head surface (NPx3)
       N : order of function
    Rmax : radius
    Xmeg : sensor coordinates
    Qmeg : sensor orientations
 [output]
       A : calculated coefficient
 [note]

 [history]
   2006.07.06 (Sako) initial version

 Copyright (C) 2011, ATR All Rights Reserved.
 License : New BSD License(see VBMEG_LICENSE.txt)

0001 function A = vb_spherical_harmo_coef(Xhead, Qhead, N, Rmax, Xmeg, Qmeg)
0002 % calculate spherical harmonic coefficient
0003 % [usage]
0004 %   A = vb_spherical_harmo_coef( Xhead, Qhead, N, Rmax, Xmeg, Qmeg)
0005 % [input]
0006 %   Xhead : head surface (NPx3)
0007 %   Qhead : normal vector of head surface (NPx3)
0008 %       N : order of function
0009 %    Rmax : radius
0010 %    Xmeg : sensor coordinates
0011 %    Qmeg : sensor orientations
0012 % [output]
0013 %       A : calculated coefficient
0014 % [note]
0015 %
0016 % [history]
0017 %   2006.07.06 (Sako) initial version
0018 %
0019 % Copyright (C) 2011, ATR All Rights Reserved.
0020 % License : New BSD License(see VBMEG_LICENSE.txt)
0021 
0022 B = vb_multi_sarvas(Xhead, Qhead, Xmeg, Qmeg);
0023 
0024 YG = vb_spherical_grad(Xhead, Qhead, N, Rmax, 'f');
0025 
0026 YY = sqrt(sum(YG .^2));
0027 YY_no_zero = find(YY < eps);
0028 YY(YY_no_zero) = eps;
0029 dist = repmat(YY,size(YG,1),1);
0030 YG_bar = YG ./ dist;
0031 
0032 x1 = diag(1 ./ YY);
0033 x2 = pinv(YG_bar' * YG_bar);
0034 x3 = YG_bar' * B;
0035 A = x1 * (x2 * x3);
0036 
0037 return;