DSP

基于TI DSP C6678的自适应LMS算法的实时处理

2019-07-13 19:58发布生成海报

站内文章 / DSP

11905 0

原文转载于：http://blog.csdn.net/hj199404182515/article/details/52613969 非常感谢。关于自适应LMS的理论基础已经非常的成熟，随便找一本关于自适应滤波器的书就会有介绍相关的内容，有的还可出了它的具体算法，但是还没有一本书有讲过怎样编写能够时实（Real Time）处理的基于C的自适应LMS算法（至少我没有见过），而这个应该说是将算法用于时实处理最为关键的一步，本次我就给大家分享一下关于我是如何用C语言来实现时实自适应LMS算法的。
说到时实处理，如果有自己写过时实FIR或者IIR滤波器程序的，可能对时实处理比较的熟悉，那么什么是时实处理呢？怎样来实现时实处理呢？

什么是实时处理？

简单来说，就是一个处理单元（一段程序代码）在某一时刻得到了输入的数据，然后立即处理得到的数据，并把要输出的数据送入输出单元进行输出，这样就完成了一次时实处理，然后等待下一时刻输入数据，然后重复上述过程。可以看到上述过程是“即刻采集即刻处理”，这就叫时实处理。而非时实处理就是先采集整个时间段的数据，采集完成后再进行处理，是“先采集后处理”所以是非时实的。

怎样实现时实处理？

上一段解释了什么是实时处理，下面来介绍怎样用处理器（比如单片机、DSP、ARM啥的）进行实时处理，通过上一段我们知道需要一个能够保证处理单元在特定的时刻进行处理的信号，这个信号可以有很多种形式，一般常用的就是利用处理器里面的定时器计数一段规定好的时间后然后产生一个溢出中断，并把处理单元的代码写到这个定时器中断程序中。这样定时器每隔一段时间就会就会进入一次中断并执行处理单元的程序，这样就实现了实时处理。

基于C的LMS算法

关于LMS相关的理论和Matlab仿真可以查看博客文章（浅谈自适应滤波器）这里不再赘述，下面将给出LMS算法的C语言代码
LMS.c文件

/*******************************************************************************************************************************************************************
 * 文件名：LMS.c
 * 作者：    IRONMAN
 * 日期：    2016.09.04
 * 版本：     V1.0
 * 说明：  此源文件基于c66x内核实现了最小均方算法(LMS)
 *
 * 算法递推形式：
 * 估计瞬时误差：        e(k) = d(k) - x'(k)w(k)                  (1)
 * 估计滤波系数矢量：w(k+1) = w(k) + 2 niu e(k)x(k)   (2)
 * 初始条件：                w(0) = 0;
 *
 * *****************************************************************************************************************************************************************/
#include 
#include 
#include 
#include "LMS.h"

#pragma DATA_ALIGN(lms_x, 8)
float lms_x[LMS_M];                                                         //输入信号向量存储数组
#pragma DATA_ALIGN(lms_w_forward, 8)
float lms_w_forward[LMS_M];                                       //用于存储未来一时刻滤波系数向量w(k+1)
#pragma DATA_ALIGN(lms_w, 8)
float lms_w[LMS_M];                                                       //用于存储当前时刻滤波系数矢量w(k)

Adaptive_Filter_In    lms_param_in;
Adaptive_Filter_Out  lms_param_out;

