【安富莱DSP教程】第37章 FIR滤波器的实现

席萌0209

特别说明：完整45期数字信号处理教程，原创高性能示波器代码全开源地址：链接

第37章 FIR滤波器的实现

    本章节讲解FIR滤波器的低通，高通，带通和带阻滤波器的实现。
    37.1 FIR滤波器介绍
    37.2 Matlab工具箱生成C头文件
    37.3 FIR低通滤波器设计
    37.4 FIR高通滤波器设计
    37.5 FIR带通滤波器设计
    37.6 FIR带阻滤波器设计
    37.7 切比雪夫窗口设计带通滤波器
    37.8 FIR滤波后的群延迟
    37.9 总结

37.1 FIR滤波器介绍

    ARM官方提供的FIR库支持Q7，Q15，Q31和浮点四种数据类型。其中Q15和Q31提供了快速算法版本。
    FIR滤波器的基本算法是一种乘法-累加（MAC）运行，输出表达式如下：

y[n] = b[0] * x[n] + b[1] * x[n-1] + b[2] * x[n-2] + ...+ b[numTaps-1] * x[n-numTaps+1]

结构图如下：

这种网络结构就是在35.2.1小节所讲的直接型结构。

席萌0209

37.2 Matlab工具箱fdatool生成C头文件


    下面我们讲解下如何通过fdatool工具生成C头文件，也就是生成滤波器系数。首先在matlab的命
窗口输入fadtool就能打开这个工具箱：

fadtool界面打开效果如下：

FIR滤波器的低通，高通，带通，带阻滤波的设置会在下面一一讲解，这里说一下设置后相应参数后如何生成滤波器系数。参数设置好以后点击如下按钮：

点击Design Filter按钮以后就生成了所需的滤波器系数，生成滤波器系数以后点击fadtool界面上的菜单Targets->Generate C header ,打开后显示如下界面：

然后点击Generate，生成如下界面：

再点击保存，并打开fdatool.h文件，可以看到生成的系数：

1. /*
2. * Filter Coefficients (C Source) generated by the Filter Design and Analysis Tool
3. *
4. * Generated by MATLAB(R) 7.14 and the Signal Processing Toolbox 6.17.
5. *
6. * Generated on: 22-Dec-2014 21:34:29
7. *
8. */
10. /*
11. * Discrete-Time FIR Filter (real)
12. * -------------------------------
13. * Filter Structure  : Direct-Form FIR
14. * Filter Length     : 51
15. * Stable            : Yes
16. * Linear Phase      : Yes (Type 1)
17. */
19. /* General type conversion for MATLAB generated C-code  */
20. #include "tmwtypes.h"
21. /*
22. * Expected path to tmwtypes.h
23. * C:Program FilesMATLABR2012aexternincludetmwtypes.h
24. */
25. /*
26. * Warning - Filter coefficients were truncated to fit specified data type.  
27. *   The resulting response may not match generated theoretical response.
28. *   Use the Filter Design & Analysis Tool to design accurate
29. *   single-precision filter coefficients.
30. */
31. const int BL = 51;
32. const real32_T B[51] = {
33.   -0.0009190982091, -0.00271769613,-0.002486952813, 0.003661438357,   0.0136509249,
34.     0.01735116541,  0.00766530633,-0.006554719061,-0.007696784101, 0.006105459295,
35.     0.01387391612,0.0003508617228, -0.01690892503,-0.008905642666,  0.01744112931,
36.     0.02074504457,  -0.0122964941, -0.03424086422,-0.001034529647,  0.04779030383,
37.     0.02736303769, -0.05937951803, -0.08230702579,  0.06718690693,   0.3100151718,
38.      0.4300478697,   0.3100151718,  0.06718690693, -0.08230702579, -0.05937951803,
39.     0.02736303769,  0.04779030383,-0.001034529647, -0.03424086422,  -0.0122964941,
40.     0.02074504457,  0.01744112931,-0.008905642666, -0.01690892503,0.0003508617228,
41.     0.01387391612, 0.006105459295,-0.007696784101,-0.006554719061,  0.00766530633,
42.     0.01735116541,   0.0136509249, 0.003661438357,-0.002486952813, -0.00271769613,
43.   -0.0009190982091
44. };

复制代码

上面数组B[51]中的数据就是滤波器系数。下面小节讲解如何使用fdatool配置FIR低通，高通，带通和带阻滤波。关于fdatool的其它用法，
大家可以在matlab命令窗口中输入help fadtool打开帮助文档进行学习。

席萌0209

37.3 FIR低通滤波器设计


    本章使用的FIR滤波器函数是arm_fir_f32。下面使用此函数设计FIR低通，高通，带通和带阻滤波器。

37.3.1 函数arm_fir_f32说明


函数定义如下：
    void arm_fir_f32(
      const arm_fir_instance_f32 * S,
      float32_t * pSrc,
      float32_t * pDst,
      uint32_t blockSize)
参数定义：
      [in]   *S       points to an instance of the floating-point FIR filter structure.  
       [in]  *pSrc    points to the block of input data.  
     [out] *pDst    points to the block of output data.  
     [in]  blockSize number of samples to process per call.  
    return     none.  
注意事项：
结构arm_fir_instance_f32的定义如下（在文件arm_math.h文件）：
      typedef struct
      {
          uint16_t numTaps;     /**< number of filter coefficients in the filter. */
      float32_t *pState;      /**< points to the state variable array. The array is of length numTaps+blockSize-1. */
          float32_t *pCoeffs;    /**< points to the coefficient array. The array is of length numTaps. */
       } arm_fir_instance_f32;
1. 参数pCoeffs指向滤波因数，滤波因数数组长度为numTaps。但要注意pCoeffs指向的滤波因数应该按照如下的逆序进行排列：
      {b[numTaps-1], b[numTaps-2], b[N-2], ..., b[1], b[0]}  
    但满足线性相位特性的FIR滤波器具有奇对称或者偶对称的系数，偶对称时逆序排列还是他本身。
