本文程式碼的實現嚴重依賴前面的一篇文章：

一，小波閾值去噪基本理論

1，小波閾值處理基本理論

小波閾值收縮法是Donoho和Johnstone在1995年提出的，以下便是養活不少學者的三篇基礎論文，引得無數學者在此基礎上優化，或者應用到自己的工程中然後發表相關的論文：

【1】 Donoho D L. De-noising by soft-thresholding. IEEE Trans- actions on Information Theory, 1995, 41(3): 613−627
【2】 Donoho D L, Johnstone I M. Adapting to unknown smooth- ness via wavelet shrinkage. Journal of the American Statistic
Association, 1995, 90(432): 1200−1224
【3】 Donoho D L, Johnstone I M, Kerkyacharian G, Picard D. Wavelet shrinkage: asymptopia? Journal of Royal Statisti-
cal Society Series B, 1995, 57(2): 301−369

所謂閾值去噪簡而言之就是對訊號進行分解，然後對分解後的係數進行閾值處理，最後重構得到去噪訊號。該演算法其主要理論依據是：小波變換具有很強的去資料相關性,它能夠使訊號的能量在小波域集中在一些大的小波係數中;而噪聲的能量卻分佈於整個小波域內.因此,經小波分解後,訊號的小波係數幅值要大於噪聲的係數幅值.可以認為,幅值比較大的小波係數一般以訊號為主，而幅值比較小的係數在很大程度上是噪聲.於是,採用閾值的辦法可以把訊號係數保留,而使大部分噪聲係數減小至零.小波閾值收縮法去噪的具體處理過程為:將含噪訊號在各尺度上進行小波分解,設定一個閾值，幅值低於該閾值的小波係數置為0,高於該閾值的小波係數或者完全保留,或者做相應的“收縮(shrinkage)”處理.最後將處理後獲得的小波係數用逆小波變換進行重構,得到去噪後的訊號.

2，閾值函式的選取

  小波分解閾值去噪中,閾值函式體現了對超過和低於閾值的小波係數不同處理策略，是閾值去噪中關鍵的一步。設w表示小波係數,T為給定閾值,sign(*)為符號函式,常見的閾值函式有：

硬閾值函式：     （小波係數的絕對值低於閾值的置零，高於的保留不變）

軟閾值函式：   （小波係數的絕對值低於閾值的置零，高於的係數shrinkage處理）

值得注意的是：

1) 硬閾值函式在閾值點是不連續的，在下圖中已經用黑線標出。不連續會帶來振鈴，偽吉布斯效應等。

2) 軟閾值函式，原係數和分解得到的小波係數總存在著恆定的偏差,這將影響重構的精度

同時這兩種函式不能表達出分解後係數的能量分佈，半閾值函式是一種簡單而經典的改進方案。見下圖：

                                                                               圖1

3，閾值的確定

選取的閾值最好剛好大於噪聲的最大水平，可以證明的是噪聲的最大限度以非常高的概率低於（此閾值是Donoho提出的），其中根號右邊的這個引數（叫做sigma）就是估計出來的噪聲標準偏差（根據第一級分解出的小波細節係數，即整個det1絕對值係數中間位置的值），本文將用此閾值去處理各尺度上的細節係數，注意所謂全域性閾值就是近似係數不做任何閾值處理外，其他均閾值處理。

最後吐槽一下這個“絕對值係數中間位置的值”

1）如果det1的長度為偶數那麼，這個“中值”便是中間位置的兩個數之和的平均值，比如【2,2,3,5】，中值即是2.5而不是3

2）如果det1的長度為奇數那麼，這個中值就是中間位置的那個數，比如【2,3,5】，中值即3

4，閾值策略

以前寫的ppt挪用過來：

5，一維訊號的多級分解與重構

二，matlab實現一維小波閾值去噪

1，核心庫函式說明

1）wnoisest函式

作用：估計一維小波高頻係數中的噪聲偏差

這個估計值使用的是絕對值中間位置的值（估計的噪聲偏差值）除以0.6745（Median Absolute Deviation / 0.6745），適合0均值的高斯白噪聲

2）wavedec函式

一維訊號的多尺度分解，將返回諸多細節係數和每個係數的長度，在matlab中鍵入“doc wavedec”具體功能一目瞭然

3）waverec函式

一維訊號小波分解係數的重構，將返回重構後的訊號在matlab中鍵入“doc waverec”具體功能一目瞭然，也可以鍵入“open waverec”檢視matlab具體是怎麼做的。

4）wdencmp函式

