所以，我尝试在 C# 中实现 scipy.signal.filtfilt 函数，但是我没有得到我在 python 中得到的结果，我看到大的瞬变或非常微弱的信号。我查看了 scipy 函数的来源并试图模仿它，但没有成功。
我直接从我的 python 脚本中获取了滤波器系数 a 和 b 以及初始状态，所以没有区别。与 filtfilt 中的方法类似，我通过将信号旋转 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;
}
最佳答案
因此，事实证明 IIR 滤波器对给定的参数和数值精度非常敏感。我从一个文本文件中导入了我的过滤器系数，我认为它已经足够好了，但事实并非如此。从二进制文件导入给出了我想要的结果，我的算法产生与 scipy filtfilt 相同的输出。
为了完整起见，用于从 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;
}
}
关于与 numpy.filtfilt 相当的 C# 数字滤波器实现，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/45076210/