2. pState指向状态变量数组，这个数组用于函数内部计算数据的缓存。
3. blockSize 这个参数的大小没有特殊要求，用户只需保证大于1小于等于采样点个数即可。

37.3.2 fdatool获取低通滤波器系数


    设计一个如下的例子：
    信号由50Hz正弦波和200Hz正弦波组成，采样率1Kbps，现设计一个低通滤波器，截止频率125Hz，采样320个数据，采用函数fir1进行设计（注意这个函数是基于窗口的方法设计FIR滤波，默认是hamming窗），滤波器阶数设置为28。fadtool的配置如下：

配置好低通滤波器后，具体滤波器系数的生成大家参考本章第二小节的方法即可。

37.3.3 低通滤波器实现


    通过工具箱fdatool获得低通滤波器系数后在开发板上运行函数arm_fir_f32 来测试低通滤波器的效果。

1. #define TEST_LENGTH_SAMPLES  320    /* 采样点数 */
2. #define BLOCK_SIZE           32     /* 调用一次arm_fir_f32处理的采样点个数 */
3. #define NUM_TAPS             29     /* 滤波器系数个数 */
5. uint32_t blockSize = BLOCK_SIZE;
6. uint32_t numBlocks = TEST_LENGTH_SAMPLES/BLOCK_SIZE;            /* 需要调用arm_fir_f32的次数 */
8. static float32_t testInput_f32_50Hz_200Hz[TEST_LENGTH_SAMPLES]; /* 采样点 */
9. static float32_t testOutput[TEST_LENGTH_SAMPLES];               /* 滤波后的输出 */
10. static float32_t firStateF32[BLOCK_SIZE + NUM_TAPS - 1];        /* 状态缓存，大小numTaps + blockSize - 1*/
12. /* 低通滤波器系数 通过fadtool获取*/
13. const float32_t firCoeffs32LP[NUM_TAPS] = {
14.   -0.001822523074f,  -0.001587929321f,  1.226008847e-18f,  0.003697750857f,  0.008075430058f,
15.   0.008530221879f,   -4.273456581e-18f, -0.01739769801f,   -0.03414586186f,  -0.03335915506f,
16.   8.073562366e-18f,  0.06763084233f,    0.1522061825f,     0.2229246944f,    0.2504960895f,
17.   0.2229246944f,     0.1522061825f,     0.06763084233f,    8.073562366e-18f, -0.03335915506f,
18.   -0.03414586186f,   -0.01739769801f,   -4.273456581e-18f, 0.008530221879f,  0.008075430058f,
19.   0.003697750857f,   1.226008847e-18f,  -0.001587929321f,  -0.001822523074f
20. };
21. /*
22. *********************************************************************************************************
23. *    函 数 名: arm_fir_f32_lp
24. *    功能说明: 调用函数arm_fir_f32_lp实现低通滤波器
25. *    形    参：无
26. *    返 回 值: 无
27. *********************************************************************************************************
28. */
29. static void arm_fir_f32_lp(void)
30. {
31. uint32_t i;
32. arm_fir_instance_f32 S;
33. float32_t  *inputF32, *outputF32;
35. /* 初始化输入输出缓存指针 */
36. inputF32 = &testInput_f32_50Hz_200Hz[0];
37. outputF32 = &testOutput[0];
39. /* 初始化结构体S */
40. arm_fir_init_f32(&S, NUM_TAPS, (float32_t *)&firCoeffs32LP[0], &firStateF32[0], blockSize);
42. /* 实现FIR滤波 */
43. for(i=0; i < numBlocks; i++)
44. {
45. arm_fir_f32(&S, inputF32 + (i * blockSize), outputF32 + (i * blockSize), blockSize);
46. }
48. /* 打印滤波后结果 */
49. for(i=0; i<TEST_LENGTH_SAMPLES; i++)
50. {
51. printf("%frn", testOutput[i]);
52. }
53. }

复制代码

运行如上函数可以通过串口打印出函数arm_fir_f32滤波后的波形数据，下面通过Matlab绘制波形来对比Matlab计算的结果和ARM官方库计算的结果。
    对比前需要先将串口打印出的一组数据加载到Matlab中， arm_fir_f32的计算结果起名sampledata，加载方法在前面的教程中已经讲解过，这里不做赘述了。Matlab中运行的代码如下：

1. %****************************************************************************************
2. %                             FIR低通滤波器设计
3. %***************************************************************************************
4. fs=1000;               %设置采样频率 1K
5. N=320;                %采样点数      
6. n=0:N-1;
7. t=n/fs;                 %时间序列
8. f=n*fs/N;               %频率序列
10. x=sin(2*pi*50*t)+sin(2*pi*200*t);     %50Hz和200Hz正弦波混合
11. b=fir1(28, 0.25);
12. y=filter(b, 1, x);
13. subplot(211);
14. plot(t, y);
15. title('Matlab FIR滤波后的波形');
16. grid on;
18. subplot(212);
19. plot(t, sampledata);
20. title('ARM官方库滤波后的波形');
21. grid on;

复制代码

Matlab运行结果如下：

从上面的波形对比来看，matlab和函数arm_fir_f32计算的结果基本是一致的。为了更好的说明滤波效果，下面从频域的角度来说明这个问题，Matlab上面运行如下代码：

