Savitzky-Golay filtering

[alexreg]

Do you plan to implement it at some point? It's one of the most well-known and widely used filters, of course.

[dylanmeysmans]

@alexreg @cdrnet: I have both a naive C# and a naive F# implementation of the Savitzky-Golay filter using MathNet.Numerics.LinearAlgebra:

using System;
using MathNet.Numerics.LinearAlgebra;

namespace Filtering
{
    /// <summary>
    /// <para>Implements a Savitzky-Golay smoothing filter, as found in [1].</para>
    /// <para>[1] Sophocles J.Orfanidis. 1995. Introduction to Signal Processing. Prentice-Hall, Inc., Upper Saddle River, NJ, USA.</para>
    /// </summary>
    public sealed class SavitzkyGolayFilter
    {
        private readonly int sidePoints;

        private Matrix<double> coefficients;

        public SavitzkyGolayFilter(int sidePoints, int polynomialOrder)
        {
            this.sidePoints = sidePoints;
            Design(polynomialOrder);
        }

        /// <summary>
        /// Smoothes the input samples.
        /// </summary>
        /// <param name="samples"></param>
        /// <returns></returns>
        public double[] Process(double[] samples)
        {
            int length = samples.Length;
            double[] output = new double[length];
            int frameSize = (sidePoints << 1) + 1;
            double[] frame = new double[frameSize];

            Array.Copy(samples, frame, frameSize);

            for (int i = 0; i < sidePoints; ++i)
            {
                output[i] = coefficients.Column(i).DotProduct(Vector<double>.Build.DenseOfArray(frame));
            }

            for (int n = sidePoints; n < length - sidePoints; ++n)
            {
                Array.ConstrainedCopy(samples, n - sidePoints, frame, 0, frameSize);
                output[n] = coefficients.Column(sidePoints).DotProduct(Vector<double>.Build.DenseOfArray(frame));
            }

            Array.ConstrainedCopy(samples, length - frameSize, frame, 0, frameSize);

            for (int i = 0; i < sidePoints; ++i)
            {
                output[length - sidePoints + i] = coefficients.Column(sidePoints + 1 + i).DotProduct(Vector<double>.Build.Dense(frame));
            }

            return output;
        }

        private void Design(int polynomialOrder)
        {
            double[,] a = new double[(sidePoints << 1) + 1, polynomialOrder + 1];

            for (int m = -sidePoints; m <= sidePoints; ++m)
            {
                for (int i = 0; i <= polynomialOrder; ++i)
                {
                    a[m + sidePoints, i] = Math.Pow(m, i);
                }
            }

            Matrix<double> s = Matrix<double>.Build.DenseOfArray(a);
            coefficients = s.Multiply(s.TransposeThisAndMultiply(s).Inverse()).Multiply(s.Transpose());
        }
    }
}

open MathNet.Numerics.LinearAlgebra
open System

