C 语言-FIR 滤波器设计与应用

本文中将介绍在 C 语言如何部署 FIR 滤波器，并展示如何设计 FIR 滤波器来滤除心电信号 50Hz 工频…

需要包含的头文件

需要包含的头文件如下所示。

#include <stdio.h>
#include <math.h>
#include <malloc.h>
#include <stdlib.h>

实现 FIR 滤波器

FIR 滤波器的差分方程如下所示。

\(y(n)=\sum_{m=0}^{M}h(m)x(n-m)\)

其中 M 被称为 FIR 滤波器的阶数，\(h(m)\) 为 FIR 滤波器系数。

通过上述差分方程，可以很容易实现 FIR 滤波器，如下所示。

/*********************************************************************************************************
* 函数名称： FirFilter
* 函数功能： FIR 滤波器
* 输入参数： h ：双精度实型一维数组，长度为（m+1）。存放滤波器的系数。
*            px：双精度实型一维数组，长度为（m+1）。存放历史输入数据，初始状态必须为零。
*            m ：滤波器的阶数。
*            x ：新的采样值。
* 输出参数： void
* 返 回 值： 滤波后的采样值。
* 创建日期： 2023年10月09日
* 注    意：
*********************************************************************************************************/
double FirFilter(double* h, double* px, int m, double x)
{
  int i;
  double sum;

  //将新的采样数据存放到状态数组
  for (i = m; i >= 1; i--)
  {
    px[i] = px[i - 1];
  }
  px[0] = x;

  //根据 FIR 滤波器的定义求解滤波后的结果
  sum = 0;
  for (i = 0; i <= m; i++)
  {
    sum = sum + h[i] * px[i];
  }

  //返回滤波后的采样值
  return sum;
}

使用窗函数设计 FIR 数字滤波器

设 N-1 阶 FIR 数字滤波器的单位冲激响应为 h(n)，则传递函数 H(z) 为

\(H(z)=\sum_{n=0}^{N-1}h(n)z^{-n}\)

窗函数法的设计步骤如下：

1、根据给定的理想频率响应 \(H_{d}(e^{jw})\)，利用傅里叶反变换，求出单位冲激响应 \(h_{d}(n)\)

\(h_{d}(n)=\frac{1}{2\pi }\int_{-\pi }^{\pi }H_{d}(e^{jw})e^{jwn}dw\)

2、将 \(h_{d}(n)\) 乘以窗函数 \(w(n)\)，得到所要求的 FIR 滤波器的系数 \(h(n)\)

\(h(n)=w(n)h_{d}(n)\)

常用的窗函数有：

矩形窗

\(w(n)=1,0\leqslant n\leqslant (N-2)/10\)

图基（Tukey）窗

\(w(n)=0.5(1-cos\left ( \frac{10\pi n}{N+8}) \right ),0\leqslant n\leqslant (N-2)/10\)

\(w(n)=1,(N-2)/10\leqslant n\leqslant 9(N-2)/10\)

\(w(n)=0.5\left ( 1-cos\left ( \frac{10\pi (N-n-1)}{N+8} \right ) \right ),9(N-2)/10\leqslant n\leqslant N-1\)

三角窗

\(w(n)=1-\left | 1-\frac{2n}{N-1} \right |,0\leq n\leq N-1\)

汉宁（Hanning）窗

\(w(n)=0.5-0.5cos\frac{2\pi n}{N-1},0\leq n\leq N-1\)

海明（Hamming）窗

\(w(n)=0.54-0.46cos\frac{2\pi n}{N-1},0\leq n\leq N-1\)

布拉克曼（Blackman）窗

\(w(n)=0.42-0.5cos\frac{2\pi n}{N-1}+0.08cos\frac{4\pi n}{N-1},0\leq n\leq N-1\)

凯塞（Kaiser）窗

\(w(n)=\frac{I_{0}(\beta \sqrt{1-(1-2n/(N-1))^{2}})}{I_{0}(\beta )},0\leq n\leq N-1\)

其中 \(I_{0}(\beta )\) 是第一类零阶修正贝塞耳函数。\(\beta\) 是控制窗函数形状的参数，\(\beta\) 越大，\(w(n)\) 窗越窄，频谱的旁瓣越小，但主瓣也相应加宽。\(\beta\) 的典型值为 \(4\leqslant \beta \leqslant 9\)。\(\beta =0\) 时，凯塞窗变为矩形窗；\(\beta =5.44\) 时，凯塞窗与海明窗接近；\(\beta =8.5\) 时，凯塞窗与布拉克曼窗接近。

实现代码如下。