1. %****************************************************************************************
2. %                             FIR低通滤波器设计
3. %***************************************************************************************
4. fs=1000;                  %设置采样频率 1K
5. N=320;                    %采样点数      
6. n=0:N-1;
7. t=n/fs;                      %时间序列
8. f=n*fs/N;                    %频率序列
10. x = sin(2*pi*50*t) + sin(2*pi*200*t);      %50Hz和200Hz正弦波合成
12. subplot(221);
13. plot(t, x);   %绘制信号Mix_Signal的波形                                                
14. xlabel('时间');
15. ylabel('幅值');
16. title('原始信号');
17. grid on;
18.   
19. subplot(222);
20. y=fft(x, N);     %对信号 Mix_Signal做FFT   
21. plot(f,abs(y));
22. xlabel('频率/Hz');
23. ylabel('振幅');
24. title('原始信号FFT');
25. grid on;
27. y3=fft(sampledata, N);            %经过FIR滤波器后得到的信号做FFT
28. subplot(223);                              
29. plot(f,abs(y3));
30. xlabel('频率/Hz');
31. ylabel('振幅');
32. title('滤波后信号FFT');
33. grid on;
35. b=fir1(28, 0.25);          %28阶FIR低通滤波器，截止频率125Hz
36. [H,F]=freqz(b,1,512);      %通过fir1设计的FIR系统的频率响应
37. subplot(224);
38. plot(F/pi,abs(H));          %绘制幅频响应
39. xlabel('归一化频率');        
40. title(['Order=',int2str(30)]);
41. grid on;

复制代码

Matlab显示效果如下:

上面波形变换前的FFT和变换后FFT可以看出，200Hz的正弦波基本被滤除。

席萌0209

37.4 FIR高通滤波器设计



37.4.1 fdatool获取高通滤波器系数


    设计一个如下的例子：
    信号由50Hz正弦波和200Hz正弦波组成，采样率1Kbps，现设计一个高通滤波器，截止频率125Hz，采样320个数据，采用函数fir1进行设计（注意这个函数是基于窗口的方法设计FIR滤波，默认是hamming窗），滤波器阶数设置为28。fadtool的配置如下：

配置好高通滤波器后，具体滤波器系数的生成大家参考本章第二小节的方法即可。

37.4.2 高通滤波器实现


    通过工具箱fdatool获得高通滤波器系数后在开发板上运行函数arm_fir_f32 来测试高通滤波器的效果。

1. #define TEST_LENGTH_SAMPLES  320    /* 采样点数 */
2. #define BLOCK_SIZE           32     /* 调用一次arm_fir_f32处理的采样点个数 */
3. #define NUM_TAPS             29     /* 滤波器系数个数 */
5. uint32_t blockSize = BLOCK_SIZE;
6. uint32_t numBlocks = TEST_LENGTH_SAMPLES/BLOCK_SIZE;            /* 需要调用arm_fir_f32的次数 */
8. static float32_t testInput_f32_50Hz_200Hz[TEST_LENGTH_SAMPLES]; /* 采样点 */
9. static float32_t testOutput[TEST_LENGTH_SAMPLES];               /* 滤波后的输出 */
10. static float32_t firStateF32[BLOCK_SIZE + NUM_TAPS - 1];        /* 状态缓存，大小numTaps + blockSize - 1*/
12. /* 高通滤波器其系数 通过fadtool获取*/
13. const float32_t firCoeffs32HP[NUM_TAPS] = {
14. 0.0018157335f,     0.001582013792f,    -6.107207639e-18f,  -0.003683975432f,   -0.008045346476f,
15. -0.008498443291f,  -1.277260999e-17f,  0.01733288541f,     0.03401865438f,     0.0332348831f,
16. -4.021742543e-17f, -0.06737889349f,    -0.1516391635f,     -0.2220942229f,     0.7486887574f,
17. -0.2220942229f,    -0.1516391635f,     -0.06737889349f,    -4.021742543e-17f,  0.0332348831f,
18. 0.03401865438f,    0.01733288541f,     -1.277260999e-17f,  -0.008498443291f,   -0.008045346476f,
19. -0.003683975432f,  -6.107207639e-18f,  0.001582013792f,    0.0018157335f
20. };
22. /*
23. *********************************************************************************************************
24. *    函 数 名: arm_fir_f32_hp
25. *    功能说明: 调用函数arm_fir_f32_hp实现高通滤波器
26. *    形    参：无
27. *    返 回 值: 无
28. *********************************************************************************************************
29. */
30. static void arm_fir_f32_hp(void)
31. {
32. uint32_t i;
33. arm_fir_instance_f32 S;
34. float32_t  *inputF32, *outputF32;
35.   
36. /* 初始化输入输出缓存指针 */
37. inputF32 = &testInput_f32_50Hz_200Hz[0];
38. outputF32 = &testOutput[0];
40. /* 初始化结构体S */
41. arm_fir_init_f32(&S, NUM_TAPS, (float32_t *)&firCoeffs32HP[0], &firStateF32[0], blockSize);
43. /* 实现FIR滤波 */
44. for(i=0; i < numBlocks; i++)
45. {
46. arm_fir_f32(&S, inputF32 + (i * blockSize), outputF32 + (i * blockSize), blockSize);
47. }
49. /* 打印滤波后结果 */
50. for(i=0; i<TEST_LENGTH_SAMPLES; i++)
51. {
52. printf("%frn", testOutput[i]);
53. }
55. }

复制代码

运行如上函数可以通过串口打印出函数arm_fir_f32滤波后的波形数据，下面通过Matlab绘制波形来对比Matlab计算的结果和ARM官方库计算的结果。
    对比前需要先将串口打印出的一组数据加载到Matlab中， arm_fir_f32的计算结果起名sampledata，加载方法在前面的教程中已经讲解过，这里不做赘述了。Matlab中运行的代码如下：