void LMS_Gradient_Instantaneous_Estimates(Adaptive_Filter_In *lms_in, Adaptive_Filter_Out* lms_out)
{
    int i;
    static int FIR_order;
    static unsigned char First_in_flag = 1;
    static float *w_ptr, *w_forward_ptr, *Temp_w_ptr;
    static float *x_ptr;
    float  temp;

    if(First_in_flag )
    {
        First_in_flag = 0;
        FIR_order = lms_in->length_x;

        memset((void *)lms_w_forward, 0, FIR_order*4);
        memset((void *)lms_w, 0, FIR_order*4);

        w_forward_ptr = lms_w_forward;
        w_ptr = lms_w;

        x_ptr = lms_in->x_ptr;
    }

/******************************************************************************************************************************/
/**************************************(1)式计算部分****************************************************************************/
/******************************************************************************************************************************/
    lms_out->y =  DSPF_sp_dotprod(x_ptr,w_ptr,FIR_order);               //一共需要82个时钟周期
    lms_out->error = lms_in->d - lms_out->y;
/******************************************************************************************************************************/
/**************************************(2)式计算部分****************************************************************************/
/******************************************************************************************************************************/
    temp = 2*LMS_NIU*lms_out->error;                                 //第二部分一共需要380个时钟周期

    for(i=0; i
    {
        w_forward_ptr[i]      = w_ptr[i]     + temp*x_ptr[i];
        w_forward_ptr[i+1] = w_ptr[i+1] + temp*x_ptr[i+1];
        w_forward_ptr[i+2] = w_ptr[i+2] + temp*x_ptr[i+2];
        w_forward_ptr[i+3] = w_ptr[i+3] + temp*x_ptr[i+3];
        w_forward_ptr[i+4] = w_ptr[i+4] + temp*x_ptr[i+4];
        w_forward_ptr[i+5] = w_ptr[i+5] + temp*x_ptr[i+5];
        w_forward_ptr[i+6] = w_ptr[i+6] + temp*x_ptr[i+6];
        w_forward_ptr[i+7] = w_ptr[i+7] + temp*x_ptr[i+7];
    }
/**************************************************************************************************************************************/
/************************************************************变量更替*******************************************************************/
/**************************************************************************************************************************************/
    Temp_w_ptr = w_forward_ptr;        //新旧滤波系数矢量指针交换
    w_forward_ptr = w_ptr;
    w_ptr = Temp_w_ptr;
}

LMS.h文件

#ifndef LMS_H_
#define LMS_H_

#include "Adaptive_filter.h"

#define LMS_M                                          16                //定义FIR滤波器的阶数   为保证高速运算须为8的倍数
#define LMS_NIU                                      0.0001              //定义收敛因子


extern float lms_x[LMS_M];
extern Adaptive_Filter_In    lms_param_in;
extern Adaptive_Filter_Out  lms_param_out;

extern void LMS_Gradient_Instantaneous_Estimates(Adaptive_Filter_In *lms_in, Adaptive_Filter_Out* lms_out);

#endif

Adaptive_filter.h文件

#ifndef ADAPTIVE_FITLER_H_
#define ADAPTIVE_FITLER_H_

//定义输入参数结构体
typedef struct
{
    float *x_ptr;         //指向输入数组矢量X的指针
    int length_x;        //矢量X的长度，也即FIR滤波器的长度
    float d;                //输入参考数据
}Adaptive_Filter_In;
//定义输出参数结构体
typedef struct
{
    float y;                //滤波器的输出
    float error;        //误差输出
}Adaptive_Filter_Out;


#endif /* ADAPTIVE_FITLER_H_ */

一共三个文件，一个C文件，两个h文件。LMS.c文件是算法的主要部分，LMS.h文件主要起到参数传递作用方便其它的模块调用，Adaptive_filter.h主要定义了自适应滤波器的输入和输出数据的结构体。
运行的处理器是TMS320C6678 DSP，下面将给出在上面运行的结果。

基于TMS320C6678运行结果

首先给出在主函数中调用和产生仿真信号的代码
实验的类型是信号增强，信号的类型和Matlab上的一样以便做对比。程序中调用了dsplib.h中的向量点乘运算函数，如果不是在DSP上或者没有DSPlib的可以自己编写该函数，也不难。

#include 
#include 
#include 
#include 
#include 
#include "KeyStone_common.h"
#include "Keystone_DDR_Init.h"
#include "LMS.h"

#define PI                            3.141592654
#define Fs                            10000
#define Ts                            1/Fs
#define F1                           2
#define F2                           10
#define F3                           20
#define F4                           500
#define NUM_OF_POINT    (32*1024)

#pragma DATA_SECTION(signal, ".DDR_buffer");
float signal[NUM_OF_POINT];
#pragma DATA_SECTION(noise, ".DDR_buffer");
float noise[NUM_OF_POINT];
#pragma DATA_SECTION(signal_noise, ".DDR_buffer");
float signal_noise[NUM_OF_POINT];
#pragma DATA_SECTION(out_error, ".DDR_buffer");
float out_error[NUM_OF_POINT];
#pragma DATA_SECTION(out_y, ".DDR_buffer");
float out_y[NUM_OF_POINT];

//init the system
void System_init(void)
{
    KeyStone_common_CPU_init();/*enable TSC, memory protection interrupts, EDC for internal RAM;clear cache; protect L1 as cache*/
    KeyStone_common_device_init();/*print device information. Enable memory protection interrupts, EDC for MSMC RAM*/
    KeyStone_Exception_cfg(TRUE);//enable exception handling
    CACHE_setL1PSize(CACHE_L1_32KCACHE);
    CACHE_setL1DSize(CACHE_L1_32KCACHE);
    CACHE_setL2Size(CACHE_0KCACHE);
    KeyStone_main_PLL_init(100, 12, 1);   //DSP core speed: 100*12/1=1200MHz
    KeyStone_PASS_PLL_init(100, 21, 2);
    KeyStone_DDR_init(66.66667, 20, 1, NULL);       //DDR init
}

void make_data(void)
{
    int i;
    for(i=0;isignal[i] = sinsp(2*PI*F1*i*Ts) + 0.5*sinsp(2*PI*F2*i*Ts) + 0.25*sinsp(2*PI*F3*i*Ts);
        noise[i] = 5*sinsp(2*PI*F4*i*Ts+PI/2);
        signal_noise[i] = signal[i] + 0.2*noise[i];
        if(i>0)
            signal_noise[i] += 0.15*noise[i-1];
        if(i>1)
            signal_noise[i] += 0.1*noise[i-2];
    }
}

void  main(void)
{
    int i, j;
    System_init();
    make_data();
    memset((void *)lms_x, 0, sizeof(lms_x)*sizeof(float));
    for(i=0;ifor(j=0;j<=i;j++)
            {
                lms_x[j] = noise[i-j];
            }
        }
        else
        {
            for(j=0;j0];
        lms_param_in.length_x = LMS_M;

        LMS_Gradient_Instantaneous_Estimates(&lms_param_in, &lms_param_out);                    //运行瞬时梯度估计LMS算法 耗时514个时钟周期

        out_error[i] = lms_param_out.error;
        out_y[i] = lms_param_out.y;
    }
    while(1);
}