One strange way to generate b-spline base

Question

I am trying to use curve fitting under Matlab. There are two kinds of spline in Matlab: piecewise polynomials and b-spline. For b-spline, we know that the basic functions can be derived by means of a recursion formula. You can see the below link: http://en.wikipedia.org/wiki/B-spline

From matlab file exchange, one provide one function called SPLINEFIT which can perform curve fitting based on B-splines. I have written one program to use the function and try to read the implementation. I can't understand the way that basis functions are produced. The author provides one way to produce basis function and even the author say in his reply to others' comments that I admit that the splinebase function is a bit cryptic (also for the author). You can find the function from the link:http://www.mathworks.com/matlabcentral/fileexchange/13812-splinefit.

The function that produce the basis functions is followed.

    %--------------------------------------------------------------------------
    function pp = splinebase(breaks,n)
    %SPLINEBASE Generate B-spline base PP of order N for breaks BREAKS

    breaks = breaks(:);     % Breaks
    breaks0 = breaks';      % Initial breaks
    h = diff(breaks);       % Spacing
    pieces = numel(h);      % Number of pieces
    deg = n - 1;            % Polynomial degree

    % Extend breaks periodically
    if deg > 0
        if deg <= pieces
            hcopy = h;
        else
            hcopy = repmat(h,ceil(deg/pieces),1);
        end
        % to the left
        hl = hcopy(end:-1:end-deg+1);
        bl = breaks(1) - cumsum(hl);
        % and to the right
        hr = hcopy(1:deg);
        br = breaks(end) + cumsum(hr);
        % Add breaks
        breaks = [bl(deg:-1:1); breaks; br];
        h = diff(breaks);
        pieces = numel(h);
    end

    % Initiate polynomial coefficients
    coefs = zeros(n*pieces,n);
    coefs(1:n:end,1) = 1;

    % Expand h
    ii = [1:pieces; ones(deg,pieces)];
    ii = cumsum(ii,1);
    ii = min(ii,pieces);
    H = h(ii(:));

    % Recursive generation of B-splines
    for k = 2:n
        % Antiderivatives of splines
        for j = 1:k-1
            coefs(:,j) = coefs(:,j).*H/(k-j);
        end
        Q = sum(coefs,2);
        Q = reshape(Q,n,pieces);
        Q = cumsum(Q,1);
        c0 = [zeros(1,pieces); Q(1:deg,:)];
        coefs(:,k) = c0(:);
        % Normalize antiderivatives by max value
        fmax = repmat(Q(n,:),n,1);
        fmax = fmax(:);
        for j = 1:k
            coefs(:,j) = coefs(:,j)./fmax;
        end
        % Diff of adjacent antiderivatives
        coefs(1:end-deg,1:k) = coefs(1:end-deg,1:k) - coefs(n:end,1:k);
        coefs(1:n:end,k) = 0;
    end

    % Scale coefficients
    scale = ones(size(H));
    for k = 1:n-1
        scale = scale./H;
        coefs(:,n-k) = scale.*coefs(:,n-k);
    end

    % Reduce number of pieces
    pieces = pieces - 2*deg;

    % Sort coefficients by interval number
    ii = [n*(1:pieces); deg*ones(deg,pieces)];
    i = cumsum(ii,1);
    coefs = coefs(ii(:),:);

    % Make piecewise polynomial
    pp = mkpp(breaks0,coefs,n);

Answer 1

I'm not Matlab literate, so I can't read your code. But I know that much of the Matlab spline code was written by Carl deBoor. It's just the code from his book "A Practical Guide to Splines" converted into Matlab language. Reading the book will probably help.

After a bit more reading, it looks like the basic b-spline recursion is in these two statements

   coefs(:,j) = coefs(:,j).*H/(k-j);
   coefs(1:end-deg,1:k) = coefs(1:end-deg,1:k) - coefs(n:end,1:k);