1. %****************************************************************************************
2. %                             FIR高通滤波器设计
3. %***************************************************************************************
4. fs=1000;                  %设置采样频率 1K
5. N=320;                   %采样点数      
6. n=0:N-1;
7. t=n/fs;                    %时间序列
8. f=n*fs/N;                  %频率序列
10. x=sin(2*pi*50*t)+sin(2*pi*200*t); %50Hz和200Hz正弦波混合           
11. b=fir1(28, 125/500, 'high');     %获得滤波器系数，截止频率125Hz，高通滤波。
12. y=filter(b, 1, x);                %获得滤波后的波形
13. subplot(211);
14. plot(t, y);
15. title('Matlab FIR滤波后的实际波形');
16. grid on;
18. subplot(212);
19. plot(t, sampledata);        %绘制ARM官方库滤波后的波形。
20. title('ARM官方库滤波后的实际波形');
21. grid on;

复制代码

Matlab显示效果如下：

从上面的波形对比来看，matlab和函数arm_fir_f32计算的结果基本是一致的。为了更好的说明滤波效果，下面从频域的角度来说明这个问题，Matlab上面运行如下代码：

1. %****************************************************************************************
2. %                             FIR高通滤波器设计
3. %***************************************************************************************
4. fs=1000;                  %设置采样频率 1K
5. N=320;                   %采样点数      
6. n=0:N-1;
7. t=n/fs;                   %时间序列
8. f=n*fs/N;                 %频率序列
10. x=sin(2*pi*50*t)+sin(2*pi*200*t); %50Hz和200Hz正弦波混合           
11. subplot(221);
12. plot(t, x);                 %绘制信号x的波形                                                
13. xlabel('时间');
14. ylabel('幅值');
15. title('原始信号');
16. grid on;
17.   
18. subplot(222);
19. y=fft(x, N);                %对信号x做FFT   
20. plot(f,abs(y));
21. xlabel('频率/Hz');
22. ylabel('振幅');
23. title('原始信号FFT');
24. grid on;
26. y3=fft(sampledata, N);     %经过FIR滤波器后得到的信号做FFT
27. subplot(223);                              
28. plot(f,abs(y3));
29. xlabel('频率/Hz');
30. ylabel('振幅');
31. title('滤波后信号FFT');
32. grid on;
34. b=fir1(28, 125/500, 'high');   %获得滤波器系数，截止频率125Hz，高通滤波。
35. [H,F]=freqz(b,1,512);         %通过fir1设计的FIR系统的频率响应
36. subplot(224);
37. plot(F/pi,abs(H));            %绘制幅频响应
38. xlabel('归一化频率');        
39. title(['Order=',int2str(30)]);
40. grid on;

复制代码

Matlab显示效果如下：

上面波形变换前的FFT和变换后FFT可以看出，50Hz的正弦波基本被滤除。

席萌0209

37.5 FIR带通滤波器设计



37.5.1 fdatool获取带通滤波器系数


    设计一个如下的例子：
     信号由50Hz正弦波和200Hz正弦波组成，采样率1Kbps，现设计一个带通滤波器，截止频率125Hz和300Hz，采样320个数据，采用函数fir1进行设计（注意这个函数是基于窗口的方法设计FIR滤波，默认是hamming窗），滤波器阶数设置为28。fadtool的配置如下：

配置好带通滤波器后，具体滤波器系数的生成大家参考本章第二小节的方法即可。

37.5.2 带通滤波器实现


    通过工具箱fdatool获得带通滤波器系数后在开发板上运行函数arm_fir_f32 来测试带通滤波器的效果。

1. #define TEST_LENGTH_SAMPLES  320    /* 采样点数 */
2. #define BLOCK_SIZE           32     /* 调用一次arm_fir_f32处理的采样点个数 */
3. #define NUM_TAPS             29     /* 滤波器系数个数 */
5. uint32_t blockSize = BLOCK_SIZE;
6. uint32_t numBlocks = TEST_LENGTH_SAMPLES/BLOCK_SIZE;            /* 需要调用arm_fir_f32的次数 */
8. static float32_t testInput_f32_50Hz_200Hz[TEST_LENGTH_SAMPLES]; /* 采样点 */
9. static float32_t testOutput[TEST_LENGTH_SAMPLES];               /* 滤波后的输出 */
10. static float32_t firStateF32[BLOCK_SIZE + NUM_TAPS - 1];        /* 状态缓存，大小numTaps + blockSize - 1*/
12. /* 带通滤波器系数 通过fadtool获取*/
13. const float32_t firCoeffs32BP[NUM_TAPS] = {
14. 0.003531039227f,    0.0002660876198f,   -0.001947779674f,  0.001266813371f,  -0.008019094355f,
15. -0.01986379735f,    0.01018396299f,     0.03163734451f,    0.00165955862f,   0.03312643617f,
16. 0.0622616075f,      -0.1229852438f,     -0.2399847955f,    0.07637182623f,   0.3482480049f,
17. 0.07637182623f,     -0.2399847955f,     -0.1229852438f,    0.0622616075f,    0.03312643617f,
18. 0.00165955862f,     0.03163734451f,     0.01018396299f,    -0.01986379735f,  -0.008019094355f,
19. 0.001266813371f,   -0.001947779674f,    0.0002660876198f,  0.003531039227f
20. };
22. /*
23. *********************************************************************************************************
24. *    函 数 名: arm_fir_f32_bp
25. *    功能说明: 调用函数arm_fir_f32_bp实现带通滤波器
26. *    形    参：无
27. *    返 回 值: 无
28. *********************************************************************************************************
29. */
30. static void arm_fir_f32_bp(void)
31. {
32. uint32_t i;
33. arm_fir_instance_f32 S;
34. float32_t  *inputF32, *outputF32;
36. /* 初始化输入输出缓存指针 */
37. inputF32 = &testInput_f32_50Hz_200Hz[0];
38. outputF32 = &testOutput[0];
40. /* 初始化结构体S */
41. arm_fir_init_f32(&S, NUM_TAPS, (float32_t *)&firCoeffs32BP[0], &firStateF32[0], blockSize);
43. /* 实现FIR滤波 */
44. for(i=0; i < numBlocks; i++)
45. {
46. arm_fir_f32(&S, inputF32 + (i * blockSize), outputF32 + (i * blockSize), blockSize);
47. }
49. /* 打印滤波后结果 */
50. for(i=0; i<TEST_LENGTH_SAMPLES; i++)
51. {
52. printf("%frn", testOutput[i]);
53. }
54. }

