一維信號小波閾值去噪的C++實現（matlab驗證） · 圖像處理算法

本文代碼的實現嚴重依賴前面的一篇文章： [小波變換Mallat算法的C++實現](http://blog.csdn.net/ebowtang/article/details/40433861) # 基本理論小波閾值收縮法是Donoho和Johnstone提出的,其主要理論依據是,小波變換具有很強的去數據相關性,它能夠使信號的能量在小波域集中在一些大的小波系數中;而噪聲的能量卻分布于整個小波域內.因此,經小波分解后,信號的小波系數幅值要大于噪聲的系數幅值.可以認為,幅值比較大的小波系數一般以信號為主，而幅值比較小的系數在很大程度上是噪聲.于是,采用閾值的辦法可以把信號系數保留,而使大部分噪聲系數減小至零.小波閾值收縮法去噪的具體處理過程為:將含噪信號在各尺度上進行小波分解,設定一個閾值，幅值低于該閾值的小波系數置為0,高于該閾值的小波系數或者完全保留,或者做相應的“收縮(shrinkage)”處理.最后將處理后獲得的小波系數用逆小波變換進行重構,得到去噪后的信號. ### 1，閾值函數的選取 ? 閾值去噪中,閾值函數體現了對超過和低于閾值的小波系數不同處理策略，是閾值去噪中關鍵的一步。設w表示小波系數,T為給定閾值,sign(*)為符號函數,常見的閾值函數有：硬閾值函數： ? ? （小波系數的絕對值低于閾值的置零，高于的保留不變） ? ? ?![](https://box.kancloud.cn/2016-03-02_56d65265edd51.jpg) 軟閾值函數： ?? ? ? ??![](https://box.kancloud.cn/2016-03-02_56d65266099d9.jpg) 值得注意的是： 1) 硬閾值函數在閾值點是不連續的，在下圖中已經用黑線標出。不連續會帶來振鈴，偽吉布斯效應等。 2) 軟閾值函數，原系數和分解得到的小波系數總存在著恒定的偏差,這將影響重構的精度使得重構圖像的邊緣模糊等現象. 同時這兩種函數不能表達出分解后系數的能量分布。見下圖： ![](https://box.kancloud.cn/2016-03-02_56d652661a5a8.jpg) ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?圖1 于是不少文章出現了折衷方案，一種新的閾值函數（非本文重點）： ![](https://box.kancloud.cn/2016-03-02_56d652662f779.jpg) 閾值線參考代碼： ~~~ %Generate?signal?and?set?threshold.??? y?=?linspace(-1,1,100);??? subplot(311);?? plot(y);title('原始線');?? thr?=?0.5;?? %?Perform?hard?thresholding.??? ythard?=?wthresh(y,'h',thr);?? subplot(312);?? plot(ythard);title('硬閾值線');?? %?Perform?soft?thresholding.??? ytsoft?=?wthresh(y,'s',thr);?? subplot(313);?? plot(ytsoft);title('軟閾值線'); ? ~~~ ###2，閾值的確定選取的閾值最好剛好大于噪聲的最大水平，可以證明的是噪聲的最大限度以非常高的概率低于![](https://box.kancloud.cn/2016-03-02_56d652664b3b4.jpg) （此閾值是Donoho提出的），其中根號右邊的這個參數（叫做sigma）就是估計出來的噪聲標準偏差（根據第一級分解出的小波細節系數，即整個det1絕對值系數中間位置的值），本文將用此閾值去處理各尺度上的細節系數，注意所謂全局閾值就是近似系數不做任何閾值處理外，其他均閾值處理。最后吐槽一下這個“絕對值系數中間位置的值” 1）如果det1的長度為偶數那么，這個“中值”便是中間位置的兩個數之和的平均值，比如【2,2,3,5】，中值即是2.5而不是3 2）如果det1的長度為奇數那么，這個中值就是中間位置的那個數，比如【2,3,5】，中值即3 ###3，閾值策略以前寫的ppt挪用過來： ![](https://box.kancloud.cn/2016-03-02_56d652665a5cb.jpg) ### 4，一維信號的分解與重構以下算法如果用簡單的文字描述，可是：先將信號對稱拓延（matlab的默認方式），然后再分別與低通分解濾波器和高通分解濾波器卷積，最后下采樣，最后可以看出最終卷積采樣的長度為floor(n-1)/2+n,如果想繼續分解下去則繼續對低頻系數CA采取同樣的方式進行分解。 ![](https://box.kancloud.cn/2016-03-02_56d6526675e25.jpg) # 一，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，閾值去噪效果從圖中可以查出，除燥效果還是比較理想的，對于噪聲比較重的地方軟閾值去噪能力更加明顯（因為沒有無噪的信號參考，這并不能代表他比硬閾值更優秀）。 ![](https://box.kancloud.cn/2016-03-02_56d6526691543.jpg) 放大其中的細節部分，便于查看細節 ![](https://box.kancloud.cn/2016-03-02_56d65266aa7ea.jpg) ### 3，完整的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進行了對比，結果完全一致（對比未完全展示） ![](https://box.kancloud.cn/2016-03-02_56d65266c0467.jpg) 以下數據就是上述中的第二行 ![](https://box.kancloud.cn/2016-03-02_56d65266d6164.jpg) 附函數測試主函數： ~~~ //測試一維數據的任意級數的分解與還原（重構，只能重構為源信號） int _tmain(int argc, _TCHAR* argv[]) { system("color 0A"); double signal[23] = { 420.2, 423.53, 423.52, 423.35, 424.52, 428, 430.79, 428.92, 420.2, 423.53, 423.52, 423.35, 424.52, 428, 430.79, 428.92, 420.2, 423.53, 423.52, 423.35, 424.52, 428, 430.79 }; int signalLen = sizeof(signal) / sizeof(double); cout << "原始信號：" << endl; for (int i = 0; i < signalLen; i++) cout << signal[i] << " "; CWavelet cw; int decdbn = 3; int declevel = 10; if (!cw.InitDecInfo(signalLen, declevel, decdbn)) { Sleep(5000); exit(1); } cout<<endl<<endl << "您的當前分解等級為" << declevel << "級" << endl; int coefLen = cw.m_msgCL1D.allSize; double *pDst = new double[coefLen]; cw.WaveDec(signal, pDst); cout <<endl<< "分解信號中的各級小波系數：" << endl; int i = cw.m_msgCL1D.msgLen.size()-1; int count = cw.m_msgCL1D.msgLen[i]; for (int j = 0; j < coefLen; j++) { cout << pDst[j] << " "; count--; if (count==0) { i--; count = cw.m_msgCL1D.msgLen[i]; cout << endl; } } cout << endl; double *srcdata = new double[signalLen]; cw.WaveRec(pDst, srcdata); cout << "重構信號：" << endl; for (int i = 0; i < signalLen; i++) cout << srcdata[i] << " " ; delete[] pDst; pDst = NULL; delete[] srcdata; srcdata = NULL; system("pause"); return 0; } ~~~ ### 3，對噪聲信號閾值去噪說明：C++處理的噪聲信號數據先被保存為txt，其是由matlab加載導出的噪聲信號，最后又將計算結果保存為txt，加載到matlab中顯示。噪聲去噪結果與matlab一致。 ![](https://box.kancloud.cn/2016-03-02_56d65266e9dfb.jpg) 附閾值去噪主測試函數： ~~~ //一維閾值去噪 const int signal_size = 3450; int main() { system("color 0A"); ifstream waveIn("noiseSignal.txt"); ofstream waveOut("denoiseSignal.txt"); double *signal = new double[signal_size]; for (int i = 0; i < signal_size; i++) waveIn >> signal[i]; double *pDen = new double[signal_size]; CWavelet cw; int scale = 3; int dbn = 3; cw.InitDecInfo(signal_size,scale,dbn); cw.thrDenoise(signal, pDen, true); for (int i = 0; i < signal_size; i++) waveOut << pDen[i] << endl; delete[] pDen; pDen = NULL; delete[] signal; signal = NULL; system("pause"); return 0; } ~~~ 注：本博文為[EbowTang](http://my.csdn.net/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. **聲明：** **本文部分文字學習并整理自網絡，部分代碼參考于網絡資源.** **如果侵犯了您的版權，請聯系本人tangyibiao520@163.com，本人將及時編輯掉！** [![Flag Counter](https://box.kancloud.cn/2016-03-02_56d6526711f3c.jpg)](http://info.flagcounter.com/4ycf) ![](https://box.kancloud.cn/2016-03-02_56d652671f2c8.jpg)