你的浏览器版本过低,可能导致网站不能正常访问!
为了你能正常使用网站功能,请使用这些浏览器。

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

[复制链接]
baiyongbin2009 发布时间:2015-4-24 11:03
特别说明:完整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]
结构图如下:
37.1.png
这种网络结构就是在35.2.1小节所讲的直接型结构。

1 收藏 2 评论11 发布时间:2015-4-24 11:03

举报

11个回答
baiyongbin2009 回答时间:2015-4-24 11:07:33
37.2 Matlab工具箱fdatool生成C头文件
    下面我们讲解下如何通过fdatool工具生成C头文件,也就是生成滤波器系数。首先在matlab的命
窗口输入fadtool就能打开这个工具箱:
37.2.png
fadtool界面打开效果如下:
37.3.png
FIR滤波器的低通,高通,带通,带阻滤波的设置会在下面一一讲解,这里说一下设置后相应参数后如何生成滤波器系数。参数设置好以后点击如下按钮:
37.4.png
点击Design Filter按钮以后就生成了所需的滤波器系数,生成滤波器系数以后点击fadtool界面上的菜单Targets->Generate C header ,打开后显示如下界面:
37.5.png
然后点击Generate,生成如下界面:
37.6.png
再点击保存,并打开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. */

  9. /*
  10. * Discrete-Time FIR Filter (real)
  11. * -------------------------------
  12. * Filter Structure  : Direct-Form FIR
  13. * Filter Length     : 51
  14. * Stable            : Yes
  15. * Linear Phase      : Yes (Type 1)
  16. */

  17. /* General type conversion for MATLAB generated C-code  */
  18. #include "tmwtypes.h"
  19. /*
  20. * Expected path to tmwtypes.h
  21. * C:\Program Files\MATLAB\R2012a\extern\include\tmwtypes.h
  22. */
  23. /*
  24. * Warning - Filter coefficients were truncated to fit specified data type.  
  25. *   The resulting response may not match generated theoretical response.
  26. *   Use the Filter Design & Analysis Tool to design accurate
  27. *   single-precision filter coefficients.
  28. */
  29. const int BL = 51;
  30. const real32_T B[51] = {
  31.   -0.0009190982091, -0.00271769613,-0.002486952813, 0.003661438357,   0.0136509249,
  32.     0.01735116541,  0.00766530633,-0.006554719061,-0.007696784101, 0.006105459295,
  33.     0.01387391612,0.0003508617228, -0.01690892503,-0.008905642666,  0.01744112931,
  34.     0.02074504457,  -0.0122964941, -0.03424086422,-0.001034529647,  0.04779030383,
  35.     0.02736303769, -0.05937951803, -0.08230702579,  0.06718690693,   0.3100151718,
  36.      0.4300478697,   0.3100151718,  0.06718690693, -0.08230702579, -0.05937951803,
  37.     0.02736303769,  0.04779030383,-0.001034529647, -0.03424086422,  -0.0122964941,
  38.     0.02074504457,  0.01744112931,-0.008905642666, -0.01690892503,0.0003508617228,
  39.     0.01387391612, 0.006105459295,-0.007696784101,-0.006554719061,  0.00766530633,
  40.     0.01735116541,   0.0136509249, 0.003661438357,-0.002486952813, -0.00271769613,
  41.   -0.0009190982091
  42. };
复制代码
上面数组B[51]中的数据就是滤波器系数。下面小节讲解如何使用fdatool配置FIR低通,高通,带通和带阻滤波。关于fdatool的其它用法,大家可以在matlab命令窗口中输入help fadtool打开帮助文档进行学习。

