|
From: | Thorsten Meyer |
Subject: | Re: Moving Least Squares, Anyone? |
Date: | Sat, 19 Nov 2005 14:15:58 +0100 |
User-agent: | Mozilla/5.0 (X11; U; Linux i686; en-US; rv:1.7.12) Gecko/20051007 Debian/1.7.12-1 |
I have been using the attached function for data smoothing and
calculation of smoothed derivatives for non-equally spaced data. It
uses moving least-squares polynomial fits of the given data. If I
remember it right, the savitzky-golay algorithm does the same, only
that the procedure is nicely vectorized for equally spaced data. My code is terrible:
regards Thorsten function dydx = deriv(x1, y1, nl, nr, m, ld, edge) % deriv - calculates numerically smoothed derivatives % % calling syntax: % % dydx = deriv(x, y, nl, nr, m, ld, edge); % % where x is a vector of x-values % y is a vector of y-values (y(i) = f(x(i))) % nl is the number of leftward data points % used for smoothing % nr is the number of rightward data points % m is the order of the polynomials fitted % to the data % ld is the order of derivative to be calculated % (default is 1). For ld=0, the function % itself is smoothed % edge is a flag that determines the way the % routine deals with the edges of the % vector y. % edge = 1 : the last full polynomial % based on the values y(1:nl+nr+1) % and y(length(y)-nl-nr:length(y)) % respectively is used to calculate % all the values dydx(1:nl) and % dydx(length(y)-nl-nr-1:length(y)) % respectively. % edge = 2 : the window used for fitting is % linearly shrunk and shifted as % the edges are approached % (default is 2) % dydx is assigned the n-th derivative of y(x) % if nargin < 7 edge = 2; if nargin < 6 ld = 1; end; end; if m < ld m = ld; end; faktor = gamma(ld + 1); xmax = max(abs(x1)); ymax = max(abs(y1)); x = x1 / xmax; y = y1 / ymax; ly = length(x); dydx = zeros(size(x)); for i = 1+nl : ly-nr xx = x(i-nl:i+nr) - x(i); yy = y(i-nl:i+nr); pp = polyfit(xx, yy, m); dydx(i) = faktor * pp(m+1-ld); end; if edge == 1 xx = x(1:(nl + nr + 1)) - x(nl+1); yy = y(1:(nl + nr + 1)); pp = polyfit(xx, yy, m); dpp = polydiff(pp, ld); dydx(1:nl) = polyval(dpp, xx(1:nl)); xx = x((ly - (nl + nr)) : ly) - x(ly - nr); yy = y((ly - (nl + nr)) : ly); pp = polyfit(xx, yy, m); dpp = polydiff(pp, ld); dydx((ly - nr + 1) : ly) = polyval(dpp, xx(nl+2:nl+1+nr)); elseif edge == 2 nn = nl + nr + 1; if nl < 2 nnn = m * 2; else nnn = round(linspace(m * 2, nl + nr, nl)); end for i = 1:nl xx = x(1:nnn(i)) - x(i); yy = y(1:nnn(i)); pp = polyfit(xx, yy, m); dydx(i) = faktor * pp(m+1-ld); end; if nr < 2 nnn = m * 2; else nnn = round(linspace(nl + nr, m * 2, nr)); end il = ly - nnn + 1; for i = 1:nr xx = x(il(i):ly) - x(ly - nr + i); yy = y(il(i):ly); pp = polyfit(xx, yy, m); dydx(ly - nr + i) = faktor * pp(m+1-ld); end; end; if ld == 0 dydx = dydx * ymax; else dydx = dydx * ymax / xmax ^ ld; end Paul Kienzle wrote:
|
[Prev in Thread] | Current Thread | [Next in Thread] |