module Savgol
Constants
- VERSION
Public Class Methods
savgol(array, window_points=11, order=4, deriv: 0, check_args: false)
click to toggle source
# File lib/savgol.rb, line 105 def savgol(array, window_points=11, order=4, deriv: 0, check_args: false) sg_check_arguments(window_points, order) if check_args half_window = (window_points -1) / 2 weights = sg_weights(half_window, order, deriv) padded_array = sg_pad_ends(array, half_window) sg_convolve(padded_array, weights) end
savgol_uneven(xvals, yvals, window_points=11, order=4, check_args: false, new_xvals: nil)
click to toggle source
Does simple least squares to fit a polynomial based on the given x values. The major speed-boost of doing savgol is lost for uneven time points. TODO: implement for different derivatives; implement with a window that is of fixed size (not based on the number of points)
# File lib/savgol.rb, line 24 def savgol_uneven(xvals, yvals, window_points=11, order=4, check_args: false, new_xvals: nil) sg_check_arguments(window_points, order) if check_args half_window = (window_points -1) / 2 yvals_padded = sg_pad_ends(yvals, half_window) xvals_padded = sg_pad_xvals(xvals, half_window) xvals_size = xvals.size if new_xvals nearest_xval_indices = sg_nearest_index(xvals, new_xvals) new_xvals.zip(nearest_xval_indices).map do |new_xval, index| sg_regress_and_find( xvals_padded[index,window_points], yvals_padded[index,window_points], order, new_xval ) end else xs_iter = xvals_padded.each_cons(window_points) yvals_padded.each_cons(window_points).map do |ys| xs = xs_iter.next sg_regress_and_find(xs, ys, order, xs[half_window]) end end end
sg_check_arguments(window_points, order)
click to toggle source
# File lib/savgol.rb, line 113 def sg_check_arguments(window_points, order) if !window_points.is_a?(Integer) || window_points.abs != window_points || window_points % 2 != 1 || window_points < 1 raise ArgumentError, "window_points size must be a positive odd integer" end if !order.is_a?(Integer) || order < 0 raise ArgumentError, "order must be an integer >= 0" end if window_points < order + 2 raise ArgumentError, "window_points is too small for the polynomials order" end end
sg_convolve(data, weights, mode=:valid)
click to toggle source
# File lib/savgol.rb, line 125 def sg_convolve(data, weights, mode=:valid) data.each_cons(weights.size).map do |ar| ar.zip(weights).map {|pair| pair[0] * pair[1] }.reduce(:+) end end
sg_nearest_index(original_vals, new_vals)
click to toggle source
returns the nearest index in original_vals for each value in new_vals (assumes both are sorted arrays). complexity: O(n + m)
# File lib/savgol.rb, line 53 def sg_nearest_index(original_vals, new_vals) newval_iter = NilEnumerator.new(new_vals.each) indices = [] index_iter = NilEnumerator.new((0...original_vals.size).each) index = index_iter.next newval=newval_iter.next last_index = original_vals.size-1 until newval >= original_vals[index] indices << index break unless newval = newval_iter.next end loop do break unless newval if index.nil? indices << last_index else until newval <= original_vals[index] index = index_iter.next if !index indices << last_index break end end if index if newval < original_vals[index] indices << ((newval - original_vals[index-1]) <= (original_vals[index] - newval) ? index-1 : index) else indices << index end end end newval = newval_iter.next end indices end
sg_pad_ends(array, half_window)
click to toggle source
pads the ends with the reverse, geometric inverse sequence
# File lib/savgol.rb, line 132 def sg_pad_ends(array, half_window) start = array[1..half_window] start.reverse! start.map! {|v| array[0] - (v - array[0]).abs } fin = array[(-half_window-1)...-1] fin.reverse! fin.map! {|v| array[-1] + (v - array[-1]).abs } start.push(*array, *fin) end
sg_pad_xvals(array, half_window)
click to toggle source
pads the ends of x vals
# File lib/savgol.rb, line 144 def sg_pad_xvals(array, half_window) deltas = array[0..half_window].each_cons(2).map {|a,b| b-a } start = array[0] prevals = deltas.map do |delta| newval = start - delta start = newval newval end prevals.reverse! deltas = array[(-half_window-1)..-1].each_cons(2).map {|a,b| b-a } start = array[-1] postvals = deltas.reverse.map do |delta| newval = start + delta start = newval newval end prevals.push(*array, *postvals) end
sg_regress_and_find(xdata, ydata, order, xval)
click to toggle source
# File lib/savgol.rb, line 95 def sg_regress_and_find(xdata, ydata, order, xval) xdata_matrix = xdata.map { |xi| (0..order).map { |pow| (xi**pow).to_f } } mx = Matrix[*xdata_matrix] my = Matrix.column_vector(ydata) ((mx.t * mx).inv * mx.t * my).transpose.to_a[0]. zip((0..order).to_a). map {|coeff, pow| coeff * (xval**pow) }. reduce(:+) end
sg_weights(half_window, order, deriv=0)
click to toggle source
returns an object that will convolve with the padded array
# File lib/savgol.rb, line 166 def sg_weights(half_window, order, deriv=0) # byebug mat = Matrix[ *(-half_window..half_window).map {|k| (0..order).map {|i| k**i }} ] # Moore-Penrose psuedo-inverse without SVD (not so precize) # A' = (A.t * A)^-1 * A.t pinv_matrix = Matrix[*(mat.transpose*mat).to_a].inverse * Matrix[*mat.to_a].transpose pinv = Matrix[*pinv_matrix.to_a] pinv.row(deriv).to_a end