小波學(xué)習(xí)之一(單尺度一維離散小波變換DWT的Mallat算法C++和MATLAB實(shí)現(xiàn))來(lái)源：/v_hyx/article/details/85570711Mallat算法離散序列的Mallat算法分解公式如下:&遵璧)=&3)=云頁(yè)柯伽-時(shí)劣_]S)=GS出A}(n)=g6U)4.(2z?-如其中，H(n)、G(n)分別表示所選取的小波函數(shù)對(duì)應(yīng)的低通和高通濾波器的抽頭系數(shù)序列。從Mallat算法的分解原理可知，分解后的序列就是原序列與濾波器序列的卷積再進(jìn)行隔點(diǎn)抽取而來(lái)。離散序列的Mallat算法重構(gòu)公式如下：W-2以j_/..上）+S￡（.n-戲）0/_J必其中，h(n)、g(n)分別表示所選取的小波函數(shù)對(duì)應(yīng)的低通和高通濾波器的抽頭系數(shù)序列。2小波變換實(shí)現(xiàn)過(guò)程(C/C++)2.1小波變換結(jié)果序列長(zhǎng)度小波的Mallat算法分解后的序列長(zhǎng)度由原序列長(zhǎng)SoureLen和濾波器長(zhǎng)FilterLen決定。從Mallat算法的分解原理可知，分解后的序列就是原序列與濾波器序列的卷積再進(jìn)行隔點(diǎn)抽取而來(lái)。即分解抽取的結(jié)果長(zhǎng)度為(SoureLen+FilterLen-1)/2。2.2獲取濾波器組對(duì)于一些通用的小波函數(shù)，簡(jiǎn)單起見(jiàn)，可以通過(guò)Matlab的wfilters（‘wavenam獲取4個(gè)濾波器；特殊的小波函數(shù)需要自行構(gòu)造獲得。下面以db1小波函數(shù)（Haar小波）為例，其變換與重構(gòu)濾波器組的結(jié)果如下：[cpp]viewplaincopyprint?//matlab輸入獲取命令>>[Lo_D,Hi_D,Lo_R,Hi_R]=wfilters（'dbl'）//獲取的結(jié)果Lo_D=0.70710.7071Hi_D=-0.70710.7071Lo_R=0.70710.7071Hi_R=0.7071-0.70712.3信號(hào)邊界延拓在Mallat算法中，假定輸入序列是無(wú)限長(zhǎng)的，而實(shí)際應(yīng)用中輸入的信號(hào)是有限的采樣序列，這就會(huì)出現(xiàn)信號(hào)邊界處理問(wèn)題。對(duì)于邊界信號(hào)的延拓一般有3種方法，即零延拓、對(duì)稱延拓和周期延拓。3種延拓方法比較情況如下：對(duì)于正交小波變換來(lái)說(shuō)，前兩種延拓方法實(shí)現(xiàn)起來(lái)比較簡(jiǎn)單，但重建時(shí)會(huì)產(chǎn)生邊界效應(yīng)，而且分解的層數(shù)越多，產(chǎn)生的邊界效應(yīng)越顯著。零延拓方法給人一種跳躍的感覺(jué)。至于對(duì)稱性延拓，由于正交小波濾波器一般都是非對(duì)稱性的（Haar小波基雖然是正交的，但它是非連續(xù)的），重建圖象給人一種錯(cuò)位的感覺(jué)。相比較而言，只有最后一種延拓方式可以得到比較精確的重建結(jié)果，它不僅能保證分解與重建正確計(jì)算，而且恢復(fù)的質(zhì)量也好。不過(guò)，周期性延拓方法雖然是常用的三種方法中比較好的方法，但會(huì)導(dǎo)致信號(hào)邊緣的非連續(xù)性，從而會(huì)使得較高頻率（子帶）層的小波系數(shù)很大，即使信號(hào)本身相當(dāng)平滑。從信號(hào)壓縮的角度看，大的系數(shù)是希望避免的。信號(hào)的對(duì)稱延拓可避免邊緣的非連續(xù)性問(wèn)題。然而，對(duì)稱延拓只能和對(duì)稱的小波濾波器一起適用。如果降低正交性要求，選擇雙正交小波變換，對(duì)稱性延拓不失為一種好的方法。周期延拓可適用于任何小波變換，但可能導(dǎo)致輸入序列邊緣的不連續(xù)，使得高頻系數(shù)較大。而對(duì)稱延拓則避免了輸入序列邊界的不連續(xù)，是當(dāng)前廣泛采用的一種延拓方法。下式中給出了最常用的對(duì)稱延拓表達(dá)式。n—1)T-1)￡n<0=]/(/I)05EEN-1[f(ZsLEN-n-l)sLEN-l<nSsLEN+filterLEN-2當(dāng)原序列長(zhǎng)sLEN為偶數(shù)時(shí)延拓后的序列長(zhǎng)為sLEN+2*(filterLEN)，而原序列長(zhǎng)為奇數(shù)時(shí)則需要在右端再延拓一個(gè)元素。注：在Matlab中默認(rèn)使用了對(duì)稱延拓。2.4小波分解在db1小波函數(shù)下，離散序列的Mallat算法分解公式展開(kāi)如下：=XM左)4(2丑-k)=g(0)4(2z?)+月⑴A^2n-1)嘰3)=￡G(k)AjQn—k)=G9)Aj5+G⑴出(2丑-1)其它的db小波，不再詳述。小波分解C++源碼如下。[cpp]viewplaincopyprint?/***@brief小波變換之分解*@paramsourceData源數(shù)據(jù)*@paramdataLen源數(shù)據(jù)長(zhǎng)度*@paramdb過(guò)濾器類型*@paramcA分解后的近似部分序列-低頻部分*@paramcD分解后的細(xì)節(jié)部分序列-高頻部分*@return*正常則返回分解后序列的數(shù)據(jù)長(zhǎng)度，錯(cuò)誤則返回-1*/intWavelet::Dwt(double*sourceData,intdataLen,Filterdb,double*cA,double*cD){if(dataLen<2)return-1;if((NULL==sourceData)||(NULL==cA)||(NULL==cD))return-1;17.m_db=db;intfilterLen=m_db.length;intn,k,p;intdecLen=(dataLen+filterLen-1)/2;doubletmp=0;cout<<"decLen="<<decLen<<endl;24.for(n=0;n<decLen;n++){cA[n]=0;cD[n]=0;for(k=0;k<filterLen;k++){p=2*n-k+1;32.//信號(hào)邊沿對(duì)稱延拓if((p<0)&&(p>=-filterLen+1))tmp=sourceData[-p-1];elseif((p>dataLen-1)&&(p<=dataLen+filterLen-2))tmp=sourceData[2*dataLen-p-1];elseif((p>=0)&&(p<dataLen-1))tmp=sourceData[p];elsetmp=0;42.//分解后的近似部分序列-低頻部分cA[n]+=m_db.lowFilterDec[k]*tmp;45.//分解后的細(xì)節(jié)部分序列-高頻部分cD[n]+=m_db.highFilterDec[k]*tmp;TOC\o"1-5"\h\z}49.}51.returndecLen;}2.5小波重構(gòu)[cpp]viewplaincopyprint?/***@brief小波變換之重構(gòu)*@paramcA分解后的近似部分序列-低頻部分*@paramcD分解后的細(xì)節(jié)部分序列-高頻部分*@paramcALength輸入數(shù)據(jù)長(zhǎng)度*@paramdb過(guò)濾器類型*@paramrecData重構(gòu)后輸出的數(shù)據(jù)*/voidWavelet::Idwt(double*cA,double*cD,intcALength,Filterdb,double*recData){if((NULL==cA)||(NULL==cD)||(NULL==recData))return;13.m_db=db;intfilterLen=m_db.length;16.17.intn,k,p;9.intrecLen=2*cALength-filterLen+1;cout<<"recLen="<<recLen<<endl;20.21.for(n=0;n<recLen;n++)22.23.recData[n]=0;24.24.for(k=0;k<cALength;k++)25.26.p=n-2*k+filterLen-1;27.28.//信號(hào)重構(gòu)29.29.if((p>=0)&&(p<filterLen))1.32.recData[n]+=m_db.lowFilterRec[p]*cA[k]+m_db.highFilterRec[p]*cD[k];//cout<<"recData["<<n<<"]="<<recData[n]<<endl;6.37.}2.6C++實(shí)現(xiàn)結(jié)果3小波變換實(shí)現(xiàn)(MATLAB)使用matlab小波工具箱實(shí)現(xiàn)db4的情況如下。1、MatlabDB4.m文件內(nèi)容。[cpp]viewplaincopyprint?%加載txt數(shù)據(jù)示例s=importdata('data2.txt');%loaddata2.txt;subplot(521);plot(s);%畫(huà)出原始信號(hào)的波形圖title('原始信號(hào)’)；5.[cA,cD]=dwt(s,'db4');%采用db4小波并對(duì)信號(hào)進(jìn)行一維離散小波分解。y=idwt(cA,cD,'db4');%一維離散小波反變換subplot(522);plot(cA);%畫(huà)出波形圖title('MATLAB低頻部分dwt-cA');11.subplot(523);plot(cD);%畫(huà)出波形圖title('MATLAB高頻部分dwt-cD');15.subplot(524);plot(y);%畫(huà)出波形圖title('MATLAB重構(gòu)idwt');2、波形圖如下。4小結(jié)在此，采用C++實(shí)現(xiàn)了dbN系列