help-octave
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: Moving polynomial fit


From: Ben Abbott
Subject: Re: Moving polynomial fit
Date: Sat, 05 Sep 2009 13:23:11 -0400

On Sep 5, 2009, at 11:09 AM, babelproofreader wrote:

I would like to write a script/function for use on a time series that will least squares fit a polynomial of a given degree to a moving window across the time series in the same way as a moving average computes an average of a
moving window, the output of the script/function being a vector of the
extreme right hand edge/last values of the fitted polynomial over the moving
window. I have tried using the Savitsky-Golay filters from the Signals
package but they are not suitable for my purposes as the end conditions mean that the values of the fitted polynomial at the right hand edges of the moving window are different from the values calculated when the data are no
longer at the right hand edge. What in-built Octave functions should I
consider using to achieve the above i.e. polyfit or some other function(s)? Also, as loops can be quite slow, is there a more efficient way of coding
the above rather than having to resort to a for loop?

There might be a better way ... Assuming you'd like to center your window about the data point of interest ...

---------------- begin: mwpolysmooth.m ---------------
function ym = mwpolysmooth (x, y, order, window)
% USAGE: ym = mwpolysmooth (x, y, order, window)

  xc = num2cell(x);
  idx = @(xc) find (abs (x - xc) < 0.5*window);
  n = cellfun (idx, xc, 'uniformoutput', false);
  mpfit = @(n) polyfit (x(n), y(n), order);
  pm = cellfun (mpfit, n, 'uniformoutput', false);
  nc = num2cell (1:numel(x));
  mpval = @(n) polyval (pm{n}, x(n));
  ym = cell2mat (cellfun (mpval, nc, 'uniformoutput', false));

end

%!demo
%! x = 0:0.1:10;
%! y = rand (size(x)) + cos (x);
%! window = 5;
%! order = 2;
%! ym = mwpolysmooth (x, y, order, window);
%! clf
%! plot (x, y, 's', x, ym)
---------------- end: mwpolysmooth.m ---------------

Ben






reply via email to

[Prev in Thread] Current Thread [Next in Thread]