這個函式用於對一維或二維訊號的壓縮或者去噪，使用方法：
1 [XC,CXC,LXC,PERF0,PERFL2] = wdencmp('gbl',X,'wname',N,THR,SORH,KEEPAPP)
2 [XC,CXC,LXC,PERF0,PERFL2] = wdencmp('lvd',X,'wname',N,THR,SORH)
3 [XC,CXC,LXC,PERF0,PERFL2] = wdencmp('lvd',C,L,'wname',N,THR,SORH)
wname是所用的小波函式，

gbl(global的縮寫)表示每層都採用同一個閾值進行處理，

lvd表示每層用不同的閾值進行處理，

N表示小波分解的層數，

THR為閾值向量，

對於格式(2)(3)每層都要求有一個閾值，因此閾值向量THR的長度為N，

SORH表示選擇軟閾值還是硬閾值(分別取為’s’和’h’)，

引數KEEPAPP取值為1時，則低頻係數不進行閾值量化處理，反之，則低頻係數進行閾值量化。

XC是消噪或壓縮後的訊號，[CXC，LXC]是XC的小波分解結構，

PERF0和PERFL2是恢復和壓縮L^2的範數百分比, 是用百分制表明降噪或壓縮所保留的能量成分。

2，閾值去噪效果

從圖中可以查出，除燥效果還是比較理想的，對於噪聲比較重的地方軟閾值去噪能力更加明顯（因為沒有無噪的訊號參考，這並不能代表他比硬閾值更優秀）。

放大其中的細節部分，便於檢視細節（抱歉，這裡用了C++的處理結果，但是可喜的是他和matlab的結果一模一樣）

3，完整的matlab程式碼

本程式碼就是上述效果的matlab程式！

clc;
clear;
% 獲取噪聲訊號
load leleccum;
indx = 1:3450;
noisez = leleccum(indx);

%訊號的分解
wname = 'db3'; 
lev = 3;
[c,l] = wavedec(noisez,lev,wname);

%求取閾值
sigma = wnoisest(c,l,1);%使用庫函式wnoisest提取第一層的細節係數來估算噪聲的標準偏差
N = numel(noisez);%整個訊號的長度
thr = sigma*sqrt(2*log(N));%最終閾值

%全域性閾值處理
keepapp = 1;%近似係數不作處理
denoisexs = wdencmp('gbl',c,l,wname,lev,thr,'s',keepapp);
denoisexh = wdencmp('gbl',c,l,wname,lev,thr,'h',keepapp);

% 作圖
subplot(311), 
plot(noisez), title('原始噪聲訊號');
subplot(312),
plot(denoisexs), title('matlab軟閾值去噪訊號') ;
subplot(313),
plot(denoisexh), title('matlab硬閾值去噪訊號') ;

三，C加加實現小波閾值去噪

說明：一維訊號的單尺度分解在前一篇文章中已經提及，這裡不再累述，這裡主要再次基礎上的多尺度分解與重構
並且在執行自己編寫的wavedec函式時必須先初始化，初始化的目的是為了獲取訊號的長度，選擇的是什麼小波，以及分解的等級等資訊，然後計算出未來的各種資訊，比如每個等級的係數的size，為了進行一維小波分解的初始化函式如下：

bool  CWavelet::InitDecInfo(
	const int signalLen,//源訊號長度
	const int decScale,//分解尺度
	const int decdbn//db濾波器的編號
	)
{
	if (decdbn != 3)
		SetFilter(decdbn);

	if (signalLen < m_dbFilter.filterLen - 1)
	{
		cerr << "錯誤資訊：濾波器長度大於訊號!" << endl;
		return false;
	}

	int srcLen = signalLen;
	m_msgCL1D.dbn = decdbn;
	m_msgCL1D.Scale = decScale;
	m_msgCL1D.msgLen.resize(decScale + 2);
	m_msgCL1D.msgLen[0] = srcLen;
	for (int i = 1; i <= decScale; i++)
	{
		int exLen = (srcLen + m_dbFilter.filterLen - 1) / 2;//對稱拓延後係數的長度
		srcLen = exLen;
		m_msgCL1D.msgLen[i] = srcLen;
	}
	m_msgCL1D.msgLen[decScale + 1] = srcLen;

	for (int i = 1; i < decScale + 2; i++)
		m_msgCL1D.allSize += m_msgCL1D.msgLen[i];

	m_bInitFlag1D = true;//設定為已經初始化
	return true;
}

1，核心函式的實現

1），訊號的多級分解

注：本函式實現了對訊號的任意級數分解（分解級數不是在此函式指定），分解的全部係數與matlab的結果完全一致