baiyongbin2009 回答时间:2015-4-24 11:15:26
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.7.png
配置好低通滤波器后,具体滤波器系数的生成大家参考本章第二小节的方法即可。

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     /* 滤波器系数个数 */

  4. uint32_t blockSize = BLOCK_SIZE;
  5. uint32_t numBlocks = TEST_LENGTH_SAMPLES/BLOCK_SIZE;            /* 需要调用arm_fir_f32的次数 */

  6. static float32_t testInput_f32_50Hz_200Hz[TEST_LENGTH_SAMPLES]; /* 采样点 */
  7. static float32_t testOutput[TEST_LENGTH_SAMPLES];               /* 滤波后的输出 */
  8. static float32_t firStateF32[BLOCK_SIZE + NUM_TAPS - 1];        /* 状态缓存,大小numTaps + blockSize - 1*/

  9. /* 低通滤波器系数 通过fadtool获取*/
  10. const float32_t firCoeffs32LP[NUM_TAPS] = {
  11.   -0.001822523074f,  -0.001587929321f,  1.226008847e-18f,  0.003697750857f,  0.008075430058f,
  12.   0.008530221879f,   -4.273456581e-18f, -0.01739769801f,   -0.03414586186f,  -0.03335915506f,
  13.   8.073562366e-18f,  0.06763084233f,    0.1522061825f,     0.2229246944f,    0.2504960895f,
  14.   0.2229246944f,     0.1522061825f,     0.06763084233f,    8.073562366e-18f, -0.03335915506f,
  15.   -0.03414586186f,   -0.01739769801f,   -4.273456581e-18f, 0.008530221879f,  0.008075430058f,
  16.   0.003697750857f,   1.226008847e-18f,  -0.001587929321f,  -0.001822523074f
  17. };
  18. /*
  19. *********************************************************************************************************
  20. *        函 数 名: arm_fir_f32_lp
  21. *        功能说明: 调用函数arm_fir_f32_lp实现低通滤波器
  22. *        形    参:无
  23. *        返 回 值: 无
  24. *********************************************************************************************************
  25. */
  26. static void arm_fir_f32_lp(void)
  27. {
  28. uint32_t i;
  29. arm_fir_instance_f32 S;
  30. float32_t  *inputF32, *outputF32;

  31. /* 初始化输入输出缓存指针 */
  32. inputF32 = &testInput_f32_50Hz_200Hz[0];
  33. outputF32 = &testOutput[0];

  34. /* 初始化结构体S */
  35. arm_fir_init_f32(&S, NUM_TAPS, (float32_t *)&firCoeffs32LP[0], &firStateF32[0], blockSize);

  36. /* 实现FIR滤波 */
  37. for(i=0; i < numBlocks; i++)
  38. {
  39. arm_fir_f32(&S, inputF32 + (i * blockSize), outputF32 + (i * blockSize), blockSize);
  40. }

  41. /* 打印滤波后结果 */
  42. for(i=0; i<TEST_LENGTH_SAMPLES; i++)
  43. {
  44. printf("%f\r\n", testOutput[i]);
  45. }
  46. }
复制代码
运行如上函数可以通过串口打印出函数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;               %频率序列

  9. x=sin(2*pi*50*t)+sin(2*pi*200*t);     %50Hz和200Hz正弦波混合
  10. b=fir1(28, 0.25);
  11. y=filter(b, 1, x);
  12. subplot(211);
  13. plot(t, y);
  14. title('Matlab FIR滤波后的波形');
  15. grid on;

  16. subplot(212);
  17. plot(t, sampledata);
  18. title('ARM官方库滤波后的波形');
  19. grid on;
复制代码
Matlab运行结果如下:
37.8.png
从上面的波形对比来看,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;                    %频率序列

  9. x = sin(2*pi*50*t) + sin(2*pi*200*t);      %50Hz和200Hz正弦波合成

  10. subplot(221);
  11. plot(t, x);   %绘制信号Mix_Signal的波形                                                
  12. xlabel('时间');
  13. ylabel('幅值');
  14. title('原始信号');
  15. grid on;
  16.   
  17. subplot(222);
  18. y=fft(x, N);     %对信号 Mix_Signal做FFT   
  19. plot(f,abs(y));
  20. xlabel('频率/Hz');
  21. ylabel('振幅');
  22. title('原始信号FFT');
  23. grid on;

  24. y3=fft(sampledata, N);            %经过FIR滤波器后得到的信号做FFT
  25. subplot(223);                              
  26. plot(f,abs(y3));
  27. xlabel('频率/Hz');
  28. ylabel('振幅');
  29. title('滤波后信号FFT');
  30. grid on;

  31. b=fir1(28, 0.25);          %28阶FIR低通滤波器,截止频率125Hz
  32. [H,F]=freqz(b,1,512);      %通过fir1设计的FIR系统的频率响应
  33. subplot(224);
  34. plot(F/pi,abs(H));          %绘制幅频响应
  35. xlabel('归一化频率');        
  36. title(['Order=',int2str(30)]);
  37. grid on;
复制代码
Matlab显示效果如下:
37.9.png
上面波形变换前的FFT和变换后FFT可以看出,200Hz的正弦波基本被滤除。

stary666 回答时间:2015-4-24 11:21:35
沙发,支持原创
baiyongbin2009 回答时间:2015-4-24 12:26:49
37.4 FIR高通滤波器设计
37.4.1 fdatool获取高通滤波器系数
    设计一个如下的例子:
    信号由50Hz正弦波和200Hz正弦波组成,采样率1Kbps,现设计一个高通滤波器,截止频率125Hz,采样320个数据,采用函数fir1进行设计(注意这个函数是基于窗口的方法设计FIR滤波,默认是hamming窗),滤波器阶数设置为28。fadtool的配置如下:
