smooth_jd

 Gaussian smoothering of function xf on grid x with smooting Gaussian
 width dxs (default: mean step size of x)
 
 by J. Decker 30/06/06

0001 function xfs = smooth_jd(x,xf,dxs)
0002 %
0003 % Gaussian smoothering of function xf on grid x with smooting Gaussian
0004 % width dxs (default: mean step size of x)
0005 %
0006 % by J. Decker 30/06/06
0007 %
0008 if nargin < 2,
0009     error('Not enough input arguments')
0010 end
0011 %
0012 Sx = size(x);
0013 if length(Sx) ~= 2 || min(Sx) ~= 1 || max(Sx) == 1,
0014     error('x must be a vector')
0015 else
0016     x = x(:);
0017     nx = length(x);
0018 end
0019 %
0020 if min(diff(x)) <= 0,
0021     error('x must be monotonic');
0022 end
0023 %
0024 if nargin < 3,
0025     dxs = mean(diff(x));
0026 end
0027 %
0028 Sf = size(xf);
0029 if Sf(1) ~= nx,
0030     error('The first dimension of xf must have the same length as x')
0031 end
0032 %
0033 Sdx = size(dxs);
0034 if Sdx(1) ~= 1,
0035     error('The first dimension of dxs must be 1')
0036 end
0037 %
0038 Sff = Sf(2:length(Sf));
0039 %
0040 if ~(length(Sdx) == length(Sf) && all(Sdx(2:end) == Sf(2:end))) 
0041     if ~(length(Sdx) == 2 && Sdx(2) == 1),
0042         error('dxs must be a scalar or of the same 2...N dimensions as xf')
0043     else
0044         dxs = repmat(dxs,[1,Sff]);
0045         Sdx = [1,Sff];
0046     end
0047 end
0048 %
0049 x_S = [x(1) - (x(2) - x(1))/2; (x(1:nx-1) + x(2:nx))/2; x(nx) + (x(nx) - x(nx-1))/2];
0050 dx = diff(x_S);%x grid step size
0051 %
0052 Xx = repmat(x,[1,nx,Sff]);
0053 Xx0 = repmat(x.',[nx,1,Sff]);
0054 Xdx0 = repmat(dx.',[nx,1,Sff]);
0055 Xf0 = repmat(reshape(xf,[1,Sf]),[nx,ones(1,length(Sf))]);
0056 Xdxs = repmat(reshape(dxs,[1,Sdx]),[nx,nx,ones(1,length(Sdx)-1)]);
0057 %
0058 xfs = reshape(sum(Xf0.*Xdx0.*exp(-(Xx - Xx0).^2./Xdxs.^2),2)./...
0059               sum(     Xdx0.*exp(-(Xx - Xx0).^2./Xdxs.^2),2),Sf);
0060 %