2019年02月19日 | 基于STM32F4的提升小波（二代小波）分解程序说明

一、主要思路 

原始信号：OrgSig

与基于MALLAT算法的小波变换不同，提升小波变换不产生数组L，只产生C数组。定义如下： 

DWT_C：[cD1 | cD2 | … cDN | cAN]，其中cDx代表第x层的细节系数，cAN代表第N层的近似系数。 

但是，信号长度必须是2的整数次幂。 

由于算法可实现原位计算，因此，每层变换后，系数仍存在原始信号的数组中，格式为：[CD，CA]。下一层再变换时，将CA作为原始信号即可，直到分解结束。 

每层变换的步骤：分裂->提升（多层预测/更新）->合并 

各层提升的系数由MATLAB中的liftwave函数计算得到。基于MATLAB中此系数的性质，在此算法中，无论预测还是更新，都采用同一个函数；无论变换还是逆变换都采用加法运算。但是，逆变换中滤波器的符号取反。 

与Mallat算法相比，不需要过多临时变量，不需要过多大的临时数组。需要的程序堆栈比较少。 

二、函数原型 

1、 提升格式更新函数DWT_lsupdate

/****************************************

**将数据进行提升格式更新/预测

//V1.00   实现基本功能 2016-10-12 14:40:59

* @原理：

1、与MATLAB中lsupdate函数的功能类似

* @return 正常则返回1，错误则返回0

*****************************************/

uint16_t DWT_lsupdate(

    float32_t* p_Sig1,               //信号1，待更新的信号

    float32_t* p_Sig2,               //信号2，另一路信号

    uint16_t SigLen,                //信号1和信号2的长度

    const float32_t* p_FilterCoef,  //滤波器系数

    const uint8_t FilterCoefLen,    //滤波器系数长度

    const int8_t DF,                //滤波时的延迟

    int8_t IsLWT                    //LWT还是ILWT？1代表LWT，-1代表ILWT

)

2、 分裂函数DWT_split

/****************************************

**分裂函数

//V1.00   实现基本功能 2016-10-12 15:04:09

* @功能：分裂后，data的前半部分为奇数序号的点（1,3,5...），后半部分为偶数序号的点（0,2,4...）

* @原理：

* @return 正常返回1，否则返回0

*****************************************/

uint8_t DWT_split(

    float32_t* data,  //原始数据

    uint16_t n        //数据长度

)

3、 合并函数

/****************************************

**合并函数

//V1.00   实现基本功能 2016-10-13 15:52:30

* @功能：分裂函数的反函数，data的前半部分为奇数序号的点（1,3,5...），后半部分为偶数序号的点（0,2,4...），

         此函数将数据按正常顺序排列

* @原理：

* @return 正常返回1，否则返回0

*****************************************/

uint8_t DWT_unite(

    float32_t* data,  //原始数据

    uint16_t n        //数据长度

)

4、 提升小波变换函数DWT_lwt

/****************************************

**提升小波变换，即1层提升小波分解

//V1.00   实现基本功能 2016-10-12 15:04:09

* @原理：

   分裂->更新/预测->合并

* @return 正常则返回变换后近似系数和细节系数的长度，错误则返回0

*****************************************/

uint16_t DWT_lwt(

    float32_t* p_OrgSig,     //原始信号

    uint16_t OrgSigLen      //信号长度

)

5、 提升小波反变换函数DWT_ilwt

/****************************************

**提升逆小波变换，即1层提升小波重构

//V1.00   实现基本功能 2016-10-13 14:52:01

* @原理：

       合并->更新/预测->分裂

* @return 正常则返回重构后的信号长度，错误则返回0

*****************************************/

uint16_t DWT_ilwt(

    float32_t* p_LWT_C,      //LWT_C数组

    uint16_t LWT_C_Len       //数组长度

)

6、 提升小波反变换函数DWT_ilwt

/****************************************

**提升逆小波变换，即1层提升小波重构

//V1.00   实现基本功能 2016-10-13 14:52:01

* @原理：

       合并->更新/预测->分裂

* @return 正常则返回重构后的信号长度，错误则返回0

*****************************************/

uint16_t DWT_ilwt(

    float32_t* p_LWT_C,      //LWT_C数组

    uint16_t LWT_C_Len       //数组长度

)

7、 提升小波分解函数DWT_lwtWaveDec

/****************************************

**提升小波分解，可实现N层提升小波分解

//V1.00   实现基本功能 2016-10-13 10:36:56

* @原理：

      1、原位分解，分解后原始数据即变为C数组

* @return 正常则返回1，错误则返回0

*****************************************/

uint16_t DWT_lwtWaveDec(

    float32_t* p_OrgSig,    //原始信号

    uint16_t OrgSigLen,     //信号长度

    uint16_t DecLevel       //分解层数

)