复制代码

运行如上函数可以通过串口打印出函数arm_fir_f32滤波后的波形数据，下面通过Matlab绘制波形来对比Matlab计算的结果和ARM官方库计算的结果。
    对比前需要先将串口打印出的一组数据加载到Matlab中， arm_fir_f32的计算结果起名sampledata，加载方法在前面的教程中已经讲解过，这里不做赘述了。Matlab中运行的代码如下：

1. %****************************************************************************************
2. %                             FIR带通滤波器设计
3. %***************************************************************************************
4. fs=1000;                   %设置采样频率 1K
5. N=320;                    %采样点数      
6. n=0:N-1;
7. t=n/fs;                     %时间序列
8. f=n*fs/N;                  %频率序列
10. x=sin(2*pi*50*t)+sin(2*pi*200*t);  %50Hz和200Hz正弦波混合           
11. b=fir1(28, [125/500 300/500]);     %获得滤波器系数，截止频率125Hz和300Hz，带通滤波。
12. y=filter(b, 1, x);                   %获得滤波后的波形
13. subplot(211);
14. plot(t, y);
15. title('Matlab FIR滤波后的实际波形');
16. grid on;
18. subplot(212);
19. plot(t, sampledata);        %绘制ARM官方库滤波后的波形。
20. title('ARM官方库滤波后的实际波形');
21. grid on;

复制代码

Matlab显示效果如下：

从上面的波形对比来看，matlab和函数arm_fir_f32计算的结果基本是一致的。为了更好的说明滤波效果，下面从频域的角度来说明这个问题，Matlab上面运行如下代码：

1. %****************************************************************************************
2. %                             FIR带通滤波器设计
3. %***************************************************************************************
4. fs=1000;                   %设置采样频率 1K
5. N=320;                    %采样点数      
6. n=0:N-1;
7. t=n/fs;                    %时间序列
8. f=n*fs/N;                  %频率序列
10. x=sin(2*pi*50*t)+sin(2*pi*200*t);  %50Hz和200Hz正弦波混合               
11. subplot(221);
12. plot(t, x);       %绘制信号x的波形                                                
13. xlabel('时间');
14. ylabel('幅值');
15. title('原始信号');
16. grid on;
17.   
18. subplot(222);
19. y=fft(x, N);     %对信号x做FFT   
20. plot(f,abs(y));
21. xlabel('频率/Hz');
22. ylabel('振幅');
23. title('原始信号FFT');
24. grid on;
26. y3=fft(sampledata, N);       %经过FIR滤波器后得到的信号做FFT
27. subplot(223);                              
28. plot(f,abs(y3));
29. xlabel('频率/Hz');
30. ylabel('振幅');
31. title('滤波后信号FFT');
32. grid on;
34. b=fir1(28, [125/500 300/500]);   %获得滤波器系数，截止频率125Hz，高通滤波。
35. [H,F]=freqz(b,1,160);            %通过fir1设计的FIR系统的频率响应
36. subplot(224);
37. plot(F/pi,abs(H));                %绘制幅频响应
38. xlabel('归一化频率');        
39. title(['Order=',int2str(28)]);
40. grid on;

复制代码

上面波形变换前的FFT和变换后FFT可以看出，50Hz的正弦波基本被滤除。

席萌0209

37.6 FIR带阻滤波器设计



37.6.1 fdatool获取带阻滤波器系数


    设计一个如下的例子：
    信号由50Hz正弦波和200Hz正弦波组成，采样率1Kbps，现设计一个带阻滤波器，截止频率125Hz和300Hz，采样320个数据，采用函数fir1进行设计（注意这个函数是基于窗口的方法设计FIR滤波，默认是hamming窗），滤波器阶数设置为28。fadtool的配置如下：

配置好带阻滤波器后，具体滤波器系数的生成大家参考本章第二小节的方法即可。

37.6.2 带阻滤波器实现


通过工具箱fdatool获得带阻滤波器系数后在开发板上运行函数arm_fir_f32 来测试带通滤波器的效果。

