C# Moving Window Filter Example

using System;
using System.IO;

using CenterSpace.NMath.Core;

namespace CenterSpace.NMath.Stats.Examples.CSharp
{
  /// <summary>
  /// A .NET example in C# showing smoothing using the Savitzky-Golay moving filter.
  /// </summary>
  class MovingWindowFilterExample
  {

    static void Main(string[] args)
    {
      int signalLength = 2000;

      Console.WriteLine();

      DoubleVector signal = NoisySignal(signalLength);
      double noisySignalVariance = NMathFunctions.Variance(signal);
      Console.WriteLine("Noisy signal variance = " + noisySignalVariance);

      // Set up a moving average filter with an asymmetric window. A moving window
      // filter replaces each input data point with a linear combination of 
      // its surrounding points.  A filter is thus described by the number of
      // of points to the left and right of the input point and the coefficients
      // of the linear combination.
      // This filter will replace each input data point with the average of 
      // its value with the 4 values to the left, and the 5 values to the right.
      // Thus the coeffients to use are all equal to one over the number of points 
      // in the window, 10 in this case.
      int numberLeft = 4;
      int numberRight = 5;
      DoubleVector filterCoefficients = MovingWindowFilter.MovingAverageCoefficients(numberLeft, numberRight);
      MovingWindowFilter movingAverageFilter = new MovingWindowFilter(numberLeft, numberRight, filterCoefficients);

      // Filter the signal. Note that we must specify a boundary option. When
      // replacing the first input value, we don't have any points to the left, similarly, when 
      // replacing the last input value, we don't have any point to the right. The "PadWithZeros"
      // boundary option prepends "numberLeft" zeros and appends "numberRight" zeros to the 
      // input vector to deal with this.
      DoubleVector filteredSignal = movingAverageFilter.Filter(signal, MovingWindowFilter.BoundaryOption.PadWithZeros);
      double filteredSignalVariance = NMathFunctions.Variance(filteredSignal);
      Console.WriteLine("Moving Average: filtered signal variance = " + filteredSignalVariance);

      // Set up a Savitzky-Golay filter. This is a smoothing filter that replaces input values with
      // the value of a polynomial of specified degree fit through the input value and its
      // surrounding points. A least squares algorithm is used to determine the fitting polynomial.
      int degree = 4;
      filterCoefficients = MovingWindowFilter.SavitzkyGolayCoefficients(numberLeft, numberRight, degree);
      MovingWindowFilter savitzkyGolayFilter = new MovingWindowFilter(numberLeft, numberRight, filterCoefficients);

      // Filter the signal. Here we use a different boundary option: "DoNotFilterBoundaryPoints". 
      // This option will not filter or replace the first "numberLeft" or last "numberRight"
      // values of the input signal.
      filteredSignal = savitzkyGolayFilter.Filter(signal, MovingWindowFilter.BoundaryOption.DoNotFilterBoundaryPoints);

      filteredSignalVariance = NMathFunctions.Variance(filteredSignal);
      Console.WriteLine("Savitzky-Golay: filtered signal variance = " + filteredSignalVariance);
      Console.WriteLine();

      // If you are filtering lots of signals with the same length, it is more
      // economic to use the "Filter" method which allows you to specify
      // the vector to place the output filtered signal in. This avoids having
      // to allocate potentially large vectors on every call to "Filter".
      int numSignals = 10;
      DoubleMatrix noisySignals = NoisySignals(signalLength, numSignals);
      DoubleVector y = new DoubleVector(signal.Length);
      for (int i = 0; i < numSignals; i++)
      {
        DoubleVector s = noisySignals.Col(i);
        Console.WriteLine(string.Format("Noisy signal {0} variance = {1}", i, NMathFunctions.Variance(s)));
        savitzkyGolayFilter.Filter(s, MovingWindowFilter.BoundaryOption.PadWithZeros, ref y);
        Console.WriteLine(string.Format("Savitzky-Golay filtered signal {0} variance = {1}", i, NMathFunctions.Variance(y)));
        Console.WriteLine();
      }

      Console.WriteLine();
      Console.WriteLine("Press Enter Key");
      Console.Read();
    }

    static DoubleVector NoisySignal(int length)
    {
      RandGenNormal rng = new RandGenNormal();
      DoubleVector signal = new DoubleVector(length);
      for (int i = 0; i < length; i++)
      {
        signal[i] = Math.Cos(.2 * i) + rng.Next();
      }
      return signal;
    }

    static DoubleMatrix NoisySignals(int rows, int columns)
    {
      RandGenNormal rng = new RandGenNormal();
      DoubleMatrix signals = new DoubleMatrix(rows, columns);
      for (int i = 0; i < rows; i++)
      {
        for (int j = 0; j < columns; j++)
        {
          signals[i, j] = Math.Cos(.2 * i * j) + rng.Next();
        }
      }
      return signals;
    }
  }
}