CMSIS-DSP 之 NLMS 自适应滤波器

一、 NLMS 简介

arm CMSIS-DSP 库提供了 LMS 最小均方自适应滤波和改进后的 NLMS 归一化最小均方自适应滤波器，本篇主要介绍使用 NLMS (Normalized Least Mean Squares 归一化最小均方) 算法的 FIR 自适应滤波器 Adaptive Filters。NLMS 滤波器是一类能够“学习”未知传递函数的自适应滤波器，它可以根据输入信号自动调整滤波系数进行数字滤波，无需 Matlab 生成滤波系数。

NLMS 滤波器由两个组件组成：第一个组成部分是标准横向滤波器或 FIR 滤波器；第二个组成部分是系数更新机制。NLMS 滤波器有两个输入信号：“输入”馈送 FIR 滤波器，而“参考输入”对应 FIR 滤波器的期望输出。即，不断更新 FIR 滤波器系数，使得 FIR 滤波器的输出与参考输入匹配。滤波器系数更新机制基于 FIR 滤波器输出和参考输入之间的差值，随着滤波器的适应，该“误差信号”趋向于零。NLMS 滤波功能接收输入和参考输入信号并生成滤波器输出和误差信号。

图 1. NLMS 结构图

二、 NLMS 使用步骤简介

2.1 CMSIS-DSP 中的 NLMS 函数接口

主要用到的是下面这两个函数。

void arm_lms_norm_init_f32 (arm_lms_norm_instance_f32 * S, uint16_t numTaps,

float32_t * pCoeffs, float32_t * pState, float32_t mu, uint32_t blockSize)

LMS 初始化。

参数 1：[in] arm_lms_norm_instance_f32 类型结构体变量。

参数 2：[in] 滤波器系数的个数。

参数 3：[in] 滤波器系数数组的地址。

参数 4：[in] 状态缓冲区地址，用于存储滤波器的历史输入和输出样本。

参数 5：[in] 控制滤波器系数更新的步长，即每次更新权重的大小。

参数 6：[in] 每次调用处理的数据个数。

void arm_lms_norm_f32(arm_lms_norm_instance_f32 * S, const float32_t * pSrc,

float32_t * pRef, float32_t * pOut, float32_t * pErr, uint32_t blockSize)

LMS 滤波。

参数 1：[in] arm_lms_norm_instance_f32 类型结构体变量。

参数 2：[in] 原始输入信号的数据地址。

参数 3：[in] 参考信号的数据地址，需要输入想要的理想输出波形作为参照。

参数 4：[out] 滤波后输出信号的数据地址。

参数 5：[out] 输出信号与参考信号的误差数据地址。

参数 6：[in] 每次调用处理的数据个数，最小为 1，最大为全部样本数。

2.2 实现 NLMS

接着调用这两个函数就可以了，例如下图，我们通过 arm_sin_f32() 创建了一个多频率的正弦波混合信号作为原始输入信号，然后对其进行低通滤波处理。

图 2. 调用 lms 函数

将滤波前后数据绘制出来，效果如下图，可以看到，最初实际输与期望输出的参考信号之间还是有偏差的，通过 lms 算法不断调整，输出信号逐渐靠近参考信号，即误差信号 (绿色) 逐渐收敛为 0。

图 3. 滤波结果

2.3步长 MU

步长控制着每次更新权重的大小，会影响误差信号的收敛速度，步长约小，收敛速度越慢，步长越大，收敛速度越快，若步长过大，则可能导致误差信号无法收敛。下图为不同步长的滤波对比 (为了方便观察误差信号，将其他信号颜色变浅)。

图 4. 不同步长对比

如果放大误差信号观察，对比更明显。

图 5. 放大误差信号

2.4 NLMS 低通、高通、带通、带阻滤波器举例

RT1010 EVK 测试代码如下。

#define TEST_LENGTH_SAMPLES 1024
#define NUMTAPS 32
#define BLOCK_SIZE 32
#define MU 0.2f
#define NUMFRAMES (TEST_LENGTH_SAMPLES / BLOCK_SIZE)

