基于C#的窗函数法低通FIR滤波器

基于C#的窗函数法低通FIR滤波器

摘要

基于Visual Studio 2015 开发环境，使用C#编程语言，运用窗函数法构造了FIR低通滤波器。并通过MATLAB对其滤波效果进行了试验。
由于凯泽窗的数学模型过于复杂，笔者能力有限，故无法在本文中给出凯泽窗的窗函数代码。

源代码
1、矩形窗函数

namespace Low_Pass_FIR_Fliter
{
    class Boxcar
    {
        public int N = 0;
        public Boxcar (double Wp, double Ws)
        {
            int i;
            double n = (0.9 * 2 * Math.PI) / (Ws - Wp);//计算滤波器长度
            for (i = 0; i < n; i++) ;
            N = i;//得到滤波器长度
        }
        public double[] Get_Win()//返回N点矩形窗序列
        {
            int n;
            double[] wd = new double[N];
            for (n = 0; n < N; n++)
            {
                wd[n] = 1;
            }
            return wd;
        }
    }
}

2、三角窗函数

namespace Low_Pass_FIR_Fliter
{
    class Bartlett
    {
        public int N = 0;
        public Bartlett (double Wp, double Ws)
        {
            int i;
            double n = (2.1 * 2 * Math.PI) / (Ws - Wp);//计算滤波器长度
            for (i = 0; i < n; i++) ;
            N = i;//得到滤波器长度
        }
        public double[] Get_Win()//返回N点三角窗序列
        {
            int n;
            double[] wd = new double[N];
            for (n = 0; n < N; n++)
            {
                if(n<=((N-1)/2))
                {
                    wd[n] = (2 * n / (N - 1));
                }
                if(n > ((N - 1) / 2))
                {
                    wd[n] = 2 - (2 * n / (N - 1));
                }
            }
            return wd;
        }
    }

3、汉宁窗

namespace Low_Pass_FIR_Fliter
{
    class Hanning
    {
        public int N = 0;
        public Hanning(double Wp, double Ws)
        {
            int i;
            double n = (3.1 * 2 * Math.PI) / (Ws - Wp);//计算滤波器长度
            for (i = 0; i < n; i++) ;
            N = i;//得到滤波器长度
        }
        public double[] Get_Win()//返回N点汉宁窗序列
        {
            int n;
            double[] wd = new double[N];
            for (n = 0; n < N; n++)
            {
                double b = Math.Cos((2 * Math.PI*n) / (N - 1));
                double res = (0.5 - 0.5 * b);
                wd[n] = res;
            }
            return wd;
        }
    }
}

4、汉明窗

namespace Low_Pass_FIR_Fliter
{
    public class Hamming
    {
        public int N = 0;
        public Hamming(double Wp,double Ws)
        {
            int i;
            double n = (3.3 * 2 * Math.PI) / (Ws-Wp);//计算滤波器长度
            for (i = 0; i < n; i++) ;
            N = i ;//得到滤波器长度
        }
        public double [] Get_Win()//返回N点汉明窗序列
        {
            int n ;
            double[] wd=new double[N];
            for(n=0;n<N;n++)
            {
                double b = Math.Cos((2*n * Math.PI) / (N-1));
                double res = (0.54 - 0.46 * b);
                wd[n] = res;
            }
            return wd;
        }
    }
}

5、布莱克曼窗

namespace Low_Pass_FIR_Fliter
{
    class Blackman
    {
        public int N = 0;
        public Blackman (double Wp, double Ws)
        {
            int i;
            double n = (5.5 * 2 * Math.PI) / (Ws - Wp);//计算滤波器长度
            for (i = 0; i < n; i++) ;
            N = i;//得到滤波器长度
        }
        public double[] Get_Win()//返回N点布莱克曼窗序列
        {
            int n;
            double[] wd = new double[N];
            for (n = 0; n < N; n++)
            {
                double a = Math.Cos((2 * n * Math.PI) / (N - 1));
                double b = Math.Cos((4 * n * Math.PI) / (N - 1));
                double res = (0.42 - 0.5 * a + 0.08 * b);
                wd[n] = res;
            }
            return wd;
        }
    }
}

6、单位冲激响应

namespace Low_Pass_FIR_Fliter
{
    class UnitImpulseReact
    {
        private double Wc;
        private int alpha;
        private int N;
        public UnitImpulseReact(double Wp,double Ws,int N)//构造单位冲激响应
        {
            this.Wc = 0.5 * (Wp + Ws);//取近似截止频率
            if (N % 2 != 0)
            {
                this.alpha = (N - 1) / 2;//计算相移常数
            }
            else
            {
                this.alpha = N / 2;
            }
            this.N = N;
        }
        public double[] Get_UIR()//获取单位冲激响应序列
        {
            double[] hd = new double[N];
            for(int n=0;n<N;n++)
            {
                double numerator = Math.Sin(Wc * (n - alpha));
                double denominator = Math.PI * (n - alpha);
                if (n == alpha)//当n=α时，取极限
                {
                    hd[n] = Wc / Math.PI;
                }
                else
                {
                    hd[n] = numerator / denominator;
                }
            }
            return hd;
        }
    }
}

使用MATLAB做正确性检验

大致思路：自定义一个15点序列，通过频谱图确定技术指标。分别通过MATLAB标准做法和自编程序对其进行滤波，若滤波结果大致一致，则可认为这个滤波器程序是正确的。

自定义序列 x[n]={ 16 12 2 -4 -1 1 -3 11 -27 18 15 -9 11 21 3} 该序列对应的频谱图为：

选取通带边界频率ωp=0.2π，阻带边界频率ωs=0.3π，阻带衰减>30dB。由此，选择汉宁窗。通过MATLAB程序构成的低通滤波器b1的频谱图为：

滤波后的序列y的频率特性为（将>0dB的部分局部放大）

同样的，在程序中基于预定的滤波器特性构建汉宁窗，并结合单位冲激响应得到滤波器的单位样值响应

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using System.IO;

namespace Low_Pass_FIR_Fliter
{
    class Program
    {
        static void Main(string[] args)
        {
            Hanning win = new Hanning(0.2 * Math.PI, 0.3 * Math.PI);//构造汉宁窗
            UnitImpulseReact UIR = new UnitImpulseReact(0.2 * Math.PI, 0.3 * Math.PI, win.N);//构造单位冲激响应
            double[] hd = new double[win.N];
            double[] w = new double[win.N];
            double[] h = new double[win.N];
            string[] res = new string[win.N];//定义数组
            w = win.Get_Win();//获取窗函数序列
            hd = UIR.Get_UIR();//获取单位冲激响应序列
            for(int i=0;i<win.N;i++)
            {
                h[i] = w[i] * hd[i];//得到滤波器序列
            }
            for (int i = 0; i < win.N; i++)
            {
                res[i] = Convert.ToString(h[i]);//将滤波器序列点数转化成字符串形式，准备写入文件
            }
            File.WriteAllLines(@"C:\VS2015\C#Programe\Low-Pass_FIR_Fliter\Low-Pass_FIR_Fliter\sequence.txt", res);//写入文件，用于MATLAB的后续处理
        }
    }
}

通过程序获得的滤波器b2的频率特性为：

程序滤波器滤波后的序列y1的频率特性为（同样将>0dB的部分局部放大）

通过波形比较可以认为程序生成的滤波器是成功的。