C# digital filter implementation comparable to numpy.filtfilt

Question

So, I try to implement the scipy.signal.filtfilt function in C#, however I do not get the results I get in python, either I see large transients or a very weak signal. I looked at the source of the scipy functions and tried to mimic it, however unsuccessfull.

I took the filter coefficents a and b and the initial state directly from my python script, so there is no difference. Similar to the approach in filtfilt I pad the data at both ends with 3 times the filter size by rotating the signal by 180°(Also done here MATLAB's filtfilt() Algorithm). The initial state is multiplied by the first element of my padded data. The filter itself calculates the new output and all the states d[] for the new input data.

So, does anybody know where my error is in my implementation?

Here is my code:

/// <remarks>
/// The filter function is implemented as a direct II transposed structure. This means that the filter implements:
///
///a[0]*y[n] = b[0]*x[n] + b[1]*x[n - 1] + ... + b[M]*x[n - M]
///              - a[1]*y[n - 1] - ... - a[N]*y[n - N]
///where M is the degree of the numerator, N is the degree of the denominator, and n is the sample number.It is implemented using the following difference equations(assuming M = N):
///
///a[0]* y[n] = b[0] * x[n] + d[0][n - 1]
///d[0][n] = b[1] * x[n] - a[1] * y[n] + d[1][n - 1]
///d[1][n] = b[2] * x[n] - a[2] * y[n] + d[2][n - 1]
///...
///d[N - 2][n] = b[N - 1]* x[n] - a[N - 1]* y[n] + d[N - 1][n - 1]
///d[N - 1][n] = b[N] * x[n] - a[N] * y[n]
///d[N] = 0</remarks>
public static double[] lfilter(double[] x, double[] a, double[] b, double[] zi)
{
    if (a.Length != b.Length)
        throw new ArgumentOutOfRangeException("A and B filter coefficents should have the same length");
    double[] y = new double[x.Length];
    int N = a.Length;
    double[] d = new double[N];
    zi.CopyTo(d, 0);
    for (int n = 0; n < x.Length; ++n)
    {
        y[n] = b[0] * x[n] + d[0];
        for(int f = 1; f < N; ++f)
        {
            d[f - 1] = b[f] * x[n] - a[f] * y[n] + d[f];
        }
    }
    return y;
}


/// <summary>
/// Apply a digital filter forwards and backwards to have a zero phase shift filter
/// </summary>
/// <param name="data">The data to filter</param>
/// <param name="a">Filter coefficents a</param>
/// <param name="b">Filter coefficents b</param>
/// <returns>The filterd data</returns>
/// <remarks>
/// In order to prevent transients at the end or start of the sequence we have to pad it
/// The padding is done by rotating the sequence by 180° at the ends and append it to the data
/// </remarks>
public static double[] FiltFilt(double[] data, double[] a, double[] b, double[] zi)
{
    //Pad the data with 3 times the filter length on both sides by rotating the data by 180°
    int additionalLength = a.Length * 3;
    int endOfData = additionalLength + data.Length;
    double[] x = new double[data.Length + additionalLength * 2];
    Array.Copy(data, 0, x, additionalLength, data.Length);
    for (int i = 0; i < additionalLength; i++)
    {
        x[additionalLength - i - 1] = (x[additionalLength] * 2) - x[additionalLength + i + 1];
        x[endOfData + i] = (x[endOfData - 1] * 2) - x[endOfData - i - 2];
    }
    //Calculate the initial values for the given sequence
    double[] zi_ = new double[zi.Length];
    for (int i = 0; i < zi.Length; i++)
        zi_[i] = zi[i] * x[0];
    double[] y = lfilter(x, a, b, zi_);
    double[] y_flipped = new double[y.Length];
    //reverse the data and filter again
    for (int i = 0; i < y_flipped.Length; i++)
    {
        y_flipped[i] = y[y.Length - i - 1];
    }
    for (int i = 0; i < zi.Length; i++)
        zi_[i] = zi[i] * y_flipped[0];
    y = lfilter(y_flipped, a, b, zi_);
    y_flipped = new double[data.Length];
    //rereverse it again and return
    for (int i = 0; i < y_flipped.Length; i++)
    {
        y_flipped[i] = y[endOfData - i - 1];
    }
    return y_flipped;
}

score 2 · Answer 1 · answered Jul 13 '17 at 11:26

So, as it turns out IIR filters are very sensitve to the parameters given and the numerical accuracy. I imported my filter coefficents from a text file, which I thought was good enough, however it wasn't. Importing from a binary file gives the result I want and my algorithm produces the same output as the scipy filtfilt.

For completness the functions used to export the parameters from python and import them into C#:

import sys
import os;
import numpy as np
import struct

import scipy.signal as signal

order = 4
lowF = 0.2
highF = 0.8

b, a = signal.butter(order, [lowF, highF], btype='band')

#b, a = signal.iirfilter(float(sys.argv[4]), float(sys.argv[5]), btype='lowpass')

path = os.path.dirname(os.path.abspath(sys.argv[0])) + '\\';

# Get the steady state of the filter's step response.
zi = signal.lfilter_zi(b, a)

f = open(path + 'parameterfile.bin',"wb") 

def writeToFile(array, file):
    file.write(struct.pack("<I", array.shape[0]))
    print('Byte length: ' + str(len(struct.pack("<I", array.shape[0]))))
    for i in range(array.shape[0]):
        file.write(bytes(array[i]));

writeToFile(a, f)
writeToFile(b, f)
writeToFile(zi, f)

f.close();

And the C# code:

private static void GetFilterAndZiFromBin(string path, out double[] a, out double[] b, out double[] zi)
{
    try
    {
        List<double> a_ = new List<double>();
        List<double> b_ = new List<double>();
        List<double> zi_ = new List<double>();
        FileStream fs = new FileStream(path,
                            FileMode.Open,
                            FileAccess.Read);
        BinaryReader br = new BinaryReader(fs);
        int length = br.ReadInt32();
        Console.WriteLine("Read " + length.ToString() + " doubles for a from file");
        while(length > 0)
        {
            a_.Add(br.ReadDouble());
            length--;
        }
        length = br.ReadInt32();
        Console.WriteLine("Read " + length.ToString() + " doubles for b from file");
        while (length > 0)
        {
            b_.Add(br.ReadDouble());
            length--;
        }
        length = br.ReadInt32();
        Console.WriteLine("Read " + length.ToString() + " doubles for zi from file");
        while (length > 0)
        {
            zi_.Add(br.ReadDouble());
            length--;
        }
        a = a_.ToArray();
        b = b_.ToArray();
        zi = zi_.ToArray();
    }
    catch (Exception e)
    {
        a = new double[0];
        b = new double[0];
        zi = new double[0];
        throw e;
    }

}

C# digital filter implementation comparable to numpy.filtfilt

1 Answers1