sym8matlab,【Matlab学习手记】sym8小波滤波

提供sym8小波&＃xff0c;四层全局软阈值滤波源代码&＃xff0c;采用Matlab语言编写&＃xff0c;可移植性强。

源代码

clear;clc;

load leleccum;

indx &＃61; 1:3450;

noisez &＃61; leleccum(indx);

wname &＃61; &＃39;sym8&＃39;;

lev &＃61; 4;

[c,l] &＃61; wavedec(noisez,lev,wname);

[Lo_D,Hi_D,Lo_R,Hi_R] &＃61; wfilters(wname);

% threshold value

sigma &＃61; wnoisest(c,l,1);%使用库函数wnoisest提取第一层的细节系数来估算噪声的标准偏差

N &＃61; numel(noisez);%整个信号的长度

thr &＃61; sigma*sqrt(2*log(N));%最终阈值

%全局阈值处理

keepapp &＃61; 1;%近似系数不作处理

denoisexs &＃61; wdencmp(&＃39;gbl&＃39;,c,l,wname,lev,thr,&＃39;s&＃39;,keepapp);

denoisexh &＃61; wdencmp(&＃39;gbl&＃39;,c,l,wname,lev,thr,&＃39;h&＃39;,keepapp);

sigOut1 &＃61; WDEN(noisez, N);

sigOut2 &＃61; wden(noisez, &＃39;sqtwolog&＃39;, &＃39;s&＃39;, &＃39;sln&＃39;, 4, &＃39;sym8&＃39;)&＃39;;

% 作图

plot(noisez)

hold on

plot(denoisexs)

plot(denoisexh)

plot(sigOut1)

plot(sigOut2)

hold off

legend(&＃39;原始噪声信号&＃39;, &＃39;matlab软阈值去噪信号&＃39;, &＃39;matlab硬阈值去噪信号&＃39;, &＃39;WDEN&＃39;, &＃39;wden&＃39;);

函数定义如下&＃xff1a;

function sigOut &＃61; WDEN(sigIn, sigLen)

% sym8小波&＃xff0c;4层全局软阈值滤波

[Lo_D, Hi_D, Lo_R, Hi_R] &＃61; wfilters(&＃39;sym8&＃39;); % 滤波器参数

filterLen &＃61; length(Lo_D); % 滤波器长度

Scale &＃61; 4;

srcLen &＃61; sigLen;

msgLen &＃61; zeros(Scale &＃43; 2, 1);

msgLen(1) &＃61; srcLen;

for i &＃61; 2 : Scale &＃43; 1

exLen &＃61; floor((srcLen &＃43; filterLen - 1) / 2);

srcLen &＃61; exLen;

msgLen(i) &＃61; exLen;

end

msgLen(Scale &＃43; 2) &＃61; srcLen;

allSize &＃61; sum(msgLen(2 : end));

dstCoef &＃61; WaveDec(sigIn, msgLen, allSize, Scale, Lo_D, Hi_D, filterLen);

pDet &＃61; dstCoef(allSize - msgLen(2) &＃43; 1 : allSize);

thr &＃61; GetThr(pDet, msgLen);

dstCoef &＃61; Wthresh(dstCoef, thr, allSize, msgLen(Scale &＃43; 1));

sigOut &＃61; WaveRec(dstCoef, msgLen, Scale, Lo_R, Hi_R, filterLen);

end

function dstCoef &＃61; WaveDec(srcData, msgLen, allSize, Scale, Lo_D, Hi_D, filterLen)

tmpSrc &＃61; srcData;

gap &＃61; msgLen(2) * 2;

dstCoef &＃61; zeros(allSize, 1);

for i &＃61; 1 : Scale

curSigLen &＃61; msgLen(i);

tmpDst &＃61; DWT(tmpSrc, curSigLen, Lo_D, Hi_D, filterLen);

dstCoef(allSize - gap &＃43; 1 : allSize - gap &＃43; 2 * msgLen(i &＃43; 1)) &＃61; tmpDst(1 : 2 * msgLen(i &＃43; 1));