1. #define TEST_LENGTH_SAMPLES  320    /* 采样点数 */
2. #define BLOCK_SIZE           32     /* 调用一次arm_fir_f32处理的采样点个数 */
3. #define NUM_TAPS             29     /* 滤波器系数个数 */
5. uint32_t blockSize = BLOCK_SIZE;
6. uint32_t numBlocks = TEST_LENGTH_SAMPLES/BLOCK_SIZE;            /* 需要调用arm_fir_f32的次数 */
8. static float32_t testInput_f32_50Hz_200Hz[TEST_LENGTH_SAMPLES]; /* 采样点 */
9. static float32_t testOutput[TEST_LENGTH_SAMPLES];               /* 滤波后的输出 */
10. static float32_t firStateF32[BLOCK_SIZE + NUM_TAPS - 1];        /* 状态缓存，大小numTaps + blockSize - 1*/
12. /* 带阻滤波器系数 通过fadtool获取*/
13. const float32_t firCoeffs32BPCheb[NUM_TAPS] = {
14. 0.01801843569f,    0.0007182828849f,  -0.004868913442f,  0.002710500965f,  -0.01462193858f,
15. -0.03147283196f,   0.01435638033f,    0.04055848345f,    0.00197162549f,   0.03706155345f,
16. 0.06650412083f,    -0.1269270927f,    -0.2418768406f,    0.07591249049f,   0.3445736468f,
17. 0.07591249049f,    -0.2418768406f,    -0.1269270927f,    0.06650412083f,   0.03706155345f,
18. 0.00197162549f,    0.04055848345f,    0.01435638033f,    -0.03147283196f,  -0.01462193858f,
19. 0.002710500965f,  -0.004868913442f,   0.0007182828849f,  0.01801843569f
20. };*/
21. const float32_t firCoeffs32BS[NUM_TAPS] = {
22. -0.003560454352f,  -0.0002683042258f,  0.001964005642f,   -0.001277366537f,   0.008085897192f,
23. 0.02002927102f,    -0.01026879996f,    -0.03190089762f,   -0.001673383522f,   -0.0334023945f,
24. -0.06278027594f,   0.1240097657f,      0.2419839799f,     -0.07700803876f,    0.6521340013f,
25. -0.07700803876f,   0.2419839799f,      0.1240097657f,     -0.06278027594f,    -0.0334023945f,
26. -0.001673383522f,  -0.03190089762f,    -0.01026879996f,   0.02002927102f,     0.008085897192f,
27. -0.001277366537f,  0.001964005642f,    -0.0002683042258f, -0.003560454352f
28. };
30. /*
31. *********************************************************************************************************
32. *    函 数 名: arm_fir_f32_bs
33. *    功能说明: 调用函数arm_fir_f32_bs实现带阻滤波器
34. *    形    参：无
35. *    返 回 值: 无
36. *********************************************************************************************************
37. */
38. static void arm_fir_f32_bs(void)
39. {
40. uint32_t i;
41. arm_fir_instance_f32 S;
42. float32_t  *inputF32, *outputF32;
44. /* 初始化输入输出缓存指针 */
45. inputF32 = &testInput_f32_50Hz_200Hz[0];
46. outputF32 = &testOutput[0];
48. /* 初始化结构体S */
49. arm_fir_init_f32(&S, NUM_TAPS, (float32_t *)&firCoeffs32BS[0], &firStateF32[0], blockSize);
51. /* 实现FIR滤波 */
52. for(i=0; i < numBlocks; i++)
53. {
54. arm_fir_f32(&S, inputF32 + (i * blockSize), outputF32 + (i * blockSize), blockSize);
55. }
57. /* 打印滤波后结果 */
58. for(i=0; i<TEST_LENGTH_SAMPLES; i++)
59. {
60. printf("%frn", testOutput[i]);
61. }
62. }

复制代码

运行如上函数可以通过串口打印出函数arm_fir_f32滤波后的波形数据，下面通过Matlab绘制波形来对比Matlab计算的结果和ARM官方库计算的结果。
对比前需要先将串口打印出的一组数据加载到Matlab中， arm_fir_f32的计算结果起名sampledata，加载方法在前面的教程中已经讲解过，这里不做赘述了。Matlab中运行的代码如下：

1. %****************************************************************************************
2. %                             FIR带阻滤波器设计
3. %***************************************************************************************
4. fs=1000;                  %设置采样频率 1K
5. N=320;                   %采样点数      
6. n=0:N-1;
7. t=n/fs;                    %时间序列
8. f=n*fs/N;                  %频率序列
10. x=sin(2*pi*50*t)+sin(2*pi*200*t);       %50Hz和200Hz正弦波混合           
11. b=fir1(28, [125/500 300/500], 'stop');   %获得滤波器系数，截止频率125Hz和300，带阻滤波。
12. y=filter(b, 1, x);                        %获得滤波后的波形
13. subplot(211);
14. plot(t, y);
15. title('Matlab FIR滤波后的实际波形');
16. grid on;
18. subplot(212);
19. plot(t, sampledata);        %绘制ARM官方库滤波后的波形。
20. title('ARM官方库滤波后的实际波形');
21. grid on;

复制代码

Matlab运行结果如下：

从上面的波形对比来看，matlab和函数arm_fir_f32计算的结果基本是一致的。为了更好的说明滤波效果，下面从频域的角度来说明这个问题，Matlab上面运行如下代码：

1. %****************************************************************************************
2. %                             FIR带阻滤波器设计
3. %***************************************************************************************
4. fs=1000;                   %设置采样频率 1K
5. N=320;                    %采样点数      
6. n=0:N-1;
7. t=n/fs;                    %时间序列
8. f=n*fs/N;                  %频率序列
10. x=sin(2*pi*50*t)+sin(2*pi*200*t);  %50Hz和200Hz正弦波混合           
11. subplot(221);
12. plot(t, x);   %绘制信号x的波形                                                
13. xlabel('时间');
14. ylabel('幅值');
15. title('原始信号');
16. grid on;
17.   
18. subplot(222);
19. y=fft(x, N);     %对信号x做FFT   
20. plot(f,abs(y));
21. xlabel('频率/Hz');
22. ylabel('振幅');
23. title('原始信号FFT');
24. grid on;
26. y3=fft(sampledata, N);       %经过FIR滤波器后得到的信号做FFT
27. subplot(223);                              
28. plot(f,abs(y3));
29. xlabel('频率/Hz');
30. ylabel('振幅');
31. title('滤波后信号FFT');
32. grid on;
34. b=fir1(28, [125/500 300/500], 'stop');  %获得滤波器系数，截止频率125Hz和300Hz，带阻滤波。     
35. [H,F]=freqz(b,1,160);                  %通过fir1设计的FIR系统的频率响应
36. subplot(224);
37. plot(F/pi,abs(H));             %绘制幅频响应
38. xlabel('归一化频率');        
39. title(['Order=',int2str(28)]);
40. grid on;

复制代码

Matlab运行效果如下：

上面波形变换前的FFT和变换后FFT可以看出，200Hz的正弦波基本被滤除。

席萌0209

37.7 切比雪夫窗口设计带通滤波器