/*********************************************************************************************************
* 函数名称： Firwin
* 函数功能： FIR 数字滤波器的设计，窗函数方法
* 输入参数： n   ：整形变量。滤波器的阶数。
*            band：整形变量。滤波器的类型。取值为 1、2、3 和 4，分别对应低通、高通、带通和带阻滤波器。
*            fln ：双精度实型变量。
*            fhn ：双精度实型变量。
*                  对于低通和高通滤波器，fln：通带边界频率；
*                  对于带通和带阻滤波器，fln：带通下边界频率，fhn：带通上边界频率。
*            wn  ：整形变量。窗函数类型。取值为 1 到 7，分别对应矩形窗、图基窗、三角窗、汉宁窗、海明窗、
*                  布拉克曼窗和凯塞窗
*            h   ：双精度实型一维数组，长度为（n+1），存放 FIR 滤波器的系数。
* 输出参数： void
* 返 回 值： void
* 创建日期： 2023年10月03日
* 注    意：
*********************************************************************************************************/
void Firwin(int n, int band, double fln, double fhn, int wn, double* h)
{
  int i, n2, mid;
  double s, pi, wc1, wc2, beta, delay;
  double window(int type, int n, int i, double beta);
  beta = 0.0;
  if (wn == 7)
  {
    printf("input beta parameter of Kaiser window (3 < beta < 10)\n");
    scanf("%lf", &beta);
  }
  pi = 4.0 * atan(1.0);
  if ((n % 2) == 0)
  {
    n2 = n / 2 - 1;
    mid = 1;
  }
  else
  {
    n2 = n / 2;
    mid = 0;
  }
  delay = n / 2.0;
  wc1 = 2.0 * pi * fln;
  if (band >= 3)
  {
    wc2 = 2.0 * pi * fhn;
  }
  switch (band)
  {
    case 1:
    {
      for (i = 0; i <= n2; i++)
      {
        s = i - delay;
        h[i] = (sin(wc1 * s) / (pi * s)) * window(wn, n + 1, i, beta);
        h[n - i] = h[i];
      }
      if (mid == 1)
      {
        h[n / 2] = wc1 / pi;
      }
      break;
    }
    case 2:
    {
      for (i = 0; i <= n2; i++)
      {
        s = i - delay;
        h[i] = (sin(pi * s) - sin(wc1 * s)) / (pi * s);
        h[i] = h[i] * window(wn, n + 1, i, beta);
        h[n - i] = h[i];
      }
      if (mid == 1)
      {
        h[n / 2] = 1.0 - wc1 / pi;
      }
      break;
    }
    case 3:
    {
      for (i = 0; i <= n2; i++)
      {
        s = i - delay;
        h[i] = (sin(wc2 * s) - sin(wc1 * s)) / (pi * s);
        h[i] = h[i] * window(wn, n + 1, i, beta);
        h[n - i] = h[i];
      }
      if (mid == 1)
      {
        h[n / 2] = (wc2 - wc1) / pi;
      }
      break;
    }
    case 4:
    {
      for (i = 0; i <= n2; i++)
      {
        s = i - delay;
        h[i] = (sin(wc1 * s) + sin(pi * s) - sin(wc2 * s)) / (pi * s);
        h[i] = h[i] * window(wn, n + 1, i, beta);
        h[n - i] = h[i];
      }
      if (mid == 1)
      {
        h[n / 2] = (wc1 + pi - wc2) / pi;
      }
      break;
    }
  }
}
static double window(int type, int n, int i, double beta)
{
  int k;
  double pi, w;
  double keiser(int i, int n, double beta);
  pi = 4.0 * atan(1.0);
  w = 1.0;
  switch (type)
  {
    case 1: 
    {
      w = 1.0; 
      break;
    }
    case 2: 
    {
      k = (n - 2) / 10;
      if (i <= k)
      {
        w = 0.5 * (1.0 - cos(i * pi / (k + 1)));
      }
      if (i > n - k - 2)
      {
        w = 0.5 * (1.0 - cos((n - i - 1) * pi / (k + 1)));
      }
      break;
    }
    case 3:
    {
      w = 1.0 - fabs(1.0 - 2 * i / (n - 1.0));
      break;
    }
    case 4:
    {
      w = 0.5 * (1.0 - cos(2 * i * pi / (n - 1)));
      break;
    }
    case 5:
    {
      w = 0.54 - 0.46 * cos(2 * i * pi / (n - 1));
      break;
    }
    case 6:
    {
      w = 0.42 - 0.5 * cos(2 * pi / (n - 1)) + 0.08 * cos(4 * i * pi / (n - 1));
      break;
    }
    case 7:
    {
      w = keiser(i, n, beta);
      break;
    }
  }
  return (w);
}
static double keiser(int i, int n, double beta)
{
  double a, w, a2, b1, b2, beta1;
  double bessel0(double x);
  b1 = bessel0(beta);
  a = 2.0 * i / (double)(n - 1) - 1.0;
  a2 = a * a;
  beta1 = beta * sqrt(1.0 - a2);
  b2 = bessel0(beta1);
  w = b2 / b1;
  return (w);
}
static double bessel0(double x)
{
  int i;
  double d, y, d2, sum;
  y = x / 2.0;
  d = 1.0;
  sum = 1.0;
  for (i = 1; i <= 25; i++)
  {
    d = d * y / i;
    d2 = d * d;
    sum = sum + d2;
    if (d2 < sum * (1.0e-8))
    {
      break;
    }
  }
  return (sum);
}

