C＃數字濾波器的實現與numpy.filtfilt相媲美

Question

所以，我嘗試在C＃中實現scipy.signal.filtfilt函數，但是我沒有得到我在Python中得到的結果，要麼我看到大的瞬態信號或非常微弱的信號。我查看了scipy函數的來源，並試圖模仿它，但不成功。

我從我的python腳本中直接獲取了過濾器係數a和b以及初始狀態，所以沒有區別。類似於filteredfilt中的方法，通過將信號旋轉180°（在此也進行MATLAB's filtfilt() Algorithm），將數據兩端的數據填充3倍濾波器大小。初始狀態乘以填充數據的第一個元素。 過濾器本身爲新輸入數據計算新輸出和所有狀態d []。

那麼，有沒有人知道我的錯誤在我的執行？

這裏是我的代碼：

/// <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; 
} 

Answer 1

所以，事實證明IIR濾波器非常sensitve給出的參數和數值精度。我從一個文本文件中導入了我的過濾器係數，我認爲它足夠好，但事實並非如此。從二進制文件導入給出我想要的結果，我的算法產生與scipy filtfilt相同的輸出。

對於completness功能用於將參數從蟒蛇導出，並將其導入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(); 

和C＃代碼：

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; 
    } 

}