// 一維多尺度小波分解,必須先初始化
//分解的尺度等資訊已經在初始化函式獲取
bool CWavelet::WaveDec(
	double *pSrcData,//要分解的訊號
	double *pDstCeof//分解出來的係數
	)
{
	if (pSrcData == NULL || pDstCeof == NULL)
		return false;

	if (!m_bInitFlag1D)
	{
		cerr << "錯誤資訊：未初始化，無法對訊號進行分解!" << endl;
		return false;
	}

	int signalLen = m_msgCL1D.msgLen[0];
	int decLevel = m_msgCL1D.Scale;

	double *pTmpSrc = new double[signalLen];
	double *pTmpDst = new double[m_msgCL1D.msgLen[1] * 2];

	for (int i = 0; i < signalLen; i++)
		pTmpSrc[i] = pSrcData[i];

	int gap = m_msgCL1D.msgLen[1] * 2;
	for (int i = 1; i <= decLevel; i++)
	{
		int curSignalLen = m_msgCL1D.msgLen[i - 1];
		DWT(pTmpSrc, curSignalLen, pTmpDst);

		for (int j = 0; j < m_msgCL1D.msgLen[i] * 2; j++)
			pDstCeof[m_msgCL1D.allSize - gap + j] = pTmpDst[j];
		for (int k = 0; k < m_msgCL1D.msgLen[i]; k++)
			pTmpSrc[k] = pTmpDst[k];
		gap -= m_msgCL1D.msgLen[i];
		gap += m_msgCL1D.msgLen[i + 1] * 2;
	}

	delete[] pTmpDst;
	pTmpDst = NULL;
	delete[] pTmpSrc;
	pTmpSrc = NULL;

	return true;
}

2），多級分解係數的重構

注：本函式只能還原成原始訊號，還不能重構到某個中間分解結果

// 重構出源訊號
bool CWavelet::WaveRec(
	double *pSrcCoef,//源被分解係數
	double *pDstData//重構出來的訊號，兩者的長度是一樣的
	)
{
	
	if (pSrcCoef == NULL || pDstData == NULL)//錯誤：無記憶體
		return false;

	//從m_msgCL1D中獲取分解資訊
	int signalLen = m_msgCL1D.msgLen[0];//訊號長度
	int decLevel = m_msgCL1D.Scale;//分解級數

	int det1Len = m_msgCL1D.msgLen[1];
	double *pTmpSrcCoef = new double[det1Len * 2];

	for (int i = 0; i < m_msgCL1D.msgLen[decLevel] * 2; i++)
		pTmpSrcCoef[i] = pSrcCoef[i];

	int gap = m_msgCL1D.msgLen[decLevel] * 2;
	for (int i = decLevel; i >= 1; i--)
	{
		int curDstLen = m_msgCL1D.msgLen[i - 1];
		IDWT(pTmpSrcCoef, curDstLen, pDstData);
		if (i != 1)
		{
			for (int j = 0; j < curDstLen; j++)
				pTmpSrcCoef[j] = pDstData[j];
			for (int k = 0; k < curDstLen; k++)
				pTmpSrcCoef[k + curDstLen] = pSrcCoef[k + gap];
			gap += m_msgCL1D.msgLen[i - 1];
		}
	}


	delete[] pTmpSrcCoef;
	pTmpSrcCoef = NULL;
	return true;
}

3），閾值的獲取

注：嚴格依照Donoho的閾值寫的程式碼

// 根據細節係數，以及訊號長度計算閾值
double CWavelet::getThr(
	double *pDetCoef,//細節係數（應該是第一級的細節係數）
	int detLen,//此段細節係數的長度
	bool is2D//當前細節係數是否來自是二維影象訊號的
	)
{
	double thr = 0.0;
	double sigma = 0.0;
	
	for (int i = 0; i < detLen; i++)
		pDetCoef[i] = abs(pDetCoef[i]);

	std::sort(pDetCoef, pDetCoef + detLen);

	if (detLen % 2 == 0 && detLen >= 2)
		sigma = (pDetCoef[detLen / 2-1] + pDetCoef[detLen / 2]) / 2 / 0.6745;
	else
		sigma = pDetCoef[detLen / 2] / 0.6745;

	if (!is2D)//一維訊號
	{
		double N = m_msgCL1D.msgLen[0];
		thr = sigma *sqrt(2.0*log(N));
	}
	else{//二維訊號
		double size = m_msgCL2D.msgHeight[0]*m_msgCL2D.msgWidth[0];
		thr = sigma *sqrt(2.0*log(size));

	}
	return thr;
}

