[置頂] 數字圖像處理,一維信號小波閾值去噪的C++實現


本文代碼的實現嚴重依賴前面的一篇文章:

小波變換Mallat算法的C++實現


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

      本博文根據小波的分解與重構原理,實現了基於硬閾值和軟閾值函數的一維小波閾值去噪的C++版本,最終結果與matlab庫函數運算結果完全一致。

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的默認方式),然后再分別與低通分解濾波器和高通分解濾波器卷積,最后下采樣,最后可以看出最終卷積采樣的長度為floor(n-1)/2+n,如果想繼續分解下去則繼續對低頻系數CA采取同樣的方式進行分解。



二,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一致。




鑒於不少人要代碼,出於不誤導人,給一些我自己的建議: 1,對程序要有自己的思考,要有自己的驗證和學習過程,寫時當時我還是個C++新手!點擊下載完整程序 2,二維小波變換不建議用,因為速度太慢了(自己徒手寫的),一維可以用。 3,推薦你另外的小波庫網站,http://wavelet2d.sourceforge.net/


聲明:

本文部分文字學習並整理自網絡,部分代碼參考於網絡資源.

如果侵犯了您的版權,請聯系本人tangyb7172@foxmail.com,本人將及時編輯掉!


注:本博文為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.


Flag Counter


注意!

本站转载的文章为个人学习借鉴使用,本站对版权不负任何法律责任。如果侵犯了您的隐私权益,请联系我们删除。



 
粤ICP备14056181号  © 2014-2021 ITdaan.com