tmpSrc(1 : msgLen(i &＃43; 1)) &＃61; tmpDst(1 : msgLen(i &＃43; 1));

gap &＃61; gap - msgLen(i &＃43; 1);

gap &＃61; gap &＃43; 2 * msgLen(i &＃43; 2);

end

end

function dstCoef &＃61; DWT(srcData, srcLen, Lo_D, Hi_D, filterLen)

decLen &＃61; floor((srcLen &＃43; filterLen - 1)/2);

dstCoef &＃61; zeros(2 * decLen, 1);

for i &＃61; 0 : decLen - 1

for j &＃61; 0 : filterLen - 1

k &＃61; 2 * i - j &＃43; 1;

if k <0 && k >&＃61; -filterLen &＃43; 1

tmp &＃61; srcData(-k);

elseif k >&＃61; 0 && k <&＃61; srcLen - 1

tmp &＃61; srcData(k &＃43; 1);

elseif k > srcLen - 1 && k <&＃61; srcLen &＃43; filterLen - 2

tmp &＃61; srcData(2 * srcLen - k);

else

tmp &＃61; 0;

end

dstCoef(i &＃43; 1) &＃61; dstCoef(i &＃43; 1) &＃43; Lo_D(j &＃43; 1) * tmp;

dstCoef(i &＃43; 1 &＃43; decLen) &＃61; dstCoef(i &＃43; 1 &＃43; decLen) &＃43; Hi_D(j &＃43; 1) * tmp;

end

end

end

function thr &＃61; GetThr(detCoef, msgLen)

detCoef &＃61; abs(detCoef);

% detCoef &＃61; sort(detCoef);

sigma &＃61; median(detCoef)/0.6745;

thr &＃61; sigma * sqrt(2 * log(msgLen(1)));

end

function dstCoef &＃61; Wthresh(dstCoef, thr, allSize, gap)

for i &＃61; gap : allSize - 1

if abs(dstCoef(i &＃43; 1))

dstCoef(i &＃43; 1) &＃61; 0;

else

if dstCoef(i &＃43; 1) <0

dstCoef(i &＃43; 1) &＃61; thr - abs(dstCoef(i &＃43; 1));

else

dstCoef(i &＃43; 1) &＃61; abs(dstCoef(i &＃43; 1)) - thr;

end

end

end

end

function dstData &＃61; WaveRec(srcCoef, msgLen, Scale, Lo_R, Hi_R, filterLen)

tmpSrcCoef &＃61; zeros(2 * msgLen(2), 1);

tmpSrcCoef(1 : 2 * msgLen(Scale &＃43; 1)) &＃61; srcCoef(1 : 2 * msgLen(Scale &＃43; 1));

gap &＃61; 2 * msgLen(Scale &＃43; 1);

for i &＃61; Scale : -1 : 1

curDstLen &＃61; msgLen(i);

dstData &＃61; IDWT(tmpSrcCoef, curDstLen, Lo_R, Hi_R, filterLen);

if i ~&＃61; 1

tmpSrcCoef(1 : curDstLen) &＃61; dstData(1 : curDstLen);

tmpSrcCoef(curDstLen &＃43; 1 : 2 * curDstLen) &＃61; srcCoef(gap &＃43; 1 : gap &＃43; curDstLen);

gap &＃61; gap &＃43; msgLen(i);

end

end

end

function recData &＃61; IDWT(srcCoef, dstLen, Lo_R, Hi_R, filterLen)

recLen &＃61; floor((dstLen &＃43; filterLen - 1) / 2);

recData &＃61; zeros(dstLen, 1);

for i &＃61; 0 : dstLen - 1

recData(i &＃43; 1) &＃61; 0;

for j &＃61; 0 : recLen - 1

k &＃61; i - 2 * j &＃43; filterLen - 2;

if k >&＃61; 0 && k

recData(i &＃43; 1) &＃61; recData(i &＃43; 1) &＃43; Lo_R(k &＃43; 1) * srcCoef(j &＃43; 1) &＃43; Hi_R(k &＃43; 1) * srcCoef(j &＃43; 1 &＃43; recLen);

end

end

end

end