c语言实现小波多级分解与重构（跟matlab效果一致）

前言

最近在做数字信号处理方面的内容，需要用到小波变换来进行降噪滤波，具体实现如下。

一、小波分解代码

#define FILTER_LEN					 (8)			//滤波器长度
#define ADS_WAVELET_LEN             （260）			//数据的长度
//多级小波分解，"db4"
double db4_Lo_D[8] = { -0.0105974017850690, 0.0328830116668852, 0.0308413818355607, -0.1870348117190931, -0.0279837694168599, 0.6308807679398587, 0.7148465705529154, 0.2303778133088964};	
double db4_Hi_D[8] = { -0.2303778133088964, 0.7148465705529154, -0.6308807679398587, -0.0279837694168599, 0.1870348117190931, 0.0308413818355607, -0.0328830116668852, -0.0105974017850690};

/*
*********************************************************************************************************
*   函 数 名: wavedec
*   功能说明: 小波多级分解函数（与matlab效果一致）
*   形    参: srcData	输入数据
			  srcLen	输入数据长度
			  level		分解级数
			  C         存储的细节系数cDi和最后一层的近似系数cA。以5层分解为例，它们的存储结构是，C=[cD1;cD2;cD3;cD4;cD5;cA5]
			  L			源数据长度和各组系数的长度。以5层分解为例，L=[Len;Len1;len2;len3;len4;len5;len5]
*   返 回 值: 无
*********************************************************************************************************
*/
void waveDec(double *srcData, int srcLen, int level, double *C, int *L )
{
	double tempsrc[ADS_WAVELET_LEN] = {0};
	int resultlen = 0;
	int tempsrclen = 0;
	//int filterLen = 8;	//滤波器长度
	int index = 0;
	int curlevel = 0;	//当前级数
	int lenIdx = 0;
	
	double temp = 0;
	
  int i,j,p;

	L[lenIdx++] = srcLen;
    for ( i = 0; i < srcLen; i++ )
        tempsrc[i] = srcData[i];

    tempsrclen = srcLen;
	
	while(1)
	{
		resultlen = ((tempsrclen + FILTER_LEN - 1) >> 1);	//通过分解和下采样获得数据量（单一部分）
		for(i = 0; i < resultlen; i ++ )		//卷积，下采样
		{
			for(j = 0; j < FILTER_LEN; j++)
			{
				//延拓处理
				p = 2 * i - j + 1;
				if((p < 0) && (p > -FILTER_LEN + 1))		
					temp = tempsrc[-p - 1];
				else if ((p > tempsrclen - 1) && (p <= tempsrclen + FILTER_LEN - 2))
					temp = tempsrc[2 * tempsrclen - p - 1];
				else if ((p >= 0) && (p < tempsrclen - 1 + 1))
					temp = tempsrc[p];
				else
					temp = 0;

				C[index] +=  temp * db4_Hi_D[j];			//高频部分（细节）
				C[index + resultlen] +=  temp * db4_Lo_D[j];		//低频部分（近似）			
			}
			index++;
		}

		curlevel++;	//级数自增1
		L[lenIdx++] = resultlen;	//记录分解后数据的长度
		tempsrclen = resultlen;		//
		
		if(curlevel < level)	//判断当前级数是否需要继续进行分解
		{
			memset(tempsrc, 0, ADS_WAVELET_LEN);
			for (i = 0; i < resultlen; i++)
			{
				tempsrc[i] = C[index + i];	//取低频部分的数据继续进行下一次的分解
				C[index + i] = 0;
			}
		}
		else
		{
			L[lenIdx] = resultlen;
			break;
		}		
	}
	return;
}
#endif

二、小波重构代码

#define FILTER_LEN					 (8)			//滤波器长度
#define ADS_WAVELET_LEN             （260）			//数据的长度

