FIR滤波器实践-基于ARM CMSIS-DSP库

FIR滤波器实践-基于ARM CMSIS-DSP库

一、CMSIS-DSP库介绍

    CMSIS-DSP库是ARM官方提供的数字信号处理库，包含了大量常用的数字信号处理算法，如快速傅里叶变换（FFT）、FIR滤波、自适应滤波器等等。由于是ARM官方提供，因此可以适配Cortex-A、Cortex-R、Cortex-M三种内核的单片机，典型的如Cortex-M内核的STM32单片机。CMSIS-DSP通常在安装Keil的时候会自动安装在keil根目录下的ARM文件夹下，采用不同的软件有不同的库包含方法，可以去百度看看。库的包含非常简单，这里也不做介绍。

二、Q15格式介绍

    Q15是最常见的浮点转定点的数据格式，Q15数据格式根据定义是长度为16位的数据，其中最高位为符号位，后面15位为小数位。Q15格式下的0x1111(十进制为4369),二级制数据为0001 0001 0001 0001，根据Q15的定义，最高位为符号位，0代表正数，后续位为小数位，计算过程如下
$$1\ast 2^{-3}+1\ast 2^{-7}+2^{-11}+2^{-15}=0.1333312988$$
    Q15数据格式本质是将小数左移15位，即扩大$2^{15}$变成一个整数进行运算，0x1111为十进制数4369，Q15对应下的小数为$4369/2^{15}=0.1333312988$ 。对于没有浮点运算单元(FPU)的单片机来说，将小数转换成浮点数可以显著提高计算效率。Q15格式下所能表示的整数数据范围为-32768 ~ 32767 ，对应的小数范围为-1 ~ 0.9999694824 。下面再举一个实际的例子，例如一个12bit的ADC采集到的电压数据为2369，其中参考电压为3.3V，则实际的电压为$2369/4096\times 3.3=1.908618164$V，现在需要对该电压缩小0.4倍，为了避免浮点运算，可以采用如下的计算方式
$$\mathrm{int}((2369\ast Q15(0.4))>>15)$$

    其中，$Q15(0.4)=0.4\ast 2^{15}=13107.2$，为了避免浮点运算，取整后为13107。由此，上式可以转换为

$$\mathrm{int}((2369\ast 13107)>>15)=948$$

    其中，int取整操作为单片机运算定点数时自动完成，现在可以计算一下948对应的真实电压值，$948/4096\times 3.3=0.7637695313$V，真实电压的0.4倍为0.7634472656V 。可以看到，尽管定点后的数据精度有所丧失，但在实际工程允许的范围内，使用是没有任何问题的，毕竟单片机没有浮点运算单元，牺牲精度提升速度有时是必须的。

    至此，我们对Q15格式已经有了一个初步的认识，就是将一个范围为-1 ~ 0.9999694824 范围内的小数转换成整数进行运算，当我们需要乘以一个固定的小数时，先将这个小数扩大$2^{15}$，就变成了Q15格式的整数，然后乘以一个数，再右移15位(即除以$2^{15}$)，这样就避免了小数运算，中间的精度损失出现在整数运算的截断中。Q15在实际运算的时候采用的int16格式的定义，int16最高位为符号位，表示的数据范围为-32768~32767，但是与传统int16格式不同的是，我们将Q15格式下的符号位后的数据解释成小数，本质还是int16格式，只是人为解释一番。下面举一个实际的例子，在CMSIS-DSP库中有Q15格式下的正弦波，调用如下所示

/**
  @brief         Fast approximation to the trigonometric sine function for Q15 data.
  @param[in]     x  Scaled input value in radians
  @return        sin(x)

  The Q15 input value is in the range [0 +0.9999] and is mapped to a radian value in the range [0 2*PI).
*/
q15_t arm_sin_q15(q15_t x)

    从代码注释中可以知道，输入的变量x为Q15格式，输入的范围为0 ~ 32767 ，对应小数范围为0 ~ 0.9999 ，对应的弧度值为0~2$\pi$ 。假设我们需要生成一个采样频率$f$=1000Hz，信号频率为50Hz的采样序列，正弦信号的采样模型为
$$y[n]=A\ast sin(2\ast \pi \ast f_0/f_s\ast n)n=0,1,\mathrm{2...}N-1$$
    在使用软件生成正弦序列时，$n$即为for循环中的变量，我们观察如下变换
