DSP定点Q31 IIR滤波器的matlab设置，使用函数arm_biquad_cascade_df1_q31

eric2013 发表于 2022-8-24 19:42:54

测试工程：
https://www.armbbs.cn/static/image/filetype/zip.gif定点Q31 IIR滤波.7z (4.70MB)

static/image/hrline/4.gif

测试波形是50Hz和200Hz正弦波，采样率设置的1KHz，IIR低通的截止频率设置的80Hz。

Matlab设置如下（详细设置看DSP教程第2版的IIR章节即可，生成方法不变）：

将生成的参数代入到工程代码，滤波输出效果如下，未做定标处理

滤波效果是没问题的，完全正常（简单做了放缩处理，方便观察）

代码如下：
由于我们继续使用的浮点生成的系数，这个系统不能直接使用，需要先缩小，防止溢出，然后再调用函数arm_float_to_q31加大。

最后再配合函数arm_biquad_cascade_df1_init_q31做调整

#define numStages5             /* 5阶IIR滤波的个数 */
#define TEST_LENGTH_SAMPLES400 /* 采样点数 */
#define BLOCK_SIZE       400 /* 调用一次arm_biquad_cascade_df1_f32处理的采样点个数 */

uint32_t blockSize = BLOCK_SIZE;
uint32_t numBlocks = TEST_LENGTH_SAMPLES/BLOCK_SIZE; /* 需要调用arm_biquad_cascade_df1_f32的次数 */

/////////////////////////////////////// 定点Q31 /////////////////////////////////////////////////////////////////////////
static q31_t testInput_f32_50Hz_200Hz1; /* 采样点 */
static q31_t testOutput1;             /* 滤波后的输出 */
static q31_t IIRStateF321;                   /* 状态缓存 */

q31_t IIRCoeffs32LP1;

/* 巴特沃斯低通滤波器系数 80Hz*/
float32_t IIRCoeffs32LP = {
   1,2,1, 1.62978785624054456349085739930160343647, -0.85983730736705732411451208463404327631,
   1,2,1, 1.438087052475096605519411241402849555016,-0.641077359359069620658999610895989462733,
   1,2,1, 1.307285028849323404998017394973430782557,-0.491812237222575276973657310009002685547,
   1,2,1, 1.226250602643673470737439856748096644878,-0.39933955836362394764549321735103148967,
   1 , 2,1, 1.187550247801674174752406543120741844177,-0.355176532195590921592298627729178406298
};

/*
*********************************************************************************************************
*    函数名: arm_iir_q31_lp
*    功能说明: 调用函数arm_iir_q31_lp实现低通滤波器
*    形参：无
*    返回值: 无
*********************************************************************************************************
*/
static void arm_iir_q31_lp(void)
{
   uint32_t i;
   arm_biquad_casd_df1_inst_q31 S;
   q63_t ScaleValue;
   q31_t *inputQ31, *outputQ31;

   /* 使用浮点系数先进行一次缩放，然后再使用函数arm_float_to_q31转换为Q31，防止溢出 */
   for(int i = 0; i < 5*numStages; i++)
   {
            IIRCoeffs32LP = IIRCoeffs32LP/8;
   }

   arm_float_to_q31(IIRCoeffs32LP, IIRCoeffs32LP1, 5*numStages);

   for(int i = 0; i < 5*numStages; i++)
   {
            //printf("%d\r\n", IIRCoeffs32LP1);
   }

   /* 初始化输入输出缓存指针 */
   inputQ31 = &testInput_f32_50Hz_200Hz1;
   outputQ31 = &testOutput1;

   for(i=0; i<TEST_LENGTH_SAMPLES; i++)
   {
            /* 50Hz正弦波+200Hz正弦波，采样率1KHz */
            testInput_f32_50Hz_200Hz1 =arm_sin_q31(2147483647/1000*50*i) / 16+ arm_sin_q31(2147483647/1000*200*i)/16;
   }

   /* 初始化，根据系数的缩放，此函数的最后一个参数要做调整 */
   arm_biquad_cascade_df1_init_q31(&S, numStages, (q31_t *)&IIRCoeffs32LP1, (q31_t *)&IIRStateF321, 1);

   /* 实现IIR滤波，这里每次处理1个点 */
   for(i=0; i < numBlocks; i++)
   {
            arm_biquad_cascade_df1_q31(&S, inputQ31 + (i * blockSize),outputQ31 + (i * blockSize),blockSize);
   }

   /* 放缩系数，testOutput1乘以放缩系数才是准确的 */
   //0.057512362781628162400338055704196449369
   //0.050747576720993219090427572837143088691
   //0.046131802093312926360546555315522709861
   //0.043272238929987591471437724521820200607
   //0.041906571098479158954397405523195629939

   /* 打印滤波后结果 */
   for(i=0; i<TEST_LENGTH_SAMPLES; i++)
   {
            printf("%d, %d\r\n", testInput_f32_50Hz_200Hz1, testOutput1);
   }
}