4），高頻係數閾值處理

注：本閾值函式只對高頻係數做處理，不對近似係數處理

// 將係數閾值處理，一維二維均適用
void CWavelet::Wthresh(
	double *pDstCoef,//細節係數（應該是除近似係數外的所有的細節係數）
	double thr,//閾值
	const int allsize,//分解出來的係數的總長度（非）
	const int gap,//跳過最後一層的近似係數
	SORH  ish//閾值函式的選取
	)
{
	//
	if (ish)//硬閾值
	{
		for (int i = gap; i < allsize; i++)
		{
			if (abs(pDstCoef[i]) < thr)//小於閾值的置零，大於的不變
				pDstCoef[i] = 0.0;
		}
	}
	else//軟閾值
	{
		for (int i = gap; i < allsize; i++)
		{
			if (abs(pDstCoef[i]) < thr)//小於閾值的置零，大於的收縮
			{
				pDstCoef[i] = 0.0;
			}
			else
			{
				if (pDstCoef[i] < 0.0)
					pDstCoef[i] = thr - abs(pDstCoef[i]);
				else
					pDstCoef[i] = abs(pDstCoef[i]) - thr;
				
			}
		}
	}

}

5），閾值去噪函式

注：本函式涉及到上面提及的多個函式，此函式是核心的對外介面

bool CWavelet::thrDenoise(
	double *pSrcNoise,//源一維噪聲訊號
	double *pDstData,//去噪後的訊號
	bool isHard//閾值函式的選取，有預設值
	)
{
	if (pSrcNoise == NULL || pDstData == NULL)
		exit(1);

	if (!m_bInitFlag1D)//錯誤：未初始化
		return false;

	double *pDstCoef = new double[m_msgCL1D.allSize];
	WaveDec(pSrcNoise, pDstCoef);//分解出係數

	int Det1Len = m_msgCL1D.msgLen[1];
	int gapDet = m_msgCL1D.allSize - Det1Len;
	double *pDet1 = new double[Det1Len];
	for (int i = gapDet, j = 0; i < m_msgCL1D.allSize; i++, j++)
		pDet1[j] = pDstCoef[i];

	int gapApp = m_msgCL1D.msgLen[m_msgCL1D.Scale];//跳過最後一層的近似係數
	double thr = getThr(pDet1, Det1Len);//獲取閾值
	Wthresh(pDstCoef, thr, m_msgCL1D.allSize, gapApp, isHard);//將細節係數閾值
	WaveRec(pDstCoef, pDstData);//重構訊號

	delete[] pDstCoef;
	pDstCoef = NULL;
	delete[] pDet1;
	pDet1 = NULL;
	return true;
}

2，函式正確性驗證

注：本次測試實現了對原始訊號的10級分解，並與matlab進行了對比，結果完全一致（對比未完全展示）

以下資料為matlab的部分結果，完全相同於上述第二行，其他資料也完全一致。

3，對噪聲訊號閾值去噪

說明：C++處理的噪聲訊號資料先被儲存為txt，該txt檔案是由matlab載入匯出的噪聲訊號，最後C++又將計算結果儲存為txt，載入到matlab中顯示。噪聲去噪結果與matlab一致。

宣告：

本文部分文字學習並整理自網路，部分程式碼參考於網路資源.

如果侵犯了您的版權，請聯絡本人[email protected]，本人將及時編輯掉！

注：本博文為EbowTang原創，後續可能繼續更新本文。本著共享、開源精神可以轉載，但請務必複製本條資訊！

原文地址：http://blog.csdn.net/ebowtang/article/details/40481393

原作者部落格：http://blog.csdn.net/ebowtang

參考資源：

【1】網友，鄒宇華，部落格地址，http://blog.csdn.net/chenyusiyuan/article/details/2862768
【2】《維基百科----小波變換》
【3】喬世傑.小波影象編碼中的對稱邊界延拓法[ J].中國影象圖形學報,2000,5(2):725-729.

【4】MALLAT S.A theory formulti-resolution signal decompo-sition: the wavelet representation[ J]. IEEE Transaction on Pattern Analysis and Machine Intelligence, 1989, 11(4):674-693.

【5】《小波十講》

【6】《小波與傅立葉分析基礎》

【7】岡薩雷斯《數字影象處理》

【8】matlab小波演算法說明文件

【9】閾值去噪鼻祖論文，Donoho, D.L. (1995), "De-noising by soft-thresholding," IEEE Trans. on Inf. Theory, 41, 3, pp. 613–627.