|
From: | Joe Koski |
Subject: | Re: Moving Least Squares, Anyone? |
Date: | Mon, 21 Nov 2005 13:22:38 -0700 |
User-agent: | Microsoft-Entourage/11.2.1.051004 |
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:
But the function has been working for me for quite some time now. And you can certainly use it to see if the algorithm does what you need.
- the function arguments x and y are ordered the wrong way round
- the calculation is slow because it is not vectorized (does anybody know a way to vectorize it?)
- the treatment of the vector edges can probably be improved a lot
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:
On Nov 18, 2005, at 5:27 PM, Dmitri A. Sergatskov wrote:
On 11/18/05, Joe Koski <address@hidden> <mailto:address@hidden> wrote:
Has anyone run across "moving least squares" code in .m, .cc, .f, etc.
formats that would be adaptable for use with octave? Apparently, moving
least squares can be used to create approximating curves for
interpolation/approximation much in the same way that spline curves can be
used. In some situations, the moving least squares approach can reduce the
need to solve large matrices, which is normally associated with the curve
fit process.
A Google of "moving least squares" shows several papers, but no code that I
can find. Any ideas?
Would Kalman filter do the same / better?
Savitsky-Golay is a moving least squares filter. You can use it for
curve smoothing if the data are equally spaced.
- Paul
-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.
Octave's home on the web: http://www.octave.org
How to fund new projects: http://www.octave.org/funding.html
Subscription information: http://www.octave.org/archive.html
-------------------------------------------------------------
[Prev in Thread] | Current Thread | [Next in Thread] |