下面给出主函数程序。它调用 Firwin 函数。通过人机对话，它可以设计低通、高通、带通和带阻这四种形式的滤波器。下面对输入参数进行说明。

band：滤波器的类型；

n ：滤波器的阶数；

fs ：采样频率；

对于低通和高通滤波器，fl：通带边界频率；

对于带通和带阻滤波器，f1：通带下边界频率，fh：通带上边界频率；

wn：窗函数的类型；

fname：幅频响应文件名。

/*********************************************************************************************************
* 函数名称： Gain
* 函数功能： 数字滤波器的频率响应
* 输入参数： b   ：双精度实型一维数组，长度为（m+1）。存放滤波器分子多项式的系数 b(i)。
*            a   ：双精度实型一维数组，长度为（n+1）。存放滤波器分母多项式的系数 a(i)。
*            m   ：整型变量。滤波器分子多项式的阶数。
*            n   ：整形变量。滤波器分母多项式的阶数。
*            x   ：双精度实型一维数组，长度为 len。
*                  当 sign = 0 时，存放滤波器频率响应的实部 Rc[H(w)]；
*                  当 sign = 1 时，存放滤波器幅频响应 |H(w)|；
*                  当 sign = 2 时，存放分贝表示的滤波器幅频响应 |H(w)|。
*            y   ：双精度实型一维数组，长度为 len。
*                  当 sign = 0 时，存放滤波器频率响应的虚部 Im[H(w)]；
*                  当 sign = 1 或 2 时，存放滤波器的相频响应 φ(w)。
*            sign：整形变量。
*                  当 sign = 0 时，计算滤波器频率响应的实部 Rc[H(w)] 和 虚部 Im[H(w)]；
*                  当 sign = 1 时，计算滤波器的幅频响应 |H(w)| 和相频响应 φ(w)；
*                  当 sign = 2 时，计算滤波器的幅频响应 |H(w)| （用 dB 表示）和相频响应 φ(w)。
* 输出参数： void
* 返 回 值： void
* 创建日期： 2023年09月29日
* 注    意：
*********************************************************************************************************/
void Gain(double* b, double* a, int m, int n, double* x, double* y, int len, int sign)
{
  int i, k;
  double ar, ai, br, bi, zr, zi, im, re, den, numr, numi, freq, temp;
  for (k = 0; k < len; k++)
  {
    freq = k * 0.5 / (len - 1);
    zr = cos(-8.0 * atan(1.0) * freq);
    zi = sin(-8.0 * atan(1.0) * freq);
    br = 0.0;
    bi = 0.0;
    for (i = m; i > 0; i--)
    {
      re = br;
      im = bi;
      br = (re + b[i]) * zr - im * zi;
      bi = (re + b[i]) * zi + im * zr;
    }
    ar = 0.0;
    ai = 0.0;
    for (i = n; i > 0; i--)
    {
      re = ar;
      im = ai;
      ar = (re + a[i]) * zr - im * zi;
      ai = (re + a[i]) * zi + im * zr;
    }
    br = br + b[0];
    ar = ar + 1.0;
    numr = ar * br + ai * bi;
    numi = ar * bi - ai * br;
    den = ar * ar + ai * ai;
    x[k] = numr / den;
    y[k] = numi / den;
    switch (sign)
    {
      case 1:
      {
        temp = sqrt(x[k] * x[k] + y[k] * y[k]);
        y[k] = atan2(y[k], x[k]);
        x[k] = temp;
        break;
      }
      case 2:
      {
        temp = x[k] * x[k] + y[k] * y[k];
        y[k] = atan2(y[k], x[k]);
        x[k] = 10.0 * log10(temp);
      }
    }
  }
}