37.10.png
配置好高通滤波器后,具体滤波器系数的生成大家参考本章第二小节的方法即可。

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     /* 滤波器系数个数 */

  4. uint32_t blockSize = BLOCK_SIZE;
  5. uint32_t numBlocks = TEST_LENGTH_SAMPLES/BLOCK_SIZE;            /* 需要调用arm_fir_f32的次数 */

  6. static float32_t testInput_f32_50Hz_200Hz[TEST_LENGTH_SAMPLES]; /* 采样点 */
  7. static float32_t testOutput[TEST_LENGTH_SAMPLES];               /* 滤波后的输出 */
  8. static float32_t firStateF32[BLOCK_SIZE + NUM_TAPS - 1];        /* 状态缓存,大小numTaps + blockSize - 1*/

  9. /* 高通滤波器其系数 通过fadtool获取*/
  10. const float32_t firCoeffs32HP[NUM_TAPS] = {
  11. 0.0018157335f,     0.001582013792f,    -6.107207639e-18f,  -0.003683975432f,   -0.008045346476f,
  12. -0.008498443291f,  -1.277260999e-17f,  0.01733288541f,     0.03401865438f,     0.0332348831f,
  13. -4.021742543e-17f, -0.06737889349f,    -0.1516391635f,     -0.2220942229f,     0.7486887574f,
  14. -0.2220942229f,    -0.1516391635f,     -0.06737889349f,    -4.021742543e-17f,  0.0332348831f,
  15. 0.03401865438f,    0.01733288541f,     -1.277260999e-17f,  -0.008498443291f,   -0.008045346476f,
  16. -0.003683975432f,  -6.107207639e-18f,  0.001582013792f,    0.0018157335f
  17. };

  18. /*
  19. *********************************************************************************************************
  20. *        函 数 名: arm_fir_f32_hp
  21. *        功能说明: 调用函数arm_fir_f32_hp实现高通滤波器
  22. *        形    参:无
  23. *        返 回 值: 无
  24. *********************************************************************************************************
  25. */
  26. static void arm_fir_f32_hp(void)
  27. {
  28. uint32_t i;
  29. arm_fir_instance_f32 S;
  30. float32_t  *inputF32, *outputF32;
  31.   
  32. /* 初始化输入输出缓存指针 */
  33. inputF32 = &testInput_f32_50Hz_200Hz[0];
  34. outputF32 = &testOutput[0];

  35. /* 初始化结构体S */
  36. arm_fir_init_f32(&S, NUM_TAPS, (float32_t *)&firCoeffs32HP[0], &firStateF32[0], blockSize);

  37. /* 实现FIR滤波 */
  38. for(i=0; i < numBlocks; i++)
  39. {
  40. arm_fir_f32(&S, inputF32 + (i * blockSize), outputF32 + (i * blockSize), blockSize);
  41. }

  42. /* 打印滤波后结果 */
  43. for(i=0; i<TEST_LENGTH_SAMPLES; i++)
  44. {
  45. printf("%f\r\n", testOutput[i]);
  46. }

  47. }
复制代码
运行如上函数可以通过串口打印出函数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;                  %频率序列

  9. x=sin(2*pi*50*t)+sin(2*pi*200*t); %50Hz和200Hz正弦波混合           
  10. b=fir1(28, 125/500, 'high');     %获得滤波器系数,截止频率125Hz,高通滤波。
  11. y=filter(b, 1, x);                %获得滤波后的波形
  12. subplot(211);
  13. plot(t, y);
  14. title('Matlab FIR滤波后的实际波形');
  15. grid on;

  16. subplot(212);
  17. plot(t, sampledata);        %绘制ARM官方库滤波后的波形。
  18. title('ARM官方库滤波后的实际波形');
  19. grid on;
复制代码
Matlab显示效果如下:
37.11.png
从上面的波形对比来看,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;                 %频率序列

  9. x=sin(2*pi*50*t)+sin(2*pi*200*t); %50Hz和200Hz正弦波混合           
  10. subplot(221);
  11. plot(t, x);                 %绘制信号x的波形                                                
  12. xlabel('时间');
  13. ylabel('幅值');
  14. title('原始信号');
  15. grid on;
  16.   
  17. subplot(222);
  18. y=fft(x, N);                %对信号x做FFT   
  19. plot(f,abs(y));
  20. xlabel('频率/Hz');
  21. ylabel('振幅');
  22. title('原始信号FFT');
  23. grid on;

  24. y3=fft(sampledata, N);     %经过FIR滤波器后得到的信号做FFT
  25. subplot(223);                              
  26. plot(f,abs(y3));
  27. xlabel('频率/Hz');
  28. ylabel('振幅');
  29. title('滤波后信号FFT');
  30. grid on;

  31. b=fir1(28, 125/500, 'high');   %获得滤波器系数,截止频率125Hz,高通滤波。
  32. [H,F]=freqz(b,1,512);         %通过fir1设计的FIR系统的频率响应
  33. subplot(224);
  34. plot(F/pi,abs(H));            %绘制幅频响应
  35. xlabel('归一化频率');        
  36. title(['Order=',int2str(30)]);
  37. grid on;