float32_t testInput_f32_50Hz_200Hz_400Hz[TEST_LENGTH_SAMPLES]; //输入信号
static float32_t lmsOutput[TEST_LENGTH_SAMPLES];
float32_t refOutput_lowpass[TEST_LENGTH_SAMPLES];
float32_t refOutput_highpass[TEST_LENGTH_SAMPLES];
float32_t refOutput_bandpass[TEST_LENGTH_SAMPLES];
float32_t refOutput_bandstop[TEST_LENGTH_SAMPLES];
float32_t errOutput[TEST_LENGTH_SAMPLES];
static float32_t lmsStateF32[NUMTAPS + BLOCK_SIZE];
float32_t lmsNormCoeff_f32[32];

int32_t main(void)
{
//------------------------------ MCU ------------------------------
BOARD_ConfigMPU();
BOARD_InitBootPins();
BOARD_InitBootClocks();
BOARD_InitDebugConsole();
PRINTF("\n\r\n\r=====================================================\n\r");
PRINTF("build time : %s %s\n\r", __DATE__, __TIME__);

//------------------------------ CMSIS DSP ------------------------------
PRINTF("\r\n-------------------- LMS test --------------------\r\n");
PRINTF(" arm_signal_converge_example_f32 \r\n");
uint32_t i;
arm_lms_norm_instance_f32 lmsNorm_S;
float32_t *inputF32, *outputF32, *inputREF, *outputERR;
uint32_t blockSize = BLOCK_SIZE;
uint32_t numBlocks = TEST_LENGTH_SAMPLES / BLOCK_SIZE;
arm_status status;
uint32_t index;
float32_t minValue;
for(uint32_t i=0; i<TEST_LENGTH_SAMPLES; i++) //原始输入信号
{
// 50Hz 正弦波 + 200Hz 正弦波 + 400Hz，采样率 1KHz
testInput_f32_50Hz_200Hz_400Hz[i] = arm_sin_f32(2*3.1415926f*50*i/1000) + arm_sin_f32(2*3.1415926f*200*i/1000) + arm_sin_f32(2*3.1415926f*400*i/1000);
// 50Hz 正弦波，采样率 1KHz  //理论 lowpass 滤波后信号
refOutput_lowpass[i] = arm_sin_f32(2*3.1415926f*50*i/1000);

// 200Hz 正弦波 + 400Hz，采样率 1KHz //理论 highpass 滤波后信号
refOutput_highpass[i] = arm_sin_f32(2*3.1415926f*200*i/1000) + arm_sin_f32(2*3.1415926f*400*i/1000);
// 200Hz 正弦波，采样率 1KHz //理论 bandpass 滤波后信号
refOutput_bandpass[i] = arm_sin_f32(2*3.1415926f*200*i/1000);
// 50Hz 正弦波 + 400Hz，采样率 1KHz //理论 bandstop 滤波后信号
refOutput_bandstop[i] = arm_sin_f32(2*3.1415926f*50*i/1000) + arm_sin_f32(2*3.1415926f*400*i/1000);
}
inputF32 = (float32_t *)&testInput_f32_50Hz_200Hz_400Hz[0]; /* 原始波形 */
outputF32 = (float32_t *)&lmsOutput[0]; /* 滤波后输出波形 */
inputREF = (float32_t *)&refOutput_lowpass[0]; /* 参考波形 */
outputERR = (float32_t *)&errOutput[0]; /* 误差数据 */

//------------------------------ low pass ------------------------------
PRINTF("\r\n---------- low pass ----------\r\n");
inputREF = (float32_t *)&refOutput_lowpass[0]; /* 参考波形 */
memset(lmsNormCoeff_f32,0,sizeof(lmsNormCoeff_f32));
memset(lmsStateF32,0,sizeof(lmsStateF32));

/* Initialize the LMSNorm data structure */
arm_lms_norm_init_f32(&lmsNorm_S, NUMTAPS, lmsNormCoeff_f32, lmsStateF32, MU, BLOCK_SIZE);
/* ----------------------------------------------------------------------
* Loop over the frames of data and execute each of the processing
* functions in the system.
* ------------------------------------------------------------------- */
for(i=0; i < NUMFRAMES; i++) //滤波
{
/* Execute the LMS Norm processing function*/
arm_lms_norm_f32(&lmsNorm_S, /* LMSNorm instance */
inputF32 + (i * blockSize), /* Input signal */
inputREF + (i * blockSize), /* Reference Signal */
outputF32 + (i * blockSize), /* Converged Signal */
outputERR + (i * blockSize), /* Error Signal, this will become small as the signal converges */
BLOCK_SIZE); /* BlockSize */
}

for(i = 0; i < TEST_LENGTH_SAMPLES; i++) //输出结果
{
PRINTF("ch: %f, %f, %f, %f\r\n", inputF32[i], inputREF[i], outputF32[i], outputERR[i]); //串口输出
}

//------------------------------ high pass ------------------------------
PRINTF("\r\n---------- high pass ----------\r\n");
inputREF = (float32_t *)&refOutput_highpass[0]; /* 参考波形 */
memset(lmsNormCoeff_f32,0,sizeof(lmsNormCoeff_f32));
memset(lmsStateF32,0,sizeof(lmsStateF32));

arm_lms_norm_init_f32(&lmsNorm_S, NUMTAPS, lmsNormCoeff_f32, lmsStateF32, MU, BLOCK_SIZE);
for(i=0; i < NUMFRAMES; i++)
{
/* Execute the LMS Norm processing function*/
arm_lms_norm_f32(&lmsNorm_S, /* LMSNorm instance */
inputF32 + (i * blockSize), /* Input signal */
inputREF + (i * blockSize), /* Reference Signal */
outputF32 + (i * blockSize), /* Converged Signal */
outputERR + (i * blockSize), /* Error Signal, this will become small as the signal converges */
BLOCK_SIZE); /* BlockSize */
}

for(i = 0; i < TEST_LENGTH_SAMPLES; i++) //输出结果
{
PRINTF("ch: %f, %f, %f, %f\r\n", inputF32[i], inputREF[i], outputF32[i], outputERR[i]); //串口输出
}

//------------------------------ band pass ------------------------------
PRINTF("\r\n---------- band pass ----------\r\n");
inputREF = (float32_t *)&refOutput_bandpass[0]; /* 参考波形 */
memset(lmsNormCoeff_f32,0,sizeof(lmsNormCoeff_f32));
memset(lmsStateF32,0,sizeof(lmsStateF32));
arm_lms_norm_init_f32(&lmsNorm_S, NUMTAPS, lmsNormCoeff_f32, lmsStateF32, MU, BLOCK_SIZE);
for(i=0; i < NUMFRAMES; i++)
{
/* Execute the LMS Norm processing function*/
arm_lms_norm_f32(&lmsNorm_S, /* LMSNorm instance */
inputF32 + (i * blockSize), /* Input signal */
inputREF + (i * blockSize), /* Reference Signal */
outputF32 + (i * blockSize), /* Converged Signal */
outputERR + (i * blockSize), /* Error Signal, this will become small as the signal converges */
BLOCK_SIZE); /* BlockSize */
}
for(i = 0; i < TEST_LENGTH_SAMPLES; i++) //输出结果
{
PRINTF("ch: %f, %f, %f, %f\r\n", inputF32[i], inputREF[i], outputF32[i], outputERR[i]); //串口输出
}

//------------------------------ band stop ------------------------------
PRINTF("\r\n---------- band stop ----------\r\n");
inputREF = (float32_t *)&refOutput_bandstop[0]; /* 参考波形 */
memset(lmsNormCoeff_f32,0,sizeof(lmsNormCoeff_f32));
memset(lmsStateF32,0,sizeof(lmsStateF32));

arm_lms_norm_init_f32(&lmsNorm_S, NUMTAPS, lmsNormCoeff_f32, lmsStateF32, MU, BLOCK_SIZE);
for(i=0; i < NUMFRAMES; i++)
{
/* Execute the LMS Norm processing function*/
arm_lms_norm_f32(&lmsNorm_S, /* LMSNorm instance */
inputF32 + (i * blockSize), /* Input signal */
inputREF + (i * blockSize), /* Reference Signal */
outputF32 + (i * blockSize), /* Converged Signal */
outputERR + (i * blockSize), /* Error Signal, this will become small as the signal converges */
BLOCK_SIZE); /* BlockSize */
}

for(i = 0; i < TEST_LENGTH_SAMPLES; i++) //输出结果
{
PRINTF("ch: %f, %f, %f, %f\r\n", inputF32[i], inputREF[i], outputF32[i], outputERR[i]); //串口输出
}

while (1)
{
}
}

运行效果如下。

图 6. 低通滤波