/*********************************************************************************************************
* 函数名称： main
* 函数功能： 主函数
* 输入参数： void
* 输出参数： void
* 返 回 值： int
* 创建日期： 2023年10月03日
* 注    意：
*********************************************************************************************************/
int main(void)
{
  int i, j, n, n2, band, wn;
  double fl, fh, fs, freq;
  static double h[10240], c[10240], x[300], y[300];
  char fname[40];
  FILE* fp;
  c[1] = 0.0;
  printf("Select one of the four types for FIR digital filter\n");
  printf("1 -- lowpass; 2 -- highpass\n");
  printf("3 -- bandpass; 4 -- bandstop\n");
  scanf("%d", &band);
  printf("Input the filter order\n");
  scanf("%d", &n);
  printf("Input low cutoff frequrncy fl\n");
  scanf("%lf", &fl);
  fh = 0.0;
  if (band >= 3)
  {
    printf("Input high cutoff frequrncy fh\n");
    scanf("%lf", &fh);
  }
  printf("Input sample frequence fs\n");
  scanf("%lf", &fs);
  printf("Select window\n");
  printf("1 -- rectangular; 2 -- tapered rectangular\n");
  printf("3 -- triangular; 4 -- Hanning\n");
  printf("5 -- Hamming; 6 -- Blackman\n");
  printf("7 -- Kaiser\n");
  scanf("%d", &wn);
  fl = fl / fs;
  fh = fh / fs;
  Firwin(n, band, fl, fh, wn, h);
  printf("FIR digital filter\n");
  printf(" * * * * impulse respomse * * * * \n\n");
  n2 = n / 2;
  for (i = 0; i <= n2; i++)
  {
    j = n - i;
    printf("h(%2d) = %12.8lf = h(%2d)\n", i, h[i], j);
  }
  printf("\n * * * * impulse respomse by C * * * * \n\n");
  printf("static const double s_arrFirH[%d] = {\n", n + 1);
  for (i = 0; i <= n; i++)
  {
    printf("%12.8lf, ", h[i]);
    if (0 == ((i + 1) % 10))
    {
      printf("\n");
    }
  }
  printf("}\n");
  
  printf("\nInput file name of frequency response\n");
  scanf("%s", fname);
  fp = fopen(fname, "w");
  if (fp == NULL)
  {
    printf("Cannot open this file\n");
    exit(1);
  }
  Gain(h, c, n, 1, x, y, 300, 2);
  for (i = 0; i < 300; i++)
  {
    freq = 0.5 * i / 299;
    fprintf(fp, "%lf,%lf\n", freq, x[i]);
  }
  fclose(fp);
  return 0;
}

下面设计一个 128 阶的 FIR 低通滤波器，采样率 2kHz，通带边界频率为 25Hz，选用汉宁窗。选择参数 n = 128，band = 1，fl = 25， fs = 2000，wn = 4。

终端的输入输出如下所示。

Select one of the four types for FIR digital filter
1 -- lowpass; 2 -- highpass
3 -- bandpass; 4 -- bandstop
1
Input the filter order
128
Input low cutoff frequrncy fl
25
Input sample frequence fs
2000
Select window
1 -- rectangular; 2 -- tapered rectangular
3 -- triangular; 4 -- Hanning
5 -- Hamming; 6 -- Blackman
7 -- Kaiser
4
FIR digital filter
 * * * * impulse respomse * * * *

h( 0) =  -0.00000000 = h(128)
h( 1) =  -0.00000296 = h(127)
h( 2) =  -0.00001221 = h(126)
h( 3) =  -0.00002815 = h(125)
h( 4) =  -0.00005097 = h(124)
h( 5) =  -0.00008059 = h(123)
h( 6) =  -0.00011670 = h(122)
h( 7) =  -0.00015871 = h(121)
h( 8) =  -0.00020575 = h(120)
h( 9) =  -0.00025668 = h(119)
h(10) =  -0.00031008 = h(118)
h(11) =  -0.00036427 = h(117)
h(12) =  -0.00041730 = h(116)
h(13) =  -0.00046699 = h(115)
h(14) =  -0.00051091 = h(114)
h(15) =  -0.00054645 = h(113)
h(16) =  -0.00057083 = h(112)
h(17) =  -0.00058112 = h(111)
h(18) =  -0.00057428 = h(110)
h(19) =  -0.00054721 = h(109)
h(20) =  -0.00049677 = h(108)
h(21) =  -0.00041984 = h(107)
h(22) =  -0.00031335 = h(106)
h(23) =  -0.00017435 = h(105)
h(24) =   0.00000000 = h(104)
h(25) =   0.00021232 = h(103)
h(26) =   0.00046500 = h(102)
h(27) =   0.00076017 = h(101)
h(28) =   0.00109963 = h(100)
h(29) =   0.00148483 = h(99)
h(30) =   0.00191684 = h(98)
h(31) =   0.00239630 = h(97)
h(32) =   0.00292340 = h(96)
h(33) =   0.00349789 = h(95)
h(34) =   0.00411901 = h(94)
h(35) =   0.00478552 = h(93)
h(36) =   0.00549567 = h(92)
h(37) =   0.00624721 = h(91)
h(38) =   0.00703742 = h(90)
h(39) =   0.00786305 = h(89)
h(40) =   0.00872043 = h(88)
h(41) =   0.00960543 = h(87)
h(42) =   0.01051350 = h(86)
h(43) =   0.01143972 = h(85)
h(44) =   0.01237883 = h(84)
h(45) =   0.01332529 = h(83)
h(46) =   0.01427329 = h(82)
h(47) =   0.01521684 = h(81)
h(48) =   0.01614980 = h(80)
h(49) =   0.01706596 = h(79)
h(50) =   0.01795909 = h(78)
h(51) =   0.01882298 = h(77)
h(52) =   0.01965152 = h(76)
h(53) =   0.02043879 = h(75)
h(54) =   0.02117906 = h(74)
h(55) =   0.02186687 = h(73)
h(56) =   0.02249711 = h(72)
h(57) =   0.02306505 = h(71)
h(58) =   0.02356640 = h(70)
h(59) =   0.02399733 = h(69)
h(60) =   0.02435454 = h(68)
h(61) =   0.02463528 = h(67)
h(62) =   0.02483737 = h(66)
h(63) =   0.02495926 = h(65)
h(64) =   0.02500000 = h(64)

 * * * * impulse respomse by C * * * *