8、 提升小波重构函数DWT_lwtWaveRec

/****************************************

**提升小波重构，可实现N层提升小波重构

//V1.00   实现基本功能 2016-10-13 16:05:13

* @原理：

   1、原位分解，分解后C数组即变为原始数据

* @return 正常则返回1，错误则返回0

*****************************************/

uint16_t DWT_lwtWaveRec(

    float32_t* p_LWT_C,     //LWT_C数组

    uint16_t LWT_C_Len,     //数组长度

    uint16_t DecLevel       //分解层数

三、移植过程 

1、 根据算法研究结果，确定需要进行小波分解的信号长度、小波函数和分解层数 

2、 修改.h文件 

因为是原位运算，不需要根据信号长度和小波分解层数去预定义数组L和C，因此不需要define信号长度和小波分解层数，直接在信号缓存和小波分解时指定即可。 

a、根据所选定小波的提升格式，修改提升次数和每次提升时滤波器的最大长度

#define LWT_LS_LEVEL          5 //提升的次数，即'd'和'p'的个数

#define LWT_Filter_Len_MAX     3 //各层算子的点数，取最大者，点数没这么多的补零

3、 修改.c文件 

a、修改提升格式中滤波器的系数

//以db4为例，在MATLAB中计算出db4的提升格式如下：

//% liftwave('db4')

//%

//% {   'd',-0.322275887997141,1;

//%     'p',[-1.11712360516059,-0.300142258748544],0;

//%     'd',[-0.0188083527262439,0.117648086798478],2;

//%     'p',[2.13181671275522,0.636428271190659],0;

//%     'd',[-0.469083478911028,0.140039237732612,-0.0247912381571950],0;

//%          0.734124527683251,1.36216672007377,[]}

//db4

const float32_t LWT_Filter_Coef[LWT_LS_LEVEL][LWT_Filter_Len_MAX] =

 {{ -0.322275887997141},                                          //1

  { -1.11712360516059, -0.300142258748544},                       //2

  { -0.0188083527262439, 0.117648086798478},                      //3

  {2.13181671275522, 0.636428271190659},                          //4

  { -0.469083478911028, 0.140039237732612, -0.0247912381571950}   //5

};

b、修改提升格式中各层滤波器系数的长度、滤波延迟，以及最终的归一化系数

//db4

const uint8_t LWT_Filter_Len[LWT_LS_LEVEL] = {1, 2, 2, 2, 3};      //滤波器长度

const int8_t LWT_DF[LWT_LS_LEVEL] = {1, 0, 2, 0, 0};               //滤波器延迟

const float32_t LWT_NormCoef[2] = {0.734124527683251, 1.36216672007377}; //归一化系数

const float32_t DWT_Lo_D[DWT_FILTER_LEN] = {  0.0352,   -0.0854,   -0.1350,    0.4599,    0.8069,    0.3327};

const float32_t DWT_Hi_D[DWT_FILTER_LEN] = { -0.3327,    0.8069,   -0.4599,   -0.1350,    0.0854,    0.0352};

const float32_t DWT_Lo_R[DWT_FILTER_LEN] = {  0.3327,    0.8069,    0.4599,   -0.1350,   -0.0854,    0.0352};

const float32_t DWT_Hi_R[DWT_FILTER_LEN] = {  0.0352,    0.0854,   -0.1350,   -0.4599,    0.8069,   -0.3327};

)

4、 使用分解和重构函数 

在程序中合适的位置，缓存或预定义一段原始数据OrgSig，分解后的系数即保存在OrgSig中，再将OrgSig作为重构系数，重构之后的信号仍然保存在OrgSig中。使用例程如下：

 float32_t OrgSig[DWT_SIG_LEN] = {1, 2, 1, -1, 1, 2, -1, -2, 1, 2, 3, 4, 5, 6, 7, 8, 1, 12, 13, 1, 3, 6, 8, 1, 4, 4, 4, 2, 8, 10};

DWT_lwtWaveDec(OrgSig2,32,3);   //32是数据长度，3是分解层数

DWT_lwtWaveRec(OrgSig2,32,3);

四、注意事项 

1、 堆栈设置问题 

小波变换需要的临时变量较大，当信号长度较大时，可能会引起HardFault，进而进入函数HardFault_Handler死循环。这时，需要在启动文件startup_stm32f40_41xxx.s中修改堆区大小。 

如：

; Stack_Size      EQU     0x00000400 //默认设置是这个

Stack_Size      EQU     0x0000FF00

五、测试结果 

1、 对于采样率为360Hz的ECG信号，利用db3对500个点进行4层小波分解后再重构，得到结果对比如下图，二者相关系数为1.0000。