C＃デジタルフィルタの実装numpy.filtfiltに匹敵

Question

私はscipy.signal.filtfilt関数をC＃で実装しようとしていますが、Pythonで得られる結果は得られません。大きなトランジェントや非常に弱い信号。私はscipy関数のソースを見て、それを模倣しようとしましたが、失敗しました。

私はpythonスクリプトからフィルタの係数aとbと初期状態を直接取​​ったので、違いはありません。 filtfiltの手法と同様に、信号を180°回転させることで、フィルターサイズの3倍のデータを両端にパッドします（MATLAB's filtfilt() Algorithmも同様）。初期状態には、埋め込まれたデータの最初の要素が乗算されます。 フィルタ自体は、新しい入力データの新しい出力とすべての状態d []を計算します。

私の実装でどこにエラーがありますか？結局のところとしてIIRフィルタは、指定されたパラメータと数値の精度に非常にsensitveあり、

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

ので：

は、ここに私のコードです。私はテキストファイルからフィルタ係数をインポートしましたが、私は十分だと思いましたが、そうではありませんでした。バイナリファイルからのインポートは私が望む結果を与え、私のアルゴリズムはscipy filtfiltと同じ出力を生成します。 completnessについては

のpythonからパラメータをエクスポートし、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; 
    } 

}