static const double s_arrFirH[129] = {
 -0.00000000,  -0.00000296,  -0.00001221,  -0.00002815,  -0.00005097,  -0.00008059,  -0.00011670,  -0.00015871,  -0.00020575,  -0.00025668,
 -0.00031008,  -0.00036427,  -0.00041730,  -0.00046699,  -0.00051091,  -0.00054645,  -0.00057083,  -0.00058112,  -0.00057428,  -0.00054721,
 -0.00049677,  -0.00041984,  -0.00031335,  -0.00017435,   0.00000000,   0.00021232,   0.00046500,   0.00076017,   0.00109963,   0.00148483,
  0.00191684,   0.00239630,   0.00292340,   0.00349789,   0.00411901,   0.00478552,   0.00549567,   0.00624721,   0.00703742,   0.00786305,
  0.00872043,   0.00960543,   0.01051350,   0.01143972,   0.01237883,   0.01332529,   0.01427329,   0.01521684,   0.01614980,   0.01706596,
  0.01795909,   0.01882298,   0.01965152,   0.02043879,   0.02117906,   0.02186687,   0.02249711,   0.02306505,   0.02356640,   0.02399733,
  0.02435454,   0.02463528,   0.02483737,   0.02495926,   0.02500000,   0.02495926,   0.02483737,   0.02463528,   0.02435454,   0.02399733,
  0.02356640,   0.02306505,   0.02249711,   0.02186687,   0.02117906,   0.02043879,   0.01965152,   0.01882298,   0.01795909,   0.01706596,
  0.01614980,   0.01521684,   0.01427329,   0.01332529,   0.01237883,   0.01143972,   0.01051350,   0.00960543,   0.00872043,   0.00786305,
  0.00703742,   0.00624721,   0.00549567,   0.00478552,   0.00411901,   0.00349789,   0.00292340,   0.00239630,   0.00191684,   0.00148483,
  0.00109963,   0.00076017,   0.00046500,   0.00021232,   0.00000000,  -0.00017435,  -0.00031335,  -0.00041984,  -0.00049677,  -0.00054721,
 -0.00057428,  -0.00058112,  -0.00057083,  -0.00054645,  -0.00051091,  -0.00046699,  -0.00041730,  -0.00036427,  -0.00031008,  -0.00025668,
 -0.00020575,  -0.00015871,  -0.00011670,  -0.00008059,  -0.00005097,  -0.00002815,  -0.00001221,  -0.00000296,  -0.00000000, }

Input file name of frequency response
result.csv

输出的幅频特性曲线如下所示。

本文提供的心电信号是纯净的，不带 50Hz 工频干扰的。我们需要先手动添加 50Hz 工频干扰，然后再对心电信号进行滤波。具体如下所示。FIR 滤波器的参数可以直接从终端里复制，不用手动输入。