/*
*********************************************************************************************************
*    函数名: main
*    功能说明: c程序入口
*    形参: 无
*    返回值: 错误代码(无需处理)
*********************************************************************************************************
*/
int main(void)
{
   uint8_t ucKeyCode;             /* 按键代码 */

   bsp_Init();             /* 硬件初始化 */
   PrintfLogo();    /* 打印例程信息到串口1 */

   PrintfHelp();    /* 打印操作提示信息 */

   bsp_StartAutoTimer(0, 100);    /* 启动1个100ms的自动重装的定时器 */


   /* 进入主程序循环体 */
   while (1)
   {
            bsp_Idle();             /* 这个函数在bsp.c文件。用户可以修改这个函数实现CPU休眠和喂狗 */


            if (bsp_CheckTimer(0))    /* 判断定时器超时时间 */
            {
                     /* 每隔100ms 进来一次 */
                     bsp_LedToggle(2);    /* 翻转LED的状态 */
            }

            ucKeyCode = bsp_GetKey();    /* 读取键值, 无键按下时返回 KEY_NONE = 0 */
            if (ucKeyCode != KEY_NONE)
            {
                     switch (ucKeyCode)
                     {
                           case KEY_DOWN_K1:                /* K1键按下 */
                                    arm_iir_q31_lp();
                                    break;

                           default:
                                    /* 其它的键值不处理 */
                                    break;
                     }
            }

   }
}

eric2013 发表于 2022-8-25 10:27:31

更新了下

eric2013 发表于 2022-8-25 16:07:31

再次更新matlab生成系数的溢出问题

KC8022 发表于 2024-2-20 18:15:34

汉哥，使用本帖内的工程移植后，我得到滤波前和滤波后得到的波形似乎不正确，移植的平台是HC32F460(M4内核，且开启了FPU加速)，IAR的DSP库版本是 V1.5.2。
main函数中的功能仅有一次arm_iir_q31_lp调用，所以没有上传main函数
滤波的结果上传在图片里了，能麻烦汉哥帮看看吗，是什么问题？

#include "arm_math.h"
#include "Ripple.h"
typedefunsigned char       BOOL;
typedefsigned char       S8;
typedefunsigned char       U8;
typedefsigned short       S16;
typedefunsigned short    U16;
typedefsigned int          S32;
typedefunsigned int       U32;
typedefunsigned long long U64;
typedefsigned long long    S64;
typedeffloat             F32;
typedefdouble             F64;

#define STAGE5             /* 5阶IIR滤波的个数 */
#define TEST_LENGTH_SAMPLES400 /* 采样点数 */
#define BLOCK_SIZE       400 /* 调用一次arm_biquad_cascade_df1_f32处理的采样点个数 */

