基于STM32F4的提升小波（二代小波）分解程序说明

注：本文是程序的说明和实现思路，代码见：https://download.youkuaiyun.com/download/hnxyxiaomeng/10301718

一、主要思路
原始信号：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。

六、说明
程序为原创，文章仅是程序说明。