int main(void)
{
  extern const double g_arrEcgWave[4000];
  static double s_arrInput[4000] = {0};
  static double s_arrOutput[4000] = {0};
  static const double s_arrFirH[129] = {
   -0.00000000,  -0.00000296,  -0.00001221,  -0.00002815,  -0.00005097,  -0.00008059,  -0.00011670,  -0.00015871,  -0.00020575,  -0.00025668,
   -0.00031008,  -0.00036427,  -0.00041730,  -0.00046699,  -0.00051091,  -0.00054645,  -0.00057083,  -0.00058112,  -0.00057428,  -0.00054721,
   -0.00049677,  -0.00041984,  -0.00031335,  -0.00017435,   0.00000000,   0.00021232,   0.00046500,   0.00076017,   0.00109963,   0.00148483,
    0.00191684,   0.00239630,   0.00292340,   0.00349789,   0.00411901,   0.00478552,   0.00549567,   0.00624721,   0.00703742,   0.00786305,
    0.00872043,   0.00960543,   0.01051350,   0.01143972,   0.01237883,   0.01332529,   0.01427329,   0.01521684,   0.01614980,   0.01706596,
    0.01795909,   0.01882298,   0.01965152,   0.02043879,   0.02117906,   0.02186687,   0.02249711,   0.02306505,   0.02356640,   0.02399733,
    0.02435454,   0.02463528,   0.02483737,   0.02495926,   0.02500000,   0.02495926,   0.02483737,   0.02463528,   0.02435454,   0.02399733,
    0.02356640,   0.02306505,   0.02249711,   0.02186687,   0.02117906,   0.02043879,   0.01965152,   0.01882298,   0.01795909,   0.01706596,
    0.01614980,   0.01521684,   0.01427329,   0.01332529,   0.01237883,   0.01143972,   0.01051350,   0.00960543,   0.00872043,   0.00786305,
    0.00703742,   0.00624721,   0.00549567,   0.00478552,   0.00411901,   0.00349789,   0.00292340,   0.00239630,   0.00191684,   0.00148483,
    0.00109963,   0.00076017,   0.00046500,   0.00021232,   0.00000000,  -0.00017435,  -0.00031335,  -0.00041984,  -0.00049677,  -0.00054721,
   -0.00057428,  -0.00058112,  -0.00057083,  -0.00054645,  -0.00051091,  -0.00046699,  -0.00041730,  -0.00036427,  -0.00031008,  -0.00025668,
   -0.00020575,  -0.00015871,  -0.00011670,  -0.00008059,  -0.00005097,  -0.00002815,  -0.00001221,  -0.00000296,  -0.00000000, };
  static double s_arrPx[sizeof(s_arrFirH) / sizeof(double)] = {0};
  int i, len;
  long seed;
  double pi, time;
  FILE* fp;

  //定义数据量
  len = 4000;

  //获取输入数据，加 50Hz 工频干扰
  seed = 13579l;
  pi = 3.1415926535;
  time = 0;
  for (i = 0; i < len; i++)
  {
    s_arrInput[i] = g_arrEcgWave[i] + 0.1 * sin(time * 2 * pi / (1.0 / 50.0));
    time = time + 0.0005;
  }

  //FIR 滤波
  for (i = 0; i < len; i++)
  {
    s_arrOutput[i] = FirFilter(s_arrFirH, s_arrPx, (sizeof(s_arrFirH) / sizeof(double)) - 1, s_arrInput[i]);
  }

  //输出波形数据，按照原始信号、参考信号、滤波后信号顺序
  fp = fopen("FirFilter.csv", "w");
  for (i = 0; i < len; i++)
  {
    fprintf(fp, "%d,%10.7lf,%10.7lf,%10.7lf\n", i, s_arrInput[i], g_arrEcgWave[i], s_arrOutput[i]);
  }
  return 0;
}

最终实验结果如下所示。蓝色为参考信号，绿色为输入信号，黄色为输出信号。

使用频域最小误差平方设计 FIR 数字滤波器

设 N-1 阶 FIR 数字滤波器的单位冲激响应为 h(n)，则传递函数 H(z) 为

\(H(z)=\sum_{n=0}^{N-1}h(n)z^{-n}\)

FIR 滤波器的实际幅度相应 \(H(w)\) 与理想幅度响应 \(H_{d}(w)\) 的误差函数定义为

\(E=\frac{1}{2\pi }\int_{-\pi }^{\pi }\left | H(w)-H_{d}(w) \right |^{2}dw\)

根据帕塞瓦定理，误差函数可以表示为

\(E=\sum_{n=-\infty }^{+\infty }\left | h(n)-h_{d}(n) \right |^{2}=\sum_{n=-M}^{M}\left | h(n)-h_{d}(n) \right |^{2}+\sum_{n=M+1}^{\infty }2h_{d}^{2}(n)\)

其中 \(M=(N-1)/2\)。该式表明，若要使 \(E\) 最小，必须选择 \(h(n)\)，使它与 \(h_{d}(n)\) 对应的 N 个值相等，即 \(h(n)=h_{d}(n),-M\leqslant n\leq M\)。

为了减少吉伯斯（Gibbs）效应，通常在通带和阻带之间加一个过渡带，\(f_{c}\leqslant f\leqslant f_{s}\)。过渡带函数常采用 \(P\) 阶样条函数，此时滤波器的系数为

\(h(n)=\left ( \frac{sin[\pi (f_{s}-f_{c})(n-M)/P]}{\pi (f_{s}-f_{c})(n-M)/P} \right )^{p}\frac{sin[\pi (f_{s}+f_{c})(n-M)]}{\pi (n-M)},0\leqslant n\leqslant N-1\)

过渡带函数也可采用升余弦函数，此时滤波器的系数为

\(h(n)=\frac{cos[\pi (f_{s}-f_{c})(n-M)]}{1-4(f_{s}-f_{c})^{2}(n-M)^{2}}\frac{sin[\pi (f_{s}+f_{c})(n-M)]}{\pi (n-M)},0\leqslant n\leq N-1\)

实现代码如下所示。

