From 9464474cd290608a33fd8e59f880baa2aa237032 Mon Sep 17 00:00:00 2001 From: Philip Constantinou Date: Sun, 10 May 2020 22:37:54 -0700 Subject: [PATCH] Initial comit --- Makefile | 13 ++++ README.md | 0 doc.go | 31 +++++++++ sg.go | 152 +++++++++++++++++++++++++++++++++++++++++++ sg_test.go | 184 +++++++++++++++++++++++++++++++++++++++++++++++++++++ 5 files changed, 380 insertions(+) create mode 100644 Makefile create mode 100644 README.md create mode 100644 doc.go create mode 100644 sg.go create mode 100644 sg_test.go diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..5d2e022 --- /dev/null +++ b/Makefile @@ -0,0 +1,13 @@ +all: get build test + +clean: + rm *.png *.out + +get: + go get . + +build: + go build + +test: + go test -coverprofile=coverage.out \ No newline at end of file diff --git a/README.md b/README.md new file mode 100644 index 0000000..e69de29 diff --git a/doc.go b/doc.go new file mode 100644 index 0000000..7ad38e7 --- /dev/null +++ b/doc.go @@ -0,0 +1,31 @@ +package savitzkygolay + +/* +savitzkygolay provides a filter on a set of data which provides +an effective way of smoothing data that generally follows +curves found in polynomials. It is particularly good alternative +to a moving average since it does not introduce a +delay proportial to about half the window length. + +Example: + + noise := 15.0 + xs := make([]float64, testSize, testSize) + ys = make([]float64, testSize, testSize) + for i := range xs { + ys[i] = 20 * math.Sin(float64(i)/math.Pi/6) + + (rand.Float64() * noise) - noise/2.0) + xs[i] = float64(i) + } + + filter, err := savitzkygolay.NewFilterWindow(11) + sgy, err := filter.Process(ys, xs) + +Filter interface may be retained to avoid the overhead of pre-computing +the polynomials however the size is proportial to the square of the +window size. + +The filter run on O(number of elements * size of window) + +Project unit tests generate outputs which illustrate the filter. +*/ diff --git a/sg.go b/sg.go new file mode 100644 index 0000000..26b785f --- /dev/null +++ b/sg.go @@ -0,0 +1,152 @@ +package savitzkygolay + +import ( + "fmt" + "math" +) + +// Filter provides an interface for filtering data based on a filter configuration +// data is an array of data being filtered +// x is the corresponding x-position of the coordinate +// errors occur only if the filter window is larger than the data +type Filter interface { + Process(data []float64, x []float64) ([]float64, error) +} + +// filterConfiguration provides configurations for the filter +type filterConfiguration struct { + // windowSize is the number of points in the window + windowSize int + // derivative + derivative int + // polynomial is the order of the polynomial used in the fit + polynomial int + + weights [][]float64 +} + +// NewFilter creates new savitzky golay based on the provided attributes +func NewFilter(windowSize int, derivitive int, polynomial int) (Filter, error) { + options := filterConfiguration{windowSize: windowSize, derivative: derivitive, polynomial: polynomial} + if options.windowSize%2 == 0 || options.windowSize < 5 { + return nil, fmt.Errorf("options.WindowSize [%d] must be odd and equal to or greater than 5", options.windowSize) + } + if options.derivative < 0 { + return nil, fmt.Errorf("options.Derivative [%d] must be euqal or greater than 0", options.derivative) + } + if options.polynomial < 0 { + return nil, fmt.Errorf("options.Polynomial [%d] must be equal or greater than 0", options.polynomial) + } + options.weights = options.computeWeights() + return options, nil +} + +// NewFilterWindow creates a new savitzky golay filter with default settings of a 3rd order polynomial +// windowSize must be odd and greater than 4 +func NewFilterWindow(windowSize int) (Filter, error) { + return NewFilter(windowSize, 0, 3) +} + +// Process executes a on the input set +func (options filterConfiguration) Process(data []float64, h []float64) ([]float64, error) { + if options.windowSize > len(data) { + return nil, fmt.Errorf("data length [%d] must be larger than options.WindowSize[%d]", len(data), options.windowSize) + } + + halfWindow := int(math.Floor(float64(options.windowSize) / 2.0)) + numPoints := len(data) + results := make([]float64, numPoints) + weights := options.weights + hs := 0.0 + + //For the borders + for i := 0; i < halfWindow; i++ { + wg1 := weights[halfWindow-i-1] + wg2 := weights[halfWindow+i+1] + d1 := 0.0 + d2 := 0.0 + for l := 0; l < options.windowSize; l++ { + d1 += wg1[l] * data[l] + d2 += wg2[l] * data[numPoints-options.windowSize+l] + } + hs = getHs(h, halfWindow-i-1, halfWindow, options.derivative) + results[halfWindow-i-1] = d1 / hs + hs = getHs(h, numPoints-halfWindow+i, halfWindow, options.derivative) + results[numPoints-halfWindow+i] = d2 / hs + } + + //For the internal points + wg := weights[halfWindow] + for i := options.windowSize; i <= numPoints; i++ { + d := 0.0 + for l := 0; l < options.windowSize; l++ { + d += wg[l] * data[l+i-options.windowSize] + } + hs = getHs(h, i-halfWindow-1, halfWindow, options.derivative) + results[i-halfWindow-1] = d / hs + } + return results, nil +} + +func getHs(h []float64, center int, half int, derivative int) float64 { + hs := 0.0 + count := 0 + for i := center - half; i < center+half; i++ { + if i >= 0 && i < len(h)-1 { + hs += h[i+1] - h[i] + count++ + } + } + return math.Pow(hs/float64(count), float64(derivative)) +} + +func gramPolynomial(i int, m int, k int, s int) float64 { + result := 0.0 + if k > 0 { + result = + float64(float64(4*k-2)/float64(k*(2*m-k+1)))* + (float64(i)*gramPolynomial(i, m, k-1, s)+float64(s)*gramPolynomial(i, m, k-1, s-1)) - + (float64((k-1)*(2*m+k))/float64(k*(2*m-k+1)))* + gramPolynomial(i, m, k-2, s) + } else { + if k == 0 && s == 0 { + result = 1 + } else { + result = 0 + } + } + return result +} + +func productOfRange(a, b int) int { + gf := 1 + if a >= b { + for j := a - b + 1; j <= a; j++ { + gf *= j + } + } + return gf +} + +func polyWeight(i, t, windowMiddle, polynomial, derivitive int) float64 { + sum := 0.0 + for k := 0; k <= polynomial; k++ { + sum += + float64(2*k+1) * + (float64(productOfRange(2*windowMiddle, k)) / float64(productOfRange(2*windowMiddle+k+1, k+1))) * + gramPolynomial(i, windowMiddle, k, 0) * gramPolynomial(t, windowMiddle, k, derivitive) + } + return sum +} + +func (options *filterConfiguration) computeWeights() [][]float64 { + weights := make([][]float64, options.windowSize) + windowMiddle := int(math.Floor(float64(options.windowSize) / 2.0)) + for row := -windowMiddle; row <= windowMiddle; row++ { + weights[row+windowMiddle] = make([]float64, options.windowSize) + for col := -windowMiddle; col <= windowMiddle; col++ { + weights[row+windowMiddle][col+windowMiddle] = polyWeight(col, row, windowMiddle, options.polynomial, options.derivative) + } + } + return weights +} diff --git a/sg_test.go b/sg_test.go new file mode 100644 index 0000000..5e92549 --- /dev/null +++ b/sg_test.go @@ -0,0 +1,184 @@ +package savitzkygolay + +import ( + "image/color" + "math" + "math/rand" + "testing" + + assert "github.com/stretchr/testify/assert" + plot "gonum.org/v1/plot" + plotter "gonum.org/v1/plot/plotter" +) + +const testSize = 500 + +type xyPairs struct { + xs []float64 + ys []float64 +} + +func (p xyPairs) Len() int { + v := p.xs + return len(v) +} + +func (p xyPairs) XY(i int) (x, y float64) { + return p.xs[i], p.ys[i] +} + +func makexyPairs(size int) xyPairs { + var r xyPairs + r.xs = make([]float64, testSize, testSize) + r.ys = make([]float64, testSize, testSize) + return r +} + +func Test_SavitzkyGolay_Args(t *testing.T) { + _, err := NewFilterWindow(0) + assert.Error(t, err, "Window size too small") + _, err = NewFilterWindow(3) + assert.Error(t, err, "Window size too small") + _, err = NewFilterWindow(6) + assert.Error(t, err, "Window size even") + + _, err = NewFilter(7, -1, 3) + assert.Error(t, err, "Derivitive must be non-negative") + + _, err = NewFilter(7, 0, -1) + assert.Error(t, err, "Polynomial must be non-negative") + + f, err := NewFilter(7, 0, 3) + assert.NoError(t, err, "Filter should be allowed") + xs := []float64{1, 2, 3, 4, 5} + _, err = f.Process(xs, xs) + assert.Error(t, err, "Window larger than data") + +} + +func Test_SavitzkyGolay_Line(t *testing.T) { + pairs := makexyPairs(testSize) + for i := range pairs.xs { + pairs.ys[i] = math.Pi + pairs.xs[i] = float64(i) + } + + filter, err := NewFilterWindow(5) + assert.NoError(t, err, "No filter initialization error expected") + sgy, err := filter.Process(pairs.ys, pairs.xs) + assert.NoError(t, err, "No error expected") + copy := pairs + copy.ys = sgy + + max, avg := pairs.difference(©) + assert.NoError(t, err, "No error expected") + assert.Less(t, avg, float64(0.1), "Small average differences") + assert.Less(t, max, float64(0.5), "Small average differences") + + visualTest(pairs, nil, copy, "line_smoothing.png") +} + +func Test_SavitzkyGolay_Sign(t *testing.T) { + pairs := makexyPairs(testSize) + for i := range pairs.xs { + pairs.ys[i] = 20 * math.Sin(float64(i)/math.Pi/6) + pairs.xs[i] = float64(i) + } + copy := pairs + + filter, err := NewFilterWindow(41) + assert.NoError(t, err, "No filter initialization error expected") + copy.ys, err = filter.Process(pairs.ys, pairs.xs) + assert.NoError(t, err, "No error expected") + + max, avg := pairs.difference(©) + assert.NoError(t, err, "No error expected") + assert.Less(t, avg, float64(0.1), "Small average differences") + assert.Less(t, max, float64(0.5), "Small average differences") + + visualTest(pairs, nil, copy, "sin_smoothing.png") +} + +func noise(size float64) float64 { + return (rand.Float64() * size) - size/2 +} + +func Test_SavitzkyGolay_SignNoise(t *testing.T) { + pairs := makexyPairs(testSize) + for i := range pairs.xs { + pairs.ys[i] = 20 * math.Sin(float64(i)/math.Pi/6) + pairs.xs[i] = float64(i) + } + + noisy := pairs.addNoise(5.0) + + filter, err := NewFilterWindow(21) + assert.NoError(t, err, "No filter initialization error expected") + sgy, err := filter.Process(noisy.ys, noisy.xs) + assert.NoError(t, err, "No error expected") + + copy := pairs + copy.ys = sgy + visualTest(pairs, noisy, copy, "sin_noise_smoothing.png") +} + +func (p *xyPairs) addNoise(n float64) xyPairs { + noisy := *p + noisy.ys = make([]float64, len(p.ys)) + for i, y := range p.ys { + noisy.ys[i] = y + noise(n) + } + return noisy +} + +func Test_SavitzkyGolay_SignNoise_1(t *testing.T) { + pairs := makexyPairs(testSize) + for i := range pairs.xs { + pairs.ys[i] = 20*math.Sin(float64(i)/math.Pi/4) + + 20*math.Sin(float64(i+10)/math.Pi/2) + pairs.xs[i] = float64(i) + } + noisy := pairs.addNoise(5.0) + + filter, err := NewFilterWindow(15) + assert.NoError(t, err, "No filter initialization error expected") + sgy, err := filter.Process(pairs.ys, pairs.xs) + assert.NoError(t, err, "No error expected") + + copy := pairs + copy.ys = sgy + visualTest(pairs, noisy, copy, "sin_noise_smoothing_order_2.png") +} + +func visualTest(original, noisy, filtered plotter.XYer, path string) { + p, _ := plot.New() + if original != nil { + line, _ := plotter.NewLine(original) + line.Color = color.Gray{Y: 128} + p.Add(line) + } + if noisy != nil { + scatter, _ := plotter.NewScatter(noisy) + p.Add(scatter) + } + + if filtered != nil { + green := color.RGBA{G: 255} + line, _ := plotter.NewLine(filtered) + line.Color = green + p.Add(line) + } + + _ = p.Save(1512, 512, path) +} + +func (p *xyPairs) difference(o *xyPairs) (max, average float64) { + for i, v := range p.ys { + oy := o.ys[i] + d := math.Abs(oy - v) + max = math.Max(max, d) + average += d + } + average = average / float64(len(p.ys)) + return max, average +}