IIR无限冲击响应滤波器

概述

实际项目中，常常需要从叠加频率的信号中获取特定频带的数据，这时候就需要使用到低通、高通、带通或陷波器来进行数字信号处理。而IIR滤波器具有结构简单、运算量小等特点，非常适用于嵌入式系统的数字信号前端处理。所以，这里分享如何使用arm_dsp库和matlab设计简单的2阶IIR直接I型滤波器。

正文

matlab设计滤波器系数矩阵

• matlab输入fdatool调出滤波器设计工具。
• 输入要设计的滤波器类型及各项参数，如下。
  
• 点击Edit->Convert Structure，选择是 Direct-Form I。
• 转换好以后再点击 File-Export，第一项选择 Coefficient File（ASCII）,第二项选择Decimal。
• 最后，点击Export即可输出系数。
• 打开输出的fcf文件，根据系数，将a0系数删除，将a1、a2系数取反（使用ARM提供给的DSP库）

ARM DSP库IIR滤波器API

这里提供自己整理的驱动文件作为参考，滤波效果还算可以。

#include "iir.h"

/**********************************USER TODO*************************************************/
#define NUMSTAGES  			2 		//2阶iir滤波器个数
#define BLOCK_SIZE          1		//数据处理块大小

//巴特沃斯低通滤波器系数 80Hz截止频率 250Hz采样率
//a1、a2要取反，例如：matlab为1.001，则应用时要为-1.001
static const float32_t IIRCoeffs32LP[5*NUMSTAGES] = 
{                                                                                 
	1.0f,  2.0f,  1.0f,  -0.63253540497219674f,  -0.48559457329907518f,        
	1.0f,  2.0f,  1.0f,  -0.46382419413091569f,  -0.089353576652349942f,                                                                                              
}; 

static const float32_t IIRScaleValue[NUMSTAGES] =
{
	0.52953249456781792f,                                        
	0.38829444269581642f,  
};

/********************************************************************************************/
static float32_t IIRStateF32[4*NUMSTAGES];
static arm_biquad_casd_df1_inst_f32 S;

/*
@brief         Initialization function for the floating-point Biquad cascade filter.
@param[in,out] S           points to an instance of the floating-point Biquad cascade structure.
@param[in]     numStages   number of 2nd order stages in the filter.
@param[in]     pCoeffs     points to the filter coefficients.
@param[in]     pState      points to the state buffer.
@return        none
*/
static void arm_biquad_cascade_df1_init_f32(arm_biquad_casd_df1_inst_f32 * S, uint8_t numStages, const float32_t * pCoeffs, float32_t * pState)
{
	/* Assign filter stages */
	S->numStages = numStages;

	/* Assign coefficient pointer */
	S->pCoeffs = pCoeffs;

	/* Clear state buffer and size is always 4 * numStages */
	memset(pState, 0, (4U * (uint32_t) numStages) * sizeof(float32_t));

	/* Assign state pointer */
	S->pState = pState;
}