37.7.1 fdatool获取滤波器系数


    设计一个如下的例子：
    信号由50Hz正弦波和200Hz正弦波组成，采样率1Kbps，现设计一个使用切比雪夫窗口的带通滤波器，截止频率125Hz和300Hz，切比雪夫波纹设置为30db，采样320个数据，采用函数fir1进行设计（注意这个函数是基于窗口的方法设计FIR滤波，默认是hamming窗），滤波器阶数设置为28。fadtool的配置如下：

配置好带通滤波器后，具体滤波器系数的生成大家参考本章第二小节的方法即可。

37.7.2 带通滤波器实现


    通过工具箱fdatool获得带通滤波器系数后在开发板上运行函数arm_fir_f32 来测试带通滤波器的效果。

1. #define TEST_LENGTH_SAMPLES  320    /* 采样点数 */
2. #define BLOCK_SIZE           32     /* 调用一次arm_fir_f32处理的采样点个数 */
3. #define NUM_TAPS             29     /* 滤波器系数个数 */
5. uint32_t blockSize = BLOCK_SIZE;
6. uint32_t numBlocks = TEST_LENGTH_SAMPLES/BLOCK_SIZE;            /* 需要调用arm_fir_f32的次数 */
8. static float32_t testInput_f32_50Hz_200Hz[TEST_LENGTH_SAMPLES]; /* 采样点 */
9. static float32_t testOutput[TEST_LENGTH_SAMPLES];               /* 滤波后的输出 */
10. static float32_t firStateF32[BLOCK_SIZE + NUM_TAPS - 1];        /* 状态缓存，大小numTaps + blockSize - 1*/
11. /* 带通滤波器系数 切比雪夫窗口 通过fadtool获取*/
12. const float32_t firCoeffs32BPCheb[NUM_TAPS] = {
13. 0.01801843569f,    0.0007182828849f,  -0.004868913442f,  0.002710500965f,  -0.01462193858f,
14. -0.03147283196f,   0.01435638033f,    0.04055848345f,    0.00197162549f,   0.03706155345f,
15. 0.06650412083f,    -0.1269270927f,    -0.2418768406f,    0.07591249049f,   0.3445736468f,
16. 0.07591249049f,    -0.2418768406f,    -0.1269270927f,    0.06650412083f,   0.03706155345f,
17. 0.00197162549f,    0.04055848345f,    0.01435638033f,    -0.03147283196f,  -0.01462193858f,
18. 0.002710500965f,  -0.004868913442f,   0.0007182828849f,  0.01801843569f
19. };
21. /*
22. *********************************************************************************************************
23. *    函 数 名: arm_fir_f32_bp
24. *    功能说明: 调用函数arm_fir_f32_bp实现带通滤波器
25. *    形    参：无
26. *    返 回 值: 无
27. *********************************************************************************************************
28. */
29. static void arm_fir_f32_bp(void)
30. {
31. uint32_t i;
32. arm_fir_instance_f32 S;
33. float32_t  *inputF32, *outputF32;
35. /* 初始化输入输出缓存指针 */
36. inputF32 = &testInput_f32_50Hz_200Hz[0];
37. outputF32 = &testOutput[0];
39. /* 初始化结构体S */
40. arm_fir_init_f32(&S, NUM_TAPS, (float32_t *)&firCoeffs32BP[0], &firStateF32[0], blockSize);
42. /* 实现FIR滤波 */
43. for(i=0; i < numBlocks; i++)
44. {
45. arm_fir_f32(&S, inputF32 + (i * blockSize), outputF32 + (i * blockSize), blockSize);
46. }
48. /* 打印滤波后结果 */
49. for(i=0; i<TEST_LENGTH_SAMPLES; i++)
50. {
51. printf("%frn", testOutput[i]);
52. }
53. }

复制代码

运行如上函数可以通过串口打印出函数arm_fir_f32滤波后的波形数据，下面通过Matlab绘制波形来对比Matlab计算的结果和ARM官方库计算的结果。
    对比前需要先将串口打印出的一组数据加载到Matlab中， arm_fir_f32的计算结果起名sampledata，加载方法在前面的教程中已经讲解过，这里不做赘述了。Matlab中运行的代码如下：

1. %****************************************************************************************
2. %                             FIR带通滤波器设计，切比雪夫窗口
3. %***************************************************************************************
4. fs=1000;                  %设置采样频率 1K
5. N=320;                   %采样点数      
6. n=0:N-1;
7. t=n/fs;                    %时间序列
8. f=n*fs/N;                 %频率序列
10. x=sin(2*pi*50*t)+sin(2*pi*200*t);          %50Hz和200Hz正弦波混合      
11. Window = chebwin(29, 30);               %30db的切比雪夫窗
12. b=fir1(28, [125/500 300/500], Window);   %获得滤波器系数，截止频率125Hz和300Hz，带通滤波。
13. y=filter(b, 1, x);                           %获得滤波后的波形
14. subplot(211);
15. plot(t, y);
16. title('Matlab FIR滤波后的实际波形');
17. grid on;
19. subplot(212);
20. plot(t, sampledata);                 %绘制ARM官方库滤波后的波形。
21. title('ARM官方库滤波后的实际波形');
22. grid on;

复制代码

Matlab的运行效果如下：

从上面的波形对比来看，matlab和函数arm_fir_f32计算的结果基本是一致的。为了更好的说明滤波效果，下面从频域的角度来说明这个问题，Matlab上面运行如下代码：