module SavitzkyGolay =
    let design polynomialOrder sidePoints =
        let width = (sidePoints * 2) + 1
        let s =
            Array.init width 
                (fun m -> 
                Array.init (polynomialOrder + 1) (fun i -> float (pown m i))) 
            |> matrix
        (width, 
         s.Multiply(s.TransposeThisAndMultiply(s).Inverse())
          .Multiply(s.Transpose()))
    
    let apply polynomialOrder sidePoints =
        let (width, coefficients) = design polynomialOrder sidePoints
        fun (signal : float []) -> 
            let start =
                Array.init (sidePoints + 1) 
                    (fun i -> 
                    coefficients.Column(i)
                                .DotProduct(vector signal.[..width - 1]))
            let end' =
                Array.init (sidePoints + 1) 
                    (fun i -> 
                    coefficients.Column(sidePoints + i)
                                .DotProduct(vector 
                                                signal.[Array.length signal 
                                                        - width..]))
            let steady =
                Array.init (Array.length signal - (2 * sidePoints + 3)) 
                    (fun i -> 
                    coefficients.Column(sidePoints + 1)
                                .DotProduct(vector 
                                                signal.[i..i + 2 * sidePoints]))
            Array.concat [| start; steady; end' |]

I ported my work from the MATLAB files sg.m and sgfilt.m found here. The Savitzky-Golay filters that ship with MATLAB and LabVIEW produce the same results as that implementation, as does mine. However, I do not know enough about this library or digital filter design to adapt my ports to it. Feel free to do so, though.

[juancaldone]

Hello,

I am trying to implement the Savitzky-Golay filter but I am having errors on the Process function, can you submit an example on how to use it please?.

Thank you!.

[dylanmeysmans]

@juancaldone You use the Savitzky-Golay filter with a signal of at least 2n + 1 samples long for a number of side points of the window n and a polynomial order less than n. For example:

using System;
using System.Linq;
					
public class Program
{
	public static void Main()
	{
	    var signal = Enumerable.Range(0, 100).Select(x => Math.Sin(x * 0.1));
	    var random = new Random();
	    var noise = Enumerable.Range(0, 100).Select(_ => random.NextDouble() * 0.01);
	    var noisySignal = signal.Zip(noise, (x, y) => x + y).ToArray();
            var filteredSignal = new SavitzkyGolayFilter(5, 1).Process(noisySignal);
	}
}

[juancaldone]

Thank you very much for your prompt response.

[xas]

@mettekou are you sure about your C# code ?
If I choose a filter with a side length of 20, I got 41 values (20 left, 20 right, 1 middle)
So your 1st for should go from index 0 to index 19 (20 values). But your i <= SidePoints will take 21 values
Same for the 2nd for, you should take the 21st column, so it should be Coefficients.Column(SidePoints) ?
Same for the 3rd for, as Coefficients.Columnshould begin at SidePoints + 1

PS: Except at the second for, you can put the Array.Copy outside the for

[dylanmeysmans]

@mettekou are you sure about your C# code ?
If I choose a filter with a side length of 20, I got 41 values (20 left, 20 right, 1 middle)
So your 1st for should go from index 0 to index 19 (20 values). But your i <= SidePoints will take 21 values
Same for the 2nd for, you should take the 21st column, so it should be Coefficients.Column(SidePoints) ?
Same for the 3rd for, as Coefficients.Columnshould begin at SidePoints + 1

@xas Indeed, this is a mistake I made when porting the filter from MATLAB to C#, as MATLAB uses one-based indexing for vectors and matrices. I discovered it myself a couple of months ago and forgot to update this version on GitHub. I have now done so; please find the updated version above.

[cc4566]

Can you add comments to the C# code to understand what are you doing? great work by the way

[dylanmeysmans]

Can you add comments to the C# code to understand what are you doing? great work by the way

@cc4566 Please consult the section on Savitzky-Golay smoothing filters in this freely available book (page 427 to 453) for an in-depth explanation, as my code was ported from that section's MATLAB code.

[KonradZaremba]

Works. Thanks a lot. Is there any plan for 2d implementation?

[dylanmeysmans]

Works. Thanks a lot. Is there any plan for 2d implementation?

@KonradZaremba No, I only use it for signal processing.

[emanuelhristea]

There is a performance penalty here:
output[n] = coefficients.Column(sidePoints).DotProduct(Vector.Build.DenseOfArray(frame));
Is there any reason to copy over each frame? it is copied earlier.

[dylanmeysmans]

There is a performance penalty here:
output[n] = coefficients.Column(sidePoints).DotProduct(Vector.Build.DenseOfArray(frame));
Is there any reason to copy over each frame? it is copied earlier.

@emanuelhristea No, I got rid of that over a year ago in my own codebase, but I do not update the snippet in this thread every time. @cdrnet A lot of people use this snippet, maybe it is time to finally add it to the library?

[gchen001]

@mettekou Thanks a lot for translating this library from Matlab codes into the C# MathNet version. Do you have any further plan for implementing the Savitzky Golay Derivatives in this code snippet? Thanks a lot

[dylanmeysmans]

@mettekou Thanks a lot for translating this library from Matlab codes into the C# MathNet version. Do you have any further plan for implementing the Savitzky Golay Derivatives in this code snippet? Thanks a lot

@gchen001 Sorry, no, I do not.

[ChristoWolf]

Hi!

@cdrnet: Is there any goal to merge this or provide this filter any other way?