U32 blockSize = BLOCK_SIZE;
U32 numBlocks = TEST_LENGTH_SAMPLES/BLOCK_SIZE; /* 需要调用arm_biquad_cascade_df1_f32的次数 */

/* 巴特沃斯低通滤波器系数 80Hz*/
F32 IIRCoeffs32LP =
{
1,2,1, 1.62978785624054456349085739930160343647, -0.85983730736705732411451208463404327631,
1,2,1, 1.438087052475096605519411241402849555016,-0.641077359359069620658999610895989462733,
1,2,1, 1.307285028849323404998017394973430782557,-0.491812237222575276973657310009002685547,
1,2,1, 1.226250602643673470737439856748096644878,-0.39933955836362394764549321735103148967,
1 , 2,1, 1.187550247801674174752406543120741844177,-0.355176532195590921592298627729178406298
};

static q31_t testInput_f32_50Hz_200Hz1; /* 采样点 */
static q31_t testOutput1;             /* 滤波后的输出 */
static q31_t IIRStateF321;                   /* 状态缓存 */

q31_t IIRCoeffs32LP1;

void arm_iir_q31_lp(void)

{
uint32_t i;
arm_biquad_casd_df1_inst_q31 S;
q63_t ScaleValue;
q31_t *inputQ31, *outputQ31;

/* 使用浮点系数先进行一次缩放，然后再使用函数arm_float_to_q31转换为Q31，防止溢出 */
for(int i = 0; i < 5*STAGE; i++)
{
   IIRCoeffs32LP = IIRCoeffs32LP/8;
}
arm_float_to_q31(IIRCoeffs32LP, IIRCoeffs32LP1, 5*STAGE);
for(int i = 0; i < 5*STAGE; i++)
{
   //printf("%d\r\n", IIRCoeffs32LP1);
}
/* 初始化输入输出缓存指针 */
inputQ31 = &testInput_f32_50Hz_200Hz1;
outputQ31 = &testOutput1;
for(i=0; i<TEST_LENGTH_SAMPLES; i++)
{
   /* 50Hz正弦波+200Hz正弦波，采样率1KHz */
   testInput_f32_50Hz_200Hz1 =arm_sin_q31(2147483647/1000*50*i) / 16+ arm_sin_q31(2147483647/1000*200*i)/16;
}

/* 初始化，根据系数的缩放，此函数的最后一个参数要做调整 */
arm_biquad_cascade_df1_init_q31(&S, STAGE, (q31_t *)&IIRCoeffs32LP1, (q31_t *)&IIRStateF321, 1);
/* 实现IIR滤波，这里每次处理1个点 */
for(i=0; i < numBlocks; i++)
{
   arm_biquad_cascade_df1_q31(&S, inputQ31 + (i * blockSize),outputQ31 + (i * blockSize),blockSize);
}

/* 放缩系数，testOutput1乘以放缩系数才是准确的 */
//0.057512362781628162400338055704196449369
//0.050747576720993219090427572837143088691
//0.046131802093312926360546555315522709861
//0.043272238929987591471437724521820200607
//0.041906571098479158954397405523195629939

/* 打印滤波后结果 */
for(i=0; i<TEST_LENGTH_SAMPLES; i++)
{
   printf("%d\r\n", testInput_f32_50Hz_200Hz1);
}
   /* 打印滤波后结果 */
for(i=0; i<TEST_LENGTH_SAMPLES; i++)
{
   printf("%d\r\n", testOutput1);
}
}

eric2013 发表于 2024-2-21 09:14:50

KC8022 发表于 2024-2-20 18:15
汉哥，使用本帖内的工程移植后，我得到滤波前和滤波后得到的波形似乎不正确，移植的平台是HC32F460(M4内核 ...

你这个是浮点API，按说不该有问题，直接跑我原始的测试代码，是否正常。

页: [1]

硬汉嵌入式论坛's Archiver

DSP定点Q31 IIR滤波器的matlab设置，使用函数arm_biquad_cascade_df1_q31