/*********************************************************************************************************
* 函数名称： Firls
* 函数功能： FIR 数字滤波器的设计。用频域最小误差平方方法设计线性相位 FIR 低通数字滤波器
* 输入参数： n ：整形变量。滤波器的阶数。
*            fc：双精度实型变量。通带边界频率。
*            fs：双精度实型变量。阻带边界频率。
*            tp：整形变量。过渡带函数的类型。tp = 0，表示升余弦函数；tp >= 1，表示 tp 阶样条函数。
*            h ：双精度实型一维数组，长度为（n+1），存放 FIR 滤波器的系数。
* 输出参数： void
* 返 回 值： void
* 创建日期： 2023年10月03日
* 注    意：
*********************************************************************************************************/
void Firls(int n, double fc, double fs, int tp, double* h)
{
  double fq, fr;
  void ls(double* h, int n, double fc);
  void wgt(double* h, int n, int tp, double fq);
  fq = fs - fc;
  fr = fs + fc;
  ls(h, n, fr);
  wgt(h, n, tp, fq);
}
static void ls(double* h, int n, double fc)
{
  int i, m, n2;
  double q, am, pi;
  pi = 4.0 * atan(1.0);
  m = n / 2;
  am = n / 2.0;
  n2 = (n - 1) / 2;
  if (m == am)
  {
    h[m] = fc;
  }
  for (i = 0; i <= n2; i++)
  {
    q = pi * (i - am);
    h[i] = sin(fc * q) / q;
  }
}
static void wgt(double* h, int n, int tp, double fq)
{
  int i;
  double q, am, pi, ql, wt;
  pi = 4.0 * atan(1.0);
  q = pi * fq;
  am = n / 2.0;
  if (fq == 0.0)
  {
    return;
  }
  if (tp != 0)
  {
    for (i = 0; i < am; i++)
    {
      ql = q * (i - am) / tp;
      wt = pow(sin(ql) / ql, tp);
      h[i] = wt * h[i];
    }
  }
  else
  {
    for (i = 0; i < am; i++)
    {
      wt = cos(q * (i - am));
      if (fabs(wt) > 1.0e-6)
      {
        wt = wt / (1 - pow((2 * fq * (i - am)), 2));
      }
      else
      {
        wt = pi / 4.0;
      }
      h[i] = wt * h[i];
    }
  }
  for (i = 0; i < am; i++)
  {
    h[n - i] = h[i];
  }
}

下面给出主函数程序，它调用 Firls 函数。通过人机对话，它可以设计线性相位 FIR 低通数字滤波器。下面对输入参数进行说明。

n：滤波器的阶数；

f：采样频率；

fc：通带边界频率；

fs：阻带边界频率；

tp：样条函数的阶数；

fname：幅频相应文件名。

int main(void)
{
  int i, j, n, n2, tp;
  double f, fc, fs, freq;
  static double h[10240], c[10240], x[300], y[300];
  char fname[40];
  FILE* fp;
  c[1] = 0.0;
  printf("Input this filter order\n");
  scanf("%d", &n);
  printf("Input passband cutoff frequrncy fc\n");
  scanf("%lf", &fc);
  printf("Input stopband edge frequrncy fs\n");
  scanf("%lf", &fs);
  printf("Input sample frequrncy f\n");
  scanf("%lf", &f);
  fc = fc / f;
  fs = fs / f;
  printf("Input the order of spline n\n");
  scanf("%d", &tp);
  Firls(n, fc, fs, tp, h);
  printf("FIR digital filter\n");
  printf(" * * * * impulse respomse * * * * \n\n");
  n2 = n / 2;
  for (i = 0; i <= n2; i++)
  {
    j = n - i;
    printf("h(%2d) = %12.8lf = h(%2d)\n", i, h[i], j);
  }
  printf("\n * * * * impulse respomse by C * * * * \n\n");
  printf("static const double s_arrFirH[%d] = {\n", n + 1);
  for (i = 0; i <= n; i++)
  {
    printf("%12.8lf, ", h[i]);
    if (0 == ((i + 1) % 10))
    {
      printf("\n");
    }
  }
  printf("}\n");
  printf("\nInput file name of frequency response\n");
  scanf("%s", fname);
  fp = fopen(fname, "w");
  if (fp == NULL)
  {
    printf("Cannot open this file\n");
    exit(1);
  }
  Gain(h, c, n, 1, x, y, 300, 2);
  for (i = 0; i < 300; i++)
  {
    freq = 0.5 * i / 299;
    fprintf(fp, "%lf,%lf\n", freq, x[i]);
  }
  fclose(fp);
  return 0;
}

设计一个 128 阶的 FIR 低通滤波器，其采样率为 2kHz，通带边界频率为 20Hz，阻带边界频率为 25Hz，采用一阶样条函数作为过渡带。选择参数 n = 128，fc = 20，fs = 25，f = 2000，tp = 1。

终端输入输出如下所示。

Input this filter order
128
Input passband cutoff frequrncy fc
20
Input stopband edge frequrncy fs
25
Input sample frequrncy f
2000
Input the order of spline n
1
FIR digital filter
 * * * * impulse respomse * * * *

