Open alexreg opened 8 years ago
@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.
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!.
@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);
}
}
Thank you very much for your prompt response.
@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.Column
should begin at SidePoints + 1
PS: Except at the second for
, you can put the Array.Copy
outside the for
@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 youri <= SidePoints
will take 21 values Same for the 2ndfor
, you should take the 21st column, so it should beCoefficients.Column(SidePoints)
? Same for the 3rdfor
, asCoefficients.Column
should begin atSidePoints + 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.
Can you add comments to the C# code to understand what are you doing? great work by the way
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.
Works. Thanks a lot. Is there any plan for 2d implementation?
Works. Thanks a lot. Is there any plan for 2d implementation?
@KonradZaremba No, I only use it for signal processing.
There is a performance penalty here:
output[n] = coefficients.Column(sidePoints).DotProduct(Vector
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?
@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
@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.
Hi!
@cdrnet: Is there any goal to merge this or provide this filter any other way?
Works. Thanks a lot. Is there any plan for 2d implementation?
Did you compare the results with matlab results? I get different
Works. Thanks a lot. Is there any plan for 2d implementation?
Did you compare the results with matlab results? I get different
I compared the results for my implementation with the implementation in LabVIEW, which yields results identical to the implementation in MATLAB.
Do you plan to implement it at some point? It's one of the most well-known and widely used filters, of course.