1. %****************************************************************************************
2. %                             FIR带通滤波器设计，切比雪夫窗口
3. %***************************************************************************************
4. fs=1000;                   %设置采样频率 1K
5. N=320;                    %采样点数      
6. n=0:N-1;
7. t=n/fs;                     %时间序列
8. f=n*fs/N;                  %频率序列
11. x=sin(2*pi*50*t)+sin(2*pi*200*t);  %50Hz和200Hz正弦波混合               
12. subplot(221);
13. plot(t, x);   %绘制信号x的波形                                                
14. xlabel('时间');
15. ylabel('幅值');
16. title('原始信号');
17. grid on;
18.   
19. subplot(222);
20. y=fft(x, N);     %对信号x做FFT   
21. plot(f,abs(y));
22. xlabel('频率/Hz');
23. ylabel('振幅');
24. title('原始信号FFT');
25. grid on;
27. y3=fft(sampledata, N);       %经过FIR滤波器后得到的信号做FFT
28. subplot(223);                              
29. plot(f,abs(y3));
30. xlabel('频率/Hz');
31. ylabel('振幅');
32. title('滤波后信号FFT');
33. grid on;
35. Window = chebwin(29, 30);       %30db的切比雪夫窗
36. b=fir1(28, [125/500 300/500], Window);  %获得滤波器系数，截止频率125Hz和300Hz，带通滤波。
37. [H,F]=freqz(b,1,160);                    %通过fir1设计的FIR系统的频率响应
38. subplot(224);
39. plot(F/pi,abs(H));                       %绘制幅频响应
40. xlabel('归一化频率');        
41. title(['Order=',int2str(28)]);
42. grid on;

复制代码

Matlab运行结果如下：

上面波形变换前的FFT和变换后FFT可以看出，50Hz的正弦波基本被滤除，在归一化频率中我们可以看到一定的波纹

席萌0209

37.8 FIR滤波后的群延迟


    波形经过FIR滤波器后，输出的波形会有一定的延迟。对于线性相位的FIR，这个群延迟就是一个常数。但是实际应用中这个群延迟是多少呢，关于群延迟的数值，fdatool工具箱会根据用户的配置计算好。
    比如咱们前面设计的28阶FIR高通，低通，带通和带阻滤波器的群延迟就是14，反应在实际的采样值上就是滤波后输出数据的第15个才是实际滤波后的波形数据起始点。
    下面是看群延迟采样点的位置：

细心的读者可能发现全面做低通，高通，带通和带阻滤波后，输出的波形前面几个点感觉有问题，其实就是群延迟造成的。
    为了更好的说明这个问题，下面再使用Matlab举一个低通和一个高通滤波的例子：信号由50Hz正弦波和200Hz正弦波组成，采样率1Kbps，截止频率125Hz，采样320个数据，采用函数fir1进行设计，滤波器阶数设置为28。下面是低通滤波器的Matlab代码，将原始信号从第一个点开始显示，而滤波后的信号从群延迟后的第15个点开始显示：

1. fs=1000;                  %设置采样频率 1K
2. N=320;                    %采样点数      
3. n=0:N-1;
4. t=n/fs;                     %时间序列
5. f=n*fs/N;                  %频率序列
7. x1=sin(2*pi*50*t);
8. x2=sin(2*pi*200*t);
9. x=sin(2*pi*50*t)+sin(2*pi*200*t);  %50Hz和200Hz正弦波混合               
11. plot(n, x1, 'b');   %绘制信号x的波形  
12. xlabel('时间');
13. ylabel('幅值');
14. title('原始信号和滤波后信号');
15. hold on;
16.   
17. b=fir1(28, 125/500);     %获得滤波器系数，截止频率125Hz.
18. y=filter(b, 1, x);
19. plot(n(1:305), y(15:319), 'r');
20. legend('原始信号','滤波后信号');
21. grid on;

复制代码

Matlab的运行结果如下：

可以看出，显示波形基本重合，这个说明14个采样点的群延迟是正确的。下面同样使用上面的那个例子实现一个高通滤波器，截止频率是125Hz，阶数同样设置为28，将原始信号从第一个点开始显示，而滤波后的信号从群延迟后的第15个点开始显示，Matlab运行代码如下：

1. fs=1000;                  %设置采样频率 1K
2. N=320;                   %采样点数      
3. n=0:N-1;
4. t=n/fs;                    %时间序列
5. f=n*fs/N;                  %频率序列
7. x1=sin(2*pi*50*t);
8. x2=sin(2*pi*200*t);
9. x=sin(2*pi*50*t)+sin(2*pi*200*t);  %50Hz和200Hz正弦波混合               
11. plot(n, x2, 'b');   %绘制信号x的波形  
12. xlabel('时间');
13. ylabel('幅值');
14. title('原始信号和滤波后信号');
15. hold on;
16.   
17. b=fir1(28, 125/500, 'high');     %获得滤波器系数，截止频率125Hz.
18. y=filter(b, 1, x);
19. plot(n(1:305), y(15:319), 'r');
20. legend('原始信号','滤波后信号');
21. grid on;

复制代码

Matlab运行结果如下：

可以看出，显示波形基本重合，这个说明14个采样点的群延迟也是是正确的。大家在使用FIR滤波器的时候一定要注意这个问题。

席萌0209

37.9 总结


    本章节主要讲解了FIR滤波器的低通，高通，带通和带阻滤波器的实现，同时一定要注意线性相位FIR滤波器的群延迟问题。

xgqxgq

您好，我想问下，能把所有的数据格式换成u16的类型吗，浮点数32位我之前有弄过，够没问题，可是改了U16的就实现不了滤波了

eric2013

xgqxgq:您好，我想问下，能把所有的数据格式换成u16的类型吗，浮点数32位我之前有弄过，够没问题，可是改了U16的就实现不了滤波了
 (2016-07-08 13:44) 

有q15的api函数，研究下。

2602082487

你好，我想用f4把10K和200Hz的混叠中滤除200Hz的信号，能做到吗？

eric2013

2602082487:你好，我想用f4把10K和200Hz的混叠中滤除200Hz的信号，能做到吗？ (2016-08-31 15:08) 

完全没问题。做一个高通滤波器就可以了。

2602082487

eric2013:完全没问题。做一个高通滤波器就可以了。 (2016-08-31 15:12) 

嗯，实时滤波f4的计算能力够哈，谢谢

单枪舞九州

matlab是个好帮手！