double db4_Lo_R[8] = { 0.230377813308855, 0.714846570552542, 0.630880767929590, -0.0279837694169839, -0.187034811718881, 0.0308413818359870, 0.0328830116669829, -0.0105974017849973 };
double db4_Hi_R[8] = { -0.0105974017849973, -0.0328830116669829, 0.0308413818359870, 0.187034811718881, -0.0279837694169839, -0.630880767929590, 0.714846570552542, -0.230377813308855 };
/*
*********************************************************************************************************
*   函 数 名: wrcoef
*   功能说明: 小波多级重构函数（与matlab效果一致）
*   形    参: recData	重构完成输出的数据
			  type		‘d’ / ‘a’ ， 对细节分量(cD)重构 / 对近似分量(cA)重构
			  C         存储的细节系数cDi和最后一层的近似系数cA。以5层分解为例，它们的存储结构是，C=[cD1;cD2;cD3;cD4;cD5;cA5]
			  L			源数据长度和各组系数的长度。以5层分解为例，L=[Len;Len1;len2;len3;len4;len5;len5]
			  level		第 level 层信号
*   返 回 值: 无
*********************************************************************************************************
*/
#define DATALEN1	(ADS_WAVELET_LEN + FILTER_LEN - 1)
#define DATALEN2	(ADS_WAVELET_LEN + FILTER_LEN - 1 + FILTER_LEN - 1)
#define DATALEN3	((ADS_WAVELET_LEN + FILTER_LEN - 1) / 2)
void wrcoef(double *recData, char type, double *C, int *L, int level)
{
	int i,j = 0;
	int curLevel = level;
	int srcLen = 0;
	int recLen = 0;
	int index = 0;
	//int filterLen = 8;	//滤波器长度
	double tempData1[DATALEN1] = {0};
	double tempData2[DATALEN2] = {0};
	double srcData[DATALEN3] = {0};
	
	if(!L[level+1])
		return;
	
	if(type == 'a')
	{
		for(i = 1; i <= level; i++)
		{
			index += L[i];				//计算cA近似量在C中的偏移
		}
		for(i = 0; i < L[level + 1]; i++)
		{
				srcData[i] = C[index + i];	//把cA近似量取出
		}
	}
	else if(type == 'd')
	{
		for(i = 1; i < level; i++)
		{
			index += L[i];				//计算cD近似量在C中的偏移
		}
		for(i = 0; i < L[level]; i++)
		{
				srcData[i] = C[index + i];	//把cD近似量取出
		}
	}
	else
		return;
	
	while(curLevel)
	{
		srcLen = L[curLevel];
		recLen = L[curLevel-1];
		
		//上采样补零
		for(i = 0; i < srcLen * 2; i++)		
		{
			if(i % 2 == 0)
			{
				tempData1[i] = 0;		
			}
			else
			{
				tempData1[i] = srcData[j];
				j++;
			}					
		}
		
		//卷积
		memset(tempData2, 0, DATALEN2 * sizeof(double));		//先清零
		for(i = 0; i < FILTER_LEN; i++)		
		{
			for(j = 0; j < srcLen * 2; j++)
			{

				if(type == 'd' && curLevel == level)
					tempData2[i + j] +=  tempData1[j] * db4_Hi_R[i];
				else
					tempData2[i + j] +=  tempData1[j] * db4_Lo_R[i];
			}
		}
		
		//截取
		for(i = FILTER_LEN - 1; i < recLen + FILTER_LEN - 1; i++)	
		{
			recData[i - (FILTER_LEN - 1)] = tempData2[i];
		}
		
		//为下一级重构做准备
		curLevel--;
		j = 0;
		for(i = 0; i < recLen; i++)
		{ 
			srcData[i] = recData[i];
		}
	}	
	return;
}

二、与 matlab 自有的 wavedec 和 wrcoef 函数作比较

以下结果以30个数据输入，3级分解与重构为例：

matlab 分解的结果：

代码分解的结果：

matlab 重构的结果：


代码重构的结果：

总结

无论从分解还是重构，代码实现的结果都与 matlab 自有函数实现的结果一致，完成。