$$\begin{array}{rlr}& sin(\theta )=sin(2\ast \pi \ast f_0/f_s\ast n)& \\ & \theta =2\ast \pi \ast \frac{f_0}{f_s}\ast n& \\ & \frac{\theta }{2\pi }=\frac{f_0}{f_s}\ast n& \\ & \frac{\theta }{2\pi }\ast 2^{15}=\frac{f_0}{f_s}\ast n\ast 2^{15}& \end{array}\text{\hspace{1em}(1)(2)(3)(4)}$$
    上式演示了Q15格式的来源，在式(3)中，将等式左边除以2$\pi$后，范围变成了[0,1]，再左移15位，范围就变成了[0,32767]。对于上述问题
$$\frac{f_0}{f_s}=50/1000=0.05$$
    Q15后的数据为$0.05\ast 2^{15}=1638.4$ ，去掉小数位，参数为1638，信号生成代码如下：

q15_t Sin_50Hz[256];
void Gen_Sin_50Hz(void)
{
	uint16_t i = 0;
	for(i=0; i<256; i++)
	{
		/* 50Hz正弦波 采样率1KHz */
		Sin_50Hz[i] = arm_sin_q15(1638*i);
	}
}

    在实时仿真环境中点击Keil Debug选项下的Function Editor(open ini file)，在.uvprojx工程所在的目录下面新建一个txt文件，并将txt文件改成后缀为ini的文件，文件名字随意。在Function Editor中打开，并编写如下代码如下图所示：

Function Editor的配置

    点击Complier，观察左下角Command窗口是否有错误，无错误则将函数名复制，将my_print()放置到Command窗口中，与C语言运行一个函数是一样的。运行后可以看到窗口数据已经打印出来了，将打印的数组数据复制。打开Matlab软件，在工作区中右击新建一个变量，双击变量后进入变量查看窗口，选中第一列，黏贴数据，如下图所示：

Matlab变量

    接着选中第一列，在绘图窗口中点击第一个Plot选项，如下图所示：

50Hz采样波形

    可以通过过零点的位置简单验证信号的频率是否正确，1kHz的采样率下，50Hz信号一个周期中会有20个采样点，因此可以查看两个过零点或者峰值点之间的间隔是否是20个点来判断信号频率是否为50Hz。另外，我们可以发现该函数输出的幅值范围为-32768 ~ 32767之间，对应Q15格式为 -1 ~ 0.999 之间 。

    至此，我们已经完成了信号波形的生成，为了完成后续滤波器实验，我们还需要生成一个200Hz的信号，代码如下：

q15_t Sin_50Hz_200Hz[256];
void Gen_Sin_50Hz_200Hz(void)
{
	uint16_t i = 0;
	for(i=0; i<256; i++)
	{
		/* 50Hz正弦波 + 200Hz正弦波，采样率1KHz */
		Sin_50Hz_200Hz[i] = (arm_sin_q15(1638*i)>>2) + (arm_sin_q15(6554*i)>>2);
	}
}

    其中>>2是为了保证相加后数据范围不超过Q15的范围。混合后的信号波形如下：

50Hz+200Hz采样波形

三、数字滤波器的设计

    本文偏向实践，因此不对原理进行过多介绍，只介绍如何使用Matlab工具与CMSIS-DSP库快速在单片机上面实现FIR滤波器。首先打开Matlab，在命令行中输入filterDesigner，可以打开如下图所示的界面。

Filter Designer界面

    从左往后看，首先在响应类型中选择低通，在设计方法中选择FIR，下拉框中可以随意选择，本文选择窗函数法。滤波器阶数中选择指定阶，设置为20。在窗选项中选择Hamming窗，频率Fs为采样频率，设置为1000Hz，Fc为截止频率，设置为100Hz，最后点击设计滤波器。最终如下图所示：

Filter Designer最终参数界面

    点击目标 ，生成C头文件，导出设置为有符号16位整数，如下图所示：

系数导出界面

    在导出的.h文件中有如下代码，默认为int16_T格式，我们将其改成q15_t即可

const int BL = 21;
const q15_t filtercoff[21] = {
        0,   -113,   -129,    237,    659,      0,  -1695,  -1659,   2802,
     9716,  13134,   9716,   2802,  -1659,  -1695,      0,    659,    237,
     -129,   -113,      0
};

