【FIR窗函数滤波器】C语言实现并与Matlab window+fir1对比

已于 2023-11-03 10:59:14 修改 · 1.3k 阅读

文章标签：

#c语言 #matlab #信号处理

于 2023-11-03 10:54:43 首次发布

文章详细介绍了C语言实现的低通、高通、带通、带阻滤波器的理想单位冲击响应函数，以及各种窗函数（如矩形、三角、汉宁、汉明、布莱克曼等）的实现，以及Matlab代码示例，展示了如何使用这些函数设计FIR滤波器。

参考文档

主要参考资料来源数字信号处理基于计算机的方法 4版（中文版）和Matlab窗函数【Matlab帮助中心】加窗法两部分

C语言实现

理想单位冲击响应函数

根据数字信号处理基于计算机的方法 4版（中文版）中关于理想单位冲击响应函数介绍（如下图）实现。

低通

// wc为截至频率，N为窗口长度
double *IdealLp(double wc, int N) //低通
{
   
   
    double M = (N - 1) / 2.0;
    double *hd = (double *)malloc(N * sizeof(double));
    for (int n = 0; n < N; n++)
    {
   
   
        if (n == M)
        {
   
   
            hd[n] = wc / PI;
        }
        else
        {
   
   
            hd[n] = sin(wc * (n - M)) / (PI * (n - M));
        }
    }

    return hd;
}

高通

// wc为截至频率，N为窗口长度
double *IdealHp(double wc, int N) //高通
{
   
   
    double M = (N - 1) / 2.0;
    double *hd = (double *)malloc(N * sizeof(double));
    for (int n = 0; n < N; n++)
    {
   
   
        if (n == M)
        {
   
   
            hd[n] = 1.0 - wc / PI;
        }
        else
        {
   
   
            hd[n] = -sin(wc * (n - M)) / (PI * (n - M));
        }
    }

    return hd;
}

带通

// wcl为低频截至频率，wch为高频截至频率，N为窗口长度
double *IdealBp(double wcl, double wch, int N) // 带通
{
   
   
    double M = (N - 1) / 2.0;
    double *hd = (double *)malloc(N * sizeof(double));

    hd = IdealLp(wcl, N);
    hd = subtract(IdealLp(wch, N), hd, N);

    return hd;
}

带通理想单位冲击响应函数可由低频和高频截至频率的低通理想单位冲击响应函数相减得到。

其中subtract为封装的double数组运算库中的减法函数，IdealLp为低通理想单位冲击响应函数。

带阻

// wcl为低频截至频率，wch为高频截至频率，N为窗口长度
double *IdealBs(double wcl, double wch, int N) // 带阻
{
   
   
    double M = (N - 1) / 2.0;
    double *hd = (double *)malloc(N * sizeof(double));
    for (int n = 0; n < N; n++)
    {
   
   
        if (n == M)
        {
   
   
            hd[n] = 1.0 - (wch - wcl) / PI;
        }
        else
        {
   
   
            hd[n] = (sin(wcl * (n - M)) - sin(wch * (n - M))) / (PI * (n - M));
        }
    }

    return hd;
}

至此，低通、高通、带通和带阻的理想单位冲击响应函数已全部实现。

窗函数

矩形窗

//全为1的数组
double *RectangleWindow(int n)
{
   
   
    double *result = (double *)malloc(n * sizeof(double));

    for (size_t i = 0; i < n; i++)
    {
   
   
        result[i] = 1.0;
    }

    return result;
}
</