Asymmetric least squares smoothing

// by tony.withers@uwo.ca, using method of Eilers, PHC and Boelens, HFM
// (2005) Baseline correction with asymmetric least squares smoothing.

// Creates (and overwrites) w_base, a baseline estimate for w_data. The
// asymmetry parameter (Eilers and Boelens' p) generally takes values
// between 0.001 and 0.1. Try varying lambda in orders of magnitude
// between 10^2 and 10^9. Not efficient for large N, try it for w_data
// with fewer than 1000 points.
function ALS(w_data, lambda, asymmetry)
    wave w_data
    variable lambda, asymmetry
   
    variable i, N=numpnts(w_data), rms=inf
    variable maxIts=20
   
    matrixOp /free  D = identity(N)
    differentiate /EP=1/METH=2/DIM=0 D
    differentiate /EP=1/METH=2/DIM=0 D

    // this step (specifically the matrix multiplication) is slow:
    matrixOp /free H = lambda * (D^t x D)

    duplicate /o/free w_data w, w_new
    w=1

    for (i=0;i<maxIts;i+=1)
        matrixOp /o/free  C = chol(diagRC(w, N, N)+H)
        matrixOp /o w_base = backwardSub(C,(forwardSub(C^t, w * w_data)))
        w_new = asymmetry * (w_data>w_base) + (1-asymmetry) * (w_data<w_base)
       
        // convergence test
        w-=w_new
        wavestats /Q w
        if (v_rms>=rms)    
            return i+1
        else
            rms=v_rms
            w=w_new
        endif
    endfor
    return 0
end

Forum

Support

Gallery

Igor Pro 8

Learn More

Igor XOP Toolkit

Learn More

Igor NIDAQ Tools MX

Learn More