h( 0) =  -0.00468235 = h(128)
h( 1) =  -0.00468693 = h(127)
h( 2) =  -0.00466754 = h(126)
h( 3) =  -0.00462343 = h(125)
h( 4) =  -0.00455392 = h(124)
h( 5) =  -0.00445842 = h(123)
h( 6) =  -0.00433640 = h(122)
h( 7) =  -0.00418745 = h(121)
h( 8) =  -0.00401123 = h(120)
h( 9) =  -0.00380749 = h(119)
h(10) =  -0.00357607 = h(118)
h(11) =  -0.00331693 = h(117)
h(12) =  -0.00303011 = h(116)
h(13) =  -0.00271576 = h(115)
h(14) =  -0.00237410 = h(114)
h(15) =  -0.00200550 = h(113)
h(16) =  -0.00161039 = h(112)
h(17) =  -0.00118932 = h(111)
h(18) =  -0.00074293 = h(110)
h(19) =  -0.00027196 = h(109)
h(20) =   0.00022274 = h(108)
h(21) =   0.00074024 = h(107)
h(22) =   0.00127951 = h(106)
h(23) =   0.00183945 = h(105)
h(24) =   0.00241883 = h(104)
h(25) =   0.00301637 = h(103)
h(26) =   0.00363072 = h(102)
h(27) =   0.00426042 = h(101)
h(28) =   0.00490395 = h(100)
h(29) =   0.00555975 = h(99)
h(30) =   0.00622616 = h(98)
h(31) =   0.00690150 = h(97)
h(32) =   0.00758400 = h(96)
h(33) =   0.00827189 = h(95)
h(34) =   0.00896332 = h(94)
h(35) =   0.00965643 = h(93)
h(36) =   0.01034934 = h(92)
h(37) =   0.01104013 = h(91)
h(38) =   0.01172688 = h(90)
h(39) =   0.01240766 = h(89)
h(40) =   0.01308055 = h(88)
h(41) =   0.01374362 = h(87)
h(42) =   0.01439497 = h(86)
h(43) =   0.01503271 = h(85)
h(44) =   0.01565498 = h(84)
h(45) =   0.01625997 = h(83)
h(46) =   0.01684589 = h(82)
h(47) =   0.01741100 = h(81)
h(48) =   0.01795362 = h(80)
h(49) =   0.01847214 = h(79)
h(50) =   0.01896500 = h(78)
h(51) =   0.01943071 = h(77)
h(52) =   0.01986787 = h(76)
h(53) =   0.02027516 = h(75)
h(54) =   0.02065133 = h(74)
h(55) =   0.02099524 = h(73)
h(56) =   0.02130585 = h(72)
h(57) =   0.02158219 = h(71)
h(58) =   0.02182344 = h(70)
h(59) =   0.02202883 = h(69)
h(60) =   0.02219775 = h(68)
h(61) =   0.02232968 = h(67)
h(62) =   0.02242421 = h(66)
h(63) =   0.02248104 = h(65)
h(64) =   0.02250000 = h(64)

 * * * * impulse respomse by C * * * *

static const double s_arrFirH[129] = {
 -0.00468235,  -0.00468693,  -0.00466754,  -0.00462343,  -0.00455392,  -0.00445842,  -0.00433640,  -0.00418745,  -0.00401123,  -0.00380749,
 -0.00357607,  -0.00331693,  -0.00303011,  -0.00271576,  -0.00237410,  -0.00200550,  -0.00161039,  -0.00118932,  -0.00074293,  -0.00027196,
  0.00022274,   0.00074024,   0.00127951,   0.00183945,   0.00241883,   0.00301637,   0.00363072,   0.00426042,   0.00490395,   0.00555975,
  0.00622616,   0.00690150,   0.00758400,   0.00827189,   0.00896332,   0.00965643,   0.01034934,   0.01104013,   0.01172688,   0.01240766,
  0.01308055,   0.01374362,   0.01439497,   0.01503271,   0.01565498,   0.01625997,   0.01684589,   0.01741100,   0.01795362,   0.01847214,
  0.01896500,   0.01943071,   0.01986787,   0.02027516,   0.02065133,   0.02099524,   0.02130585,   0.02158219,   0.02182344,   0.02202883,
  0.02219775,   0.02232968,   0.02242421,   0.02248104,   0.02250000,   0.02248104,   0.02242421,   0.02232968,   0.02219775,   0.02202883,
  0.02182344,   0.02158219,   0.02130585,   0.02099524,   0.02065133,   0.02027516,   0.01986787,   0.01943071,   0.01896500,   0.01847214,
  0.01795362,   0.01741100,   0.01684589,   0.01625997,   0.01565498,   0.01503271,   0.01439497,   0.01374362,   0.01308055,   0.01240766,
  0.01172688,   0.01104013,   0.01034934,   0.00965643,   0.00896332,   0.00827189,   0.00758400,   0.00690150,   0.00622616,   0.00555975,
  0.00490395,   0.00426042,   0.00363072,   0.00301637,   0.00241883,   0.00183945,   0.00127951,   0.00074024,   0.00022274,  -0.00027196,
 -0.00074293,  -0.00118932,  -0.00161039,  -0.00200550,  -0.00237410,  -0.00271576,  -0.00303011,  -0.00331693,  -0.00357607,  -0.00380749,
 -0.00401123,  -0.00418745,  -0.00433640,  -0.00445842,  -0.00455392,  -0.00462343,  -0.00466754,  -0.00468693,  -0.00468235, }

Input file name of frequency response
result.csv