四、CMSIS-DSP库中的arm_fir_q15函数介绍

    在DSP库中，实现Q15格式下的FIR滤波的函数为arm_fir_q15，该函数的声明如下

  /**
   * @brief Processing function for the Q15 FIR filter.
   * @param[in]  S          points to an instance of the Q15 FIR structure.
   * @param[in]  pSrc       points to the block of input data.
   * @param[out] pDst       points to the block of output data.
   * @param[in]  blockSize  number of samples to process.
   */
  void arm_fir_q15(
  const arm_fir_instance_q15 * S,
  const q15_t * pSrc,
        q15_t * pDst,
        uint32_t blockSize);

    函数的入口参数需要传递一个arm_fir_instance_q15类型的结构体、输入数据、输出输出、块运算的大小。块运算的大小其实就是一次运算多少个数据，为1则来一个数据运算一次，即实时滤波；为2就是每两个数据运算一次滤波函数，一般为1。arm_fir_instance_q15的定义如下：

  typedef struct
  {
          uint16_t numTaps;         /**< number of filter coefficients in the filter. */
          q15_t *pState;            /**< points to the state variable array. The array is of length numTaps+blockSize-1. */
    const q15_t *pCoeffs;           /**< points to the coefficient array. The array is of length numTaps.*/
  } arm_fir_instance_q15;

    arm_fir_instance_q15由三个部分组成，第一个为抽头系数，即滤波器的系数长度，本文为21；第二个参数为中间状态缓存区，大小为numTaps+blockSize-1，即21+1-1=20；第三个参数为滤波器的系数数组。对于该结构体的初始化可以采用直接赋值，也可以采用DSP库中的函数arm_fir_init_q15，初始化代码如下:

const int BL = 21;
const q15_t filtercoff[21] = {
        0,   -113,   -129,    237,    659,      0,  -1695,  -1659,   2802,
     9716,  13134,   9716,   2802,  -1659,  -1695,      0,    659,    237,
     -129,   -113,      0
};
q15_t firstate[20]; 
void arm_fir_lp(void)
{
	arm_fir_instance_q15 S;
	arm_fir_init_q15(&S,BL,&filtercoff[0],&firstate[0],1);
}

    为了测试滤波函数，我们还需要定义输出数组，完整代码如下所示：

q15_t Sin_50Hz_200Hz[256];
q15_t Fir_Out[256]; 
const int BL = 21;
const q15_t filtercoff[21] = {
        0,   -113,   -129,    237,    659,      0,  -1695,  -1659,   2802,
     9716,  13134,   9716,   2802,  -1659,  -1695,      0,    659,    237,
     -129,   -113,      0
};
q15_t firstate[20]; 
void Gen_Sin_50Hz_200Hz(void)
{
	uint16_t i = 0;
	for(i=0; i<256; i++)
	{
		/* 50Hz正弦波 + 200Hz正弦波，采样率1KHz */
		Sin_50Hz_200Hz[i] = (arm_sin_q15(1638*i)>>2) + (arm_sin_q15(6554*i)>>2);
	}
}
void arm_fir_lp(void)
{
	uint16_t i = 0;
	arm_fir_instance_q15 S;
	arm_fir_init_q15(&S,BL,&filtercoff[0],&firstate[0],1);
	for(i=0;i<256;i++)
	{
		arm_fir_q15(&S,&Sin_50Hz_200Hz[i],&Fir_Out[i],1);
	}
	
}

    最终滤波后的波形如下：

滤波后的波形

    使用Matlab里面自带的Signal Analyzer对比时域和频域如下：

原始波形与滤波后波形对比（50Hz+200Hz）

    可以看到，经过滤波器之后，原始信号中的200Hz分量被明显抑制，由于截止频率设置为100Hz，因此200Hz信号的频率并没有被抑制得非常明显，但是本实验也实现了Q15格式下FIR滤波器设计。对于混合400Hz信号分量的实验结果如下图所示：

原始波形与滤波后波形对比（50Hz+400Hz）

小结

    本实验基于CMSIS-DSP库实现了Q15格式下的FIR滤波器设计，并通过硬件在线仿真验证了FIR滤波器的有效性。本文未去对比采用浮点与定点两种滤波器频率响应的差距，由于笔者手中单片机没有浮点运算单元，因此没有进行对比实验。对于浮点数来说，不需要考虑正负，所有运算均采用真实值，因此设计过程非常简单。