C#零相位数字滤波器，改写自Matlab函数filtfilt

任何一个数字滤波器都有幅频特性和相频特性，如果对于滤波不要求实时性，我们可以设计一种滤波器，使得它的相频特性始终为0，这种数字滤波器就称为零相移数字滤波器。

在Matlab中，零相移滤波器对应的函数名称为filtfilt，其帮助中提供的信息如下，翻译一下：

FILTFILT Zero-phase forward and reverse digital filtering.

    Y = FILTFILT(B, A, X) filters the data in vector X with the filter described by vectors A and B to create the filtered data Y.  The filter is described by the difference equation:

FILTFILT 零相移前向后向数字滤波器

         Y = FILTFILT(B, A, X)通过向量A，B描述的旅欧不起对数据向量X滤波得到滤波后的数据Y。滤波器可以通过差分方程描述如下：

      y(n) = b(1)*x(n) + b(2)*x(n-1) + ... + b(nb+1)*x(n-nb)

                       - a(2)*y(n-1) - ... - a(na+1)*y(n-na)

After filtering in the forward direction, the filtered sequence is then reversed and run back through the filter; Y is the time reverse of the output of the second filtering operation. The result has precisely zero phase distortion and magnitude modified by the square of the filter's magnitude response.  Care is taken to minimize startup and ending transients by matching initial conditions.

通过前向滤波之后，滤波序列被反折并且再次通过滤波器，Y是第二次滤波器输出的时间反折序列。其结果是相位特性精确的为零，没有变形并且幅频特性被滤波器幅频特性的平方所改变。初值和终值都被选择的很小的跳变来匹配初始条件。

The length of the input x must be more than three times the filter order, defined as max(length(b)-1,length(a)-1). Note that FILTFILT should not be used with differentiator and Hilbert FIR filters, since the operation of these filters depends heavily on their phase response.

输入序列x的长度必须大于滤波器阶数的三倍,阶数定义为max(length(b)-1,length(a)-1)。注意FILTFILT不应当被用于Hilbert FIR滤波器，因为这些滤波器的相位特性是很重要的。

考虑到这个滤波器的科研应用价值，为便于将它引入到自己的分析软件中，于是调试将其matlab代码改写为了C#代码。源码如下：

public class DigitalFilter

    {

        private double[] b = new double[2];

        private double[] a = new double[2];

        private double[] x;

        private double zi;

        private double[] y = new double[1024];

        private double[] zf = new double[1];

 

        public DigitalFilter(double[] b, double[] a, double[] x, double zi)

        {

            this.b = b;

            this.a = a;

            this.x = x;

            this.zi = zi;

        }

        private double[] getY()

        {

            calc();

            return y;

        }

        private double[] getZf()

        {

            calc();

            return zf;

        }

        private void calc()

        {

            for (int i = 0; i < y.Length; i++)

            {

                if (i == 0)

                {

                    y[i] = b[0] * x[i] + zi;

                }

                else

                {

                    y[i] = b[0] * x[i] + b[1] * x[i - 1] - a[1] * y[i - 1];

                    if (i == x.Length - 1)

                    {

                        zf[0] = y[i];

                    }

                }

            }

        }

        public double[] zeroFilter()

        {

            int len = x.Length;    // length of input

            int nb = b.Length;

            int na = a.Length;

            int nfilt = Math.Max(na, nb);

            int nfact = 3 * (nfilt - 1);  // length of edge transients

            //运算初值

            double data = 1 + a[1];

            double zi;

            zi = (b[1] - a[1] * b[0]) / data;

            //首尾添数

            double[] yTemp = new double[y.Length + 2 * nfact];

            for (int i = 0; i < nfact; i++)

            {

                yTemp[i] = 2 * x[0] - x[nfact - i];

            }

            for (int i = nfact; i < y.Length + nfact; i++)

            {

                yTemp[i] = x[i - nfact];

            }

            for (int i = y.Length + nfact;i<yTemp.Length ; i++)

            {

                yTemp[i] = 2 * x[x.Length - 1] - x[yTemp.Length- 2 - i + y.Length-nfact];

            }

            //正向滤波

            this.zi = zi * yTemp[0];

            yTemp = zeroCalc(yTemp);

            //反序

            yTemp=this.reverse(yTemp);

            //反向滤波

            this.zi = zi * yTemp[0];

            yTemp = zeroCalc(yTemp);

            //反序

            yTemp = this.reverse(yTemp);

            for (int i = 0; i < y.Length; i++)

            {

                y[i] = yTemp[i + nfact];

            }

            return y;

        }

        private double[] zeroCalc(double[] xx)

        {

            double [] yy = new double [xx.Length];

            for (int i = 0; i < yy.Length; i++)

            {

                if (i == 0)

                {

                    yy[i] = b[0] * xx[i] + zi;

                }

                else

                {

                    yy[i] = b[0] * xx[i] + b[1] * xx[i - 1] - a[1] * yy[i - 1];

                }

            }

            return yy;

        }

        private double[] reverse(double[] data)

        {

            double tmp;

            for (int i = 0; i < data.Length / 2; i++)

            {

                tmp = data[data.Length - i - 1];

                data[data.Length - i - 1] = data[i];

                data[i] = tmp;

            }

            return data;

        }

    }

转载于:https://www.cnblogs.com/oyjj/archive/2009/09/15/2132977.html