/*
@brief		 Processing function for the floating-point Biquad cascade filter.
@param[in]	 S		   points to an instance of the floating-point Biquad cascade 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
@return		 none
*/
static void arm_biquad_cascade_df1_f32(const arm_biquad_casd_df1_inst_f32 * S, const float32_t * pSrc, float32_t * pDst, uint32_t blockSize)
{
  const float32_t *pIn = pSrc;						   /* Source pointer */
		float32_t *pOut = pDst; 					   /* Destination pointer */
		float32_t *pState = S->pState;				   /* pState pointer */
  const float32_t *pCoeffs = S->pCoeffs;			   /* Coefficient pointer */
		float32_t acc;								   /* Accumulator */
		float32_t b0, b1, b2, a1, a2;				   /* Filter coefficients */
		float32_t Xn1, Xn2, Yn1, Yn2;				   /* Filter pState variables */
		float32_t Xn;								   /* Temporary input */
		uint32_t sample, stage = S->numStages;		   /* Loop counters */

  do
  {
	/* Reading the coefficients */
	b0 = *pCoeffs++;
	b1 = *pCoeffs++;
	b2 = *pCoeffs++;
	a1 = *pCoeffs++;
	a2 = *pCoeffs++;

	/* Reading the pState values */
	Xn1 = pState[0];
	Xn2 = pState[1];
	Yn1 = pState[2];
	Yn2 = pState[3];

	/* Apply loop unrolling and compute 4 output values simultaneously. */
	/* Variable acc hold output values that are being computed:
	 *
	 * acc =  b0 * x[n] + b1 * x[n-1] + b2 * x[n-2] + a1 * y[n-1] + a2 * y[n-2]
	 * acc =  b0 * x[n] + b1 * x[n-1] + b2 * x[n-2] + a1 * y[n-1] + a2 * y[n-2]
	 * acc =  b0 * x[n] + b1 * x[n-1] + b2 * x[n-2] + a1 * y[n-1] + a2 * y[n-2]
	 * acc =  b0 * x[n] + b1 * x[n-1] + b2 * x[n-2] + a1 * y[n-1] + a2 * y[n-2]
	 */

	/* Loop unrolling: Compute 4 outputs at a time */
	sample = blockSize >> 2U;

	while (sample > 0U)
	{
	  /* Read the first input */
	  Xn = *pIn++;

	  /* acc =	b0 * x[n] + b1 * x[n-1] + b2 * x[n-2] + a1 * y[n-1] + a2 * y[n-2] */
	  Yn2 = (b0 * Xn) + (b1 * Xn1) + (b2 * Xn2) + (a1 * Yn1) + (a2 * Yn2);

	  /* Store output in destination buffer. */
	  *pOut++ = Yn2;

	  /* Every time after the output is computed state should be updated. */
	  /* The states should be updated as: */
	  /* Xn2 = Xn1 */
	  /* Xn1 = Xn  */
	  /* Yn2 = Yn1 */
	  /* Yn1 = acc */

	  /* Read the second input */
	  Xn2 = *pIn++;

	  /* acc =	b0 * x[n] + b1 * x[n-1] + b2 * x[n-2] + a1 * y[n-1] + a2 * y[n-2] */
	  Yn1 = (b0 * Xn2) + (b1 * Xn) + (b2 * Xn1) + (a1 * Yn2) + (a2 * Yn1);

	  /* Store output in destination buffer. */
	  *pOut++ = Yn1;

	  /* Every time after the output is computed state should be updated. */
	  /* The states should be updated as: */
	  /* Xn2 = Xn1 */
	  /* Xn1 = Xn  */
	  /* Yn2 = Yn1 */
	  /* Yn1 = acc */

	  /* Read the third input */
	  Xn1 = *pIn++;

	  /* acc =	b0 * x[n] + b1 * x[n-1] + b2 * x[n-2] + a1 * y[n-1] + a2 * y[n-2] */
	  Yn2 = (b0 * Xn1) + (b1 * Xn2) + (b2 * Xn) + (a1 * Yn1) + (a2 * Yn2);

	  /* Store output in destination buffer. */
	  *pOut++ = Yn2;

	  /* Every time after the output is computed state should be updated. */
	  /* The states should be updated as: */
	  /* Xn2 = Xn1 */
	  /* Xn1 = Xn  */
	  /* Yn2 = Yn1 */
	  /* Yn1 = acc */

	  /* Read the forth input */
	  Xn = *pIn++;

	  /* acc =	b0 * x[n] + b1 * x[n-1] + b2 * x[n-2] + a1 * y[n-1] + a2 * y[n-2] */
	  Yn1 = (b0 * Xn) + (b1 * Xn1) + (b2 * Xn2) + (a1 * Yn2) + (a2 * Yn1);

	  /* Store output in destination buffer. */
	  *pOut++ = Yn1;

	  /* Every time after the output is computed state should be updated. */
	  /* The states should be updated as: */
	  /* Xn2 = Xn1 */
	  /* Xn1 = Xn  */
	  /* Yn2 = Yn1 */
	  /* Yn1 = acc */
	  Xn2 = Xn1;
	  Xn1 = Xn;

	  /* decrement loop counter */
	  sample--;
	}

	/* Loop unrolling: Compute remaining outputs */
	sample = blockSize & 0x3U;

	while (sample > 0U)
	{
	  /* Read the input */
	  Xn = *pIn++;

	  /* acc =	b0 * x[n] + b1 * x[n-1] + b2 * x[n-2] + a1 * y[n-1] + a2 * y[n-2] */
	  acc = (b0 * Xn) + (b1 * Xn1) + (b2 * Xn2) + (a1 * Yn1) + (a2 * Yn2);

	  /* Store output in destination buffer. */
	  *pOut++ = acc;

	  /* Every time after the output is computed state should be updated. */
	  /* The states should be updated as: */
	  /* Xn2 = Xn1 */
	  /* Xn1 = Xn  */
	  /* Yn2 = Yn1 */
	  /* Yn1 = acc */
	  Xn2 = Xn1;
	  Xn1 = Xn;
	  Yn2 = Yn1;
	  Yn1 = acc;

	  /* decrement loop counter */
	  sample--;
	}

	/* Store the updated state variables back into the pState array */
	*pState++ = Xn1;
	*pState++ = Xn2;
	*pState++ = Yn1;
	*pState++ = Yn2;

	/* The first stage goes from the input buffer to the output buffer. */
	/* Subsequent numStages occur in-place in the output buffer */
	pIn = pDst;

	/* Reset output pointer */
	pOut = pDst;

	/* decrement loop counter */
	stage--;

  } while (stage > 0U);
}


void arm_iir_f32_init(void)
{
	arm_biquad_cascade_df1_init_f32(&S, NUMSTAGES, (float32_t *)&IIRCoeffs32LP[0], (float32_t *)&IIRStateF32[0]);
}

void arm_iir_f32_process(float32_t *inputF32, float32_t *outputF32)
{
	float32_t output[BLOCK_SIZE];
	arm_biquad_cascade_df1_f32(&S, inputF32,  output,  BLOCK_SIZE);
	for(uint8_t i=0;i<BLOCK_SIZE;i++)
		outputF32[i] = output[i]*IIRScaleValue[0]*IIRScaleValue[1];
}

#ifndef __IIR_H
#define __IIR_H

#include <stdio.h>
#include <stdint.h>
#include <string.h>

typedef float float32_t;

typedef struct
{
	uint32_t numStages;      		/**< number of 2nd order stages in the filter.  Overall order is 2*numStages. */
   	float32_t *pState;       		/**< Points to the array of state coefficients.  The array is of length 4*numStages. */
	const float32_t *pCoeffs;      	/**< Points to the array of coefficients.  The array is of length 5*numStages. */
} arm_biquad_casd_df1_inst_f32;

void arm_iir_f32_init(void);
void arm_iir_f32_process(float32_t *inputF32, float32_t *outputF32);
#endif

总结

高通、带通等其他类型的滤波器设计大体相同，重点是通过matlab输出系数矩阵即可。本文仅提供简单的滤波器设计，实际项目中需要根据需求进行调试优化。