复制代码
Matlab显示效果如下:
37.12.png
上面波形变换前的FFT和变换后FFT可以看出,50Hz的正弦波基本被滤除。

baiyongbin2009 回答时间:2015-4-24 12:32:00
37.5 FIR带通滤波器设计
37.5.1 fdatool获取带通滤波器系数
    设计一个如下的例子:
    信号由50Hz正弦波和200Hz正弦波组成,采样率1Kbps,现设计一个带通滤波器,截止频率125Hz和300Hz,采样320个数据,采用函数fir1进行设计(注意这个函数是基于窗口的方法设计FIR滤波,默认是hamming窗),滤波器阶数设置为28。fadtool的配置如下:
37.13.png
配置好带通滤波器后,具体滤波器系数的生成大家参考本章第二小节的方法即可。

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     /* 滤波器系数个数 */

  4. uint32_t blockSize = BLOCK_SIZE;
  5. uint32_t numBlocks = TEST_LENGTH_SAMPLES/BLOCK_SIZE;            /* 需要调用arm_fir_f32的次数 */

  6. static float32_t testInput_f32_50Hz_200Hz[TEST_LENGTH_SAMPLES]; /* 采样点 */
  7. static float32_t testOutput[TEST_LENGTH_SAMPLES];               /* 滤波后的输出 */
  8. static float32_t firStateF32[BLOCK_SIZE + NUM_TAPS - 1];        /* 状态缓存,大小numTaps + blockSize - 1*/

  9. /* 带通滤波器系数 通过fadtool获取*/
  10. const float32_t firCoeffs32BP[NUM_TAPS] = {
  11. 0.003531039227f,    0.0002660876198f,   -0.001947779674f,  0.001266813371f,  -0.008019094355f,
  12. -0.01986379735f,    0.01018396299f,     0.03163734451f,    0.00165955862f,   0.03312643617f,
  13. 0.0622616075f,      -0.1229852438f,     -0.2399847955f,    0.07637182623f,   0.3482480049f,
  14. 0.07637182623f,     -0.2399847955f,     -0.1229852438f,    0.0622616075f,    0.03312643617f,
  15. 0.00165955862f,     0.03163734451f,     0.01018396299f,    -0.01986379735f,  -0.008019094355f,
  16. 0.001266813371f,   -0.001947779674f,    0.0002660876198f,  0.003531039227f
  17. };

  18. /*
  19. *********************************************************************************************************
  20. *        函 数 名: arm_fir_f32_bp
  21. *        功能说明: 调用函数arm_fir_f32_bp实现带通滤波器
  22. *        形    参:无
  23. *        返 回 值: 无
  24. *********************************************************************************************************
  25. */
  26. static void arm_fir_f32_bp(void)
  27. {
  28. uint32_t i;
  29. arm_fir_instance_f32 S;
  30. float32_t  *inputF32, *outputF32;

  31. /* 初始化输入输出缓存指针 */
  32. inputF32 = &testInput_f32_50Hz_200Hz[0];
  33. outputF32 = &testOutput[0];

  34. /* 初始化结构体S */
  35. arm_fir_init_f32(&S, NUM_TAPS, (float32_t *)&firCoeffs32BP[0], &firStateF32[0], blockSize);

  36. /* 实现FIR滤波 */
  37. for(i=0; i < numBlocks; i++)
  38. {
  39. arm_fir_f32(&S, inputF32 + (i * blockSize), outputF32 + (i * blockSize), blockSize);
  40. }

  41. /* 打印滤波后结果 */
  42. for(i=0; i<TEST_LENGTH_SAMPLES; i++)
  43. {
  44. printf("%f\r\n", testOutput[i]);
  45. }
  46. }
复制代码

运行如上函数可以通过串口打印出函数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;                  %频率序列

  9. x=sin(2*pi*50*t)+sin(2*pi*200*t);  %50Hz和200Hz正弦波混合           
  10. b=fir1(28, [125/500 300/500]);     %获得滤波器系数,截止频率125Hz和300Hz,带通滤波。
  11. y=filter(b, 1, x);                   %获得滤波后的波形
  12. subplot(211);
  13. plot(t, y);
  14. title('Matlab FIR滤波后的实际波形');
  15. grid on;

  16. subplot(212);
  17. plot(t, sampledata);        %绘制ARM官方库滤波后的波形。
  18. title('ARM官方库滤波后的实际波形');
  19. grid on;
复制代码
Matlab显示效果如下:
37.14.png
从上面的波形对比来看,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;                  %频率序列

  9. x=sin(2*pi*50*t)+sin(2*pi*200*t);  %50Hz和200Hz正弦波混合               
  10. subplot(221);
  11. plot(t, x);       %绘制信号x的波形                                                
  12. xlabel('时间');
  13. ylabel('幅值');
  14. title('原始信号');
  15. grid on;
  16.   
  17. subplot(222);
  18. y=fft(x, N);     %对信号x做FFT   
  19. plot(f,abs(y));
  20. xlabel('频率/Hz');
  21. ylabel('振幅');
  22. title('原始信号FFT');
  23. grid on;

  24. y3=fft(sampledata, N);       %经过FIR滤波器后得到的信号做FFT
  25. subplot(223);                              
  26. plot(f,abs(y3));
  27. xlabel('频率/Hz');
  28. ylabel('振幅');
  29. title('滤波后信号FFT');
  30. grid on;

  31. b=fir1(28, [125/500 300/500]);   %获得滤波器系数,截止频率125Hz,高通滤波。
  32. [H,F]=freqz(b,1,160);            %通过fir1设计的FIR系统的频率响应
  33. subplot(224);
  34. plot(F/pi,abs(H));                %绘制幅频响应
  35. xlabel('归一化频率');        
  36. title(['Order=',int2str(28)]);
  37. grid on;
复制代码
37.15.png
上面波形变换前的FFT和变换后FFT可以看出,50Hz的正弦波基本被滤除。

baiyongbin2009 回答时间:2015-4-24 12:36:16
37.6 FIR带阻滤波器设计
37.6.1 fdatool获取带阻滤波器系数
    设计一个如下的例子:
    信号由50Hz正弦波和200Hz正弦波组成,采样率1Kbps,现设计一个带阻滤波器,截止频率125Hz和300Hz,采样320个数据,采用函数fir1进行设计(注意这个函数是基于窗口的方法设计FIR滤波,默认是hamming窗),滤波器阶数设置为28。fadtool的配置如下:
37.16.png
配置好带阻滤波器后,具体滤波器系数的生成大家参考本章第二小节的方法即可。

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     /* 滤波器系数个数 */

  4. uint32_t blockSize = BLOCK_SIZE;
  5. uint32_t numBlocks = TEST_LENGTH_SAMPLES/BLOCK_SIZE;            /* 需要调用arm_fir_f32的次数 */

  6. static float32_t testInput_f32_50Hz_200Hz[TEST_LENGTH_SAMPLES]; /* 采样点 */
  7. static float32_t testOutput[TEST_LENGTH_SAMPLES];               /* 滤波后的输出 */
  8. static float32_t firStateF32[BLOCK_SIZE + NUM_TAPS - 1];        /* 状态缓存,大小numTaps + blockSize - 1*/

  9. /* 带阻滤波器系数 通过fadtool获取*/
  10. const float32_t firCoeffs32BPCheb[NUM_TAPS] = {
  11. 0.01801843569f,    0.0007182828849f,  -0.004868913442f,  0.002710500965f,  -0.01462193858f,
  12. -0.03147283196f,   0.01435638033f,    0.04055848345f,    0.00197162549f,   0.03706155345f,
  13. 0.06650412083f,    -0.1269270927f,    -0.2418768406f,    0.07591249049f,   0.3445736468f,
  14. 0.07591249049f,    -0.2418768406f,    -0.1269270927f,    0.06650412083f,   0.03706155345f,
  15. 0.00197162549f,    0.04055848345f,    0.01435638033f,    -0.03147283196f,  -0.01462193858f,
  16. 0.002710500965f,  -0.004868913442f,   0.0007182828849f,  0.01801843569f
  17. };*/
  18. const float32_t firCoeffs32BS[NUM_TAPS] = {
  19. -0.003560454352f,  -0.0002683042258f,  0.001964005642f,   -0.001277366537f,   0.008085897192f,
  20. 0.02002927102f,    -0.01026879996f,    -0.03190089762f,   -0.001673383522f,   -0.0334023945f,
  21. -0.06278027594f,   0.1240097657f,      0.2419839799f,     -0.07700803876f,    0.6521340013f,
  22. -0.07700803876f,   0.2419839799f,      0.1240097657f,     -0.06278027594f,    -0.0334023945f,
  23. -0.001673383522f,  -0.03190089762f,    -0.01026879996f,   0.02002927102f,     0.008085897192f,
  24. -0.001277366537f,  0.001964005642f,    -0.0002683042258f, -0.003560454352f
  25. };

  26. /*
  27. *********************************************************************************************************
  28. *        函 数 名: arm_fir_f32_bs
  29. *        功能说明: 调用函数arm_fir_f32_bs实现带阻滤波器
  30. *        形    参:无
  31. *        返 回 值: 无
  32. *********************************************************************************************************
  33. */
  34. static void arm_fir_f32_bs(void)
  35. {
  36. uint32_t i;
  37. arm_fir_instance_f32 S;
  38. float32_t  *inputF32, *outputF32;

  39. /* 初始化输入输出缓存指针 */
  40. inputF32 = &testInput_f32_50Hz_200Hz[0];
  41. outputF32 = &testOutput[0];

  42. /* 初始化结构体S */
  43. arm_fir_init_f32(&S, NUM_TAPS, (float32_t *)&firCoeffs32BS[0], &firStateF32[0], blockSize);

  44. /* 实现FIR滤波 */
  45. for(i=0; i < numBlocks; i++)
  46. {
  47. arm_fir_f32(&S, inputF32 + (i * blockSize), outputF32 + (i * blockSize), blockSize);
  48. }

  49. /* 打印滤波后结果 */
  50. for(i=0; i<TEST_LENGTH_SAMPLES; i++)
  51. {
  52. printf("%f\r\n", 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;                  %频率序列

  9. x=sin(2*pi*50*t)+sin(2*pi*200*t);       %50Hz和200Hz正弦波混合           
  10. b=fir1(28, [125/500 300/500], 'stop');   %获得滤波器系数,截止频率125Hz和300,带阻滤波。
  11. y=filter(b, 1, x);                        %获得滤波后的波形
  12. subplot(211);
  13. plot(t, y);
  14. title('Matlab FIR滤波后的实际波形');
  15. grid on;

  16. subplot(212);
  17. plot(t, sampledata);        %绘制ARM官方库滤波后的波形。
  18. title('ARM官方库滤波后的实际波形');
  19. grid on;
复制代码
Matlab运行结果如下:
37.17.png
从上面的波形对比来看,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;                  %频率序列

  9. x=sin(2*pi*50*t)+sin(2*pi*200*t);  %50Hz和200Hz正弦波混合           
  10. subplot(221);
  11. plot(t, x);   %绘制信号x的波形                                                
  12. xlabel('时间');
  13. ylabel('幅值');
  14. title('原始信号');
  15. grid on;
  16.   
  17. subplot(222);
  18. y=fft(x, N);     %对信号x做FFT   
  19. plot(f,abs(y));
  20. xlabel('频率/Hz');
  21. ylabel('振幅');
  22. title('原始信号FFT');
  23. grid on;

  24. y3=fft(sampledata, N);       %经过FIR滤波器后得到的信号做FFT
  25. subplot(223);                              
  26. plot(f,abs(y3));
  27. xlabel('频率/Hz');
  28. ylabel('振幅');
  29. title('滤波后信号FFT');
  30. grid on;

  31. b=fir1(28, [125/500 300/500], 'stop');  %获得滤波器系数,截止频率125Hz和300Hz,带阻滤波。     
  32. [H,F]=freqz(b,1,160);                  %通过fir1设计的FIR系统的频率响应
  33. subplot(224);
  34. plot(F/pi,abs(H));             %绘制幅频响应
  35. xlabel('归一化频率');        
  36. title(['Order=',int2str(28)]);
  37. grid on;
复制代码
Matlab运行效果如下:
37.18.png
上面波形变换前的FFT和变换后FFT可以看出,200Hz的正弦波基本被滤除。

baiyongbin2009 回答时间:2015-4-24 12:41:00
37.7 切比雪夫窗口设计带通滤波器
37.7.1 fdatool获取滤波器系数
    设计一个如下的例子:
    信号由50Hz正弦波和200Hz正弦波组成,采样率1Kbps,现设计一个使用切比雪夫窗口的带通滤波器,截止频率125Hz和300Hz,切比雪夫波纹设置为30db,采样320个数据,采用函数fir1进行设计(注意这个函数是基于窗口的方法设计FIR滤波,默认是hamming窗),滤波器阶数设置为28。fadtool的配置如下:
37.19.png
配置好带通滤波器后,具体滤波器系数的生成大家参考本章第二小节的方法即可。

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     /* 滤波器系数个数 */

  4. uint32_t blockSize = BLOCK_SIZE;
  5. uint32_t numBlocks = TEST_LENGTH_SAMPLES/BLOCK_SIZE;            /* 需要调用arm_fir_f32的次数 */

  6. static float32_t testInput_f32_50Hz_200Hz[TEST_LENGTH_SAMPLES]; /* 采样点 */
  7. static float32_t testOutput[TEST_LENGTH_SAMPLES];               /* 滤波后的输出 */
  8. static float32_t firStateF32[BLOCK_SIZE + NUM_TAPS - 1];        /* 状态缓存,大小numTaps + blockSize - 1*/
  9. /* 带通滤波器系数 切比雪夫窗口 通过fadtool获取*/
  10. const float32_t firCoeffs32BPCheb[NUM_TAPS] = {
  11. 0.01801843569f,    0.0007182828849f,  -0.004868913442f,  0.002710500965f,  -0.01462193858f,
  12. -0.03147283196f,   0.01435638033f,    0.04055848345f,    0.00197162549f,   0.03706155345f,
  13. 0.06650412083f,    -0.1269270927f,    -0.2418768406f,    0.07591249049f,   0.3445736468f,
  14. 0.07591249049f,    -0.2418768406f,    -0.1269270927f,    0.06650412083f,   0.03706155345f,
  15. 0.00197162549f,    0.04055848345f,    0.01435638033f,    -0.03147283196f,  -0.01462193858f,
  16. 0.002710500965f,  -0.004868913442f,   0.0007182828849f,  0.01801843569f
  17. };

  18. /*
  19. *********************************************************************************************************
  20. *        函 数 名: arm_fir_f32_bp
  21. *        功能说明: 调用函数arm_fir_f32_bp实现带通滤波器
  22. *        形    参:无
  23. *        返 回 值: 无
  24. *********************************************************************************************************
  25. */
  26. static void arm_fir_f32_bp(void)
  27. {
  28. uint32_t i;
  29. arm_fir_instance_f32 S;
  30. float32_t  *inputF32, *outputF32;

  31. /* 初始化输入输出缓存指针 */
  32. inputF32 = &testInput_f32_50Hz_200Hz[0];
  33. outputF32 = &testOutput[0];

  34. /* 初始化结构体S */
  35. arm_fir_init_f32(&S, NUM_TAPS, (float32_t *)&firCoeffs32BP[0], &firStateF32[0], blockSize);

  36. /* 实现FIR滤波 */
  37. for(i=0; i < numBlocks; i++)
  38. {
  39. arm_fir_f32(&S, inputF32 + (i * blockSize), outputF32 + (i * blockSize), blockSize);
  40. }

  41. /* 打印滤波后结果 */
  42. for(i=0; i<TEST_LENGTH_SAMPLES; i++)
  43. {
  44. printf("%f\r\n", testOutput[i]);
  45. }
  46. }
复制代码
运行如上函数可以通过串口打印出函数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;                 %频率序列

  9. x=sin(2*pi*50*t)+sin(2*pi*200*t);          %50Hz和200Hz正弦波混合      
  10. Window = chebwin(29, 30);               %30db的切比雪夫窗
  11. b=fir1(28, [125/500 300/500], Window);   %获得滤波器系数,截止频率125Hz和300Hz,带通滤波。
  12. y=filter(b, 1, x);                           %获得滤波后的波形
  13. subplot(211);
  14. plot(t, y);
  15. title('Matlab FIR滤波后的实际波形');
  16. grid on;

  17. subplot(212);
  18. plot(t, sampledata);                 %绘制ARM官方库滤波后的波形。
  19. title('ARM官方库滤波后的实际波形');
  20. grid on;
复制代码
Matlab的运行效果如下:
37.20.png
从上面的波形对比来看,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;                  %频率序列


  9. x=sin(2*pi*50*t)+sin(2*pi*200*t);  %50Hz和200Hz正弦波混合               
  10. subplot(221);
  11. plot(t, x);   %绘制信号x的波形                                                
  12. xlabel('时间');
  13. ylabel('幅值');
  14. title('原始信号');
  15. grid on;
  16.   
  17. subplot(222);
  18. y=fft(x, N);     %对信号x做FFT   
  19. plot(f,abs(y));
  20. xlabel('频率/Hz');
  21. ylabel('振幅');
  22. title('原始信号FFT');
  23. grid on;

  24. y3=fft(sampledata, N);       %经过FIR滤波器后得到的信号做FFT
  25. subplot(223);                              
  26. plot(f,abs(y3));
  27. xlabel('频率/Hz');
  28. ylabel('振幅');
  29. title('滤波后信号FFT');
  30. grid on;

  31. Window = chebwin(29, 30);       %30db的切比雪夫窗
  32. b=fir1(28, [125/500 300/500], Window);  %获得滤波器系数,截止频率125Hz和300Hz,带通滤波。
  33. [H,F]=freqz(b,1,160);                    %通过fir1设计的FIR系统的频率响应
  34. subplot(224);
  35. plot(F/pi,abs(H));                       %绘制幅频响应
  36. xlabel('归一化频率');        
  37. title(['Order=',int2str(28)]);
  38. grid on;
复制代码
Matlab运行结果如下:
37.21.png
上面波形变换前的FFT和变换后FFT可以看出,50Hz的正弦波基本被滤除,在归一化频率中我们可以看到一定的波纹

baiyongbin2009 回答时间:2015-4-24 12:44:24
37.8 FIR滤波后的群延迟
    波形经过FIR滤波器后,输出的波形会有一定的延迟。对于线性相位的FIR,这个群延迟就是一个常数。但是实际应用中这个群延迟是多少呢,关于群延迟的数值,fdatool工具箱会根据用户的配置计算好。
    比如咱们前面设计的28阶FIR高通,低通,带通和带阻滤波器的群延迟就是14,反应在实际的采样值上就是滤波后输出数据的第15个才是实际滤波后的波形数据起始点。
    下面是看群延迟采样点的位置:
37.22.png
细心的读者可能发现全面做低通,高通,带通和带阻滤波后,输出的波形前面几个点感觉有问题,其实就是群延迟造成的。
    为了更好的说明这个问题,下面再使用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;                  %频率序列

  6. x1=sin(2*pi*50*t);
  7. x2=sin(2*pi*200*t);
  8. x=sin(2*pi*50*t)+sin(2*pi*200*t);  %50Hz和200Hz正弦波混合               

  9. plot(n, x1, 'b');   %绘制信号x的波形  
  10. xlabel('时间');
  11. ylabel('幅值');
  12. title('原始信号和滤波后信号');
  13. hold on;
  14.   
  15. b=fir1(28, 125/500);     %获得滤波器系数,截止频率125Hz.
  16. y=filter(b, 1, x);
  17. plot(n(1:305), y(15:319), 'r');
  18. legend('原始信号','滤波后信号');
  19. grid on;
复制代码
Matlab的运行结果如下:
37.23.png
可以看出,显示波形基本重合,这个说明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;                  %频率序列

  6. x1=sin(2*pi*50*t);
  7. x2=sin(2*pi*200*t);
  8. x=sin(2*pi*50*t)+sin(2*pi*200*t);  %50Hz和200Hz正弦波混合               

  9. plot(n, x2, 'b');   %绘制信号x的波形  
  10. xlabel('时间');
  11. ylabel('幅值');
  12. title('原始信号和滤波后信号');
  13. hold on;
  14.   
  15. b=fir1(28, 125/500, 'high');     %获得滤波器系数,截止频率125Hz.
  16. y=filter(b, 1, x);
  17. plot(n(1:305), y(15:319), 'r');
  18. legend('原始信号','滤波后信号');
  19. grid on;
复制代码
Matlab运行结果如下:
37.24.png
可以看出,显示波形基本重合,这个说明14个采样点的群延迟也是是正确的。大家在使用FIR滤波器的时候一定要注意这个问题。

baiyongbin2009 回答时间:2015-4-24 12:44:52
37.9  总结
    本章节主要讲解了FIR滤波器的低通,高通,带通和带阻滤波器的实现,同时一定要注意线性相位FIR滤波器的群延迟问题。


12下一页

所属标签

STM32团队

意法半导体微控制器和微处理器拥有广泛的产品线,包含低成本的8位单片机和基于ARM® Cortex®-M0、M0+、M3、M4、M33、M7及A7内核并具备丰富外设选择的32位微控制器及微处理器


最新内容

关于
我们是谁
投资者关系
意法半导体可持续发展举措
创新与技术
意法半导体官网
联系我们
联系ST分支机构
寻找销售人员和分销渠道
社区
媒体中心
活动与培训
隐私策略
隐私策略
Cookies管理
行使您的权利
官方最新发布
STM32N6 AI生态系统
STM32MCU,MPU高性能GUI
ST ACEPACK电源模块
意法半导体生物传感器
STM32Cube扩展软件包
关注我们
微信公众号二维码 微信公众号
手机版二维码 手机版