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

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

•  源代码

clear;clc;   
load leleccum;  
indx = 1:3450;  
noisez = leleccum(indx);   
wname = 'sym8';   
lev = 4;  
[c,l] = wavedec(noisez,lev,wname);  
[Lo_D,Hi_D,Lo_R,Hi_R] = wfilters(wname);  
% threshold value 
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);  
sigOut1 = WDEN(noisez, N);
sigOut2 = wden(noisez, 'sqtwolog', 's', 'sln', 4, 'sym8')';
% 作图    
plot(noisez) 
hold on
plot(denoisexs) 
plot(denoisexh)
plot(sigOut1)
plot(sigOut2)
hold off
legend('原始噪声信号', 'matlab软阈值去噪信号', 'matlab硬阈值去噪信号', 'WDEN', 'wden');

函数定义如下：
function sigOut = WDEN(sigIn, sigLen)
% sym8小波，4层全局软阈值滤波
[Lo_D, Hi_D, Lo_R, Hi_R] = wfilters('sym8');    % 滤波器参数
filterLen = length(Lo_D);   % 滤波器长度
Scale = 4;
srcLen = sigLen;
msgLen = zeros(Scale + 2, 1);
msgLen(1) = srcLen;
for i = 2 : Scale + 1
    exLen = floor((srcLen + filterLen - 1) / 2);
    srcLen = exLen;
    msgLen(i) = exLen;
end
msgLen(Scale + 2) = srcLen;
allSize = sum(msgLen(2 : end));
dstCoef = WaveDec(sigIn, msgLen, allSize, Scale, Lo_D, Hi_D, filterLen);
pDet = dstCoef(allSize - msgLen(2) + 1 : allSize);
thr = GetThr(pDet, msgLen);
dstCoef = Wthresh(dstCoef, thr, allSize, msgLen(Scale + 1));
sigOut = WaveRec(dstCoef, msgLen, Scale, Lo_R, Hi_R, filterLen);
end

function dstCoef = WaveDec(srcData, msgLen, allSize, Scale, Lo_D, Hi_D, filterLen)
tmpSrc = srcData;
gap = msgLen(2) * 2;
dstCoef = zeros(allSize, 1);
for i = 1 : Scale
    curSigLen = msgLen(i);
    tmpDst = DWT(tmpSrc, curSigLen, Lo_D, Hi_D, filterLen);
    dstCoef(allSize - gap + 1 : allSize - gap + 2 * msgLen(i + 1)) = tmpDst(1 : 2 * msgLen(i + 1));
    tmpSrc(1 : msgLen(i + 1)) = tmpDst(1 : msgLen(i + 1));
    gap = gap - msgLen(i + 1);
    gap = gap + 2 * msgLen(i + 2);
end
end

function dstCoef = DWT(srcData, srcLen, Lo_D, Hi_D, filterLen)
decLen = floor((srcLen + filterLen - 1)/2);
dstCoef = zeros(2 * decLen, 1);
for i = 0 : decLen - 1
    for j = 0 : filterLen - 1
        k = 2 * i - j + 1;
        if k < 0 && k >= -filterLen + 1
            tmp = srcData(-k);
        elseif k >= 0 && k <= srcLen - 1 
            tmp = srcData(k + 1);
        elseif k > srcLen - 1 && k <= srcLen + filterLen - 2 
            tmp = srcData(2 * srcLen - k);
        else
            tmp = 0;
        end
        dstCoef(i + 1) = dstCoef(i + 1) + Lo_D(j + 1) * tmp;
        dstCoef(i + 1 + decLen) = dstCoef(i + 1 + decLen) + Hi_D(j + 1) * tmp;
    end  
end
end

function thr = GetThr(detCoef, msgLen)
detCoef = abs(detCoef);
% detCoef = sort(detCoef);
sigma = median(detCoef)/0.6745;
thr = sigma * sqrt(2 * log(msgLen(1)));
end

function dstCoef = Wthresh(dstCoef, thr, allSize, gap)
for i = gap : allSize - 1
    if abs(dstCoef(i + 1)) < thr
        dstCoef(i + 1) = 0;
    else
        if dstCoef(i + 1) < 0
            dstCoef(i + 1) = thr - abs(dstCoef(i + 1));
        else
           dstCoef(i + 1) = abs(dstCoef(i + 1)) - thr; 
        end
    end
end
end

function dstData = WaveRec(srcCoef, msgLen, Scale, Lo_R, Hi_R, filterLen)
tmpSrcCoef = zeros(2 * msgLen(2), 1);
tmpSrcCoef(1 : 2 * msgLen(Scale + 1)) = srcCoef(1 : 2 * msgLen(Scale + 1));
gap = 2 * msgLen(Scale + 1);
for i = Scale : -1 : 1
    curDstLen = msgLen(i);
    dstData = IDWT(tmpSrcCoef, curDstLen, Lo_R, Hi_R, filterLen);
    if i ~= 1
        tmpSrcCoef(1 : curDstLen) = dstData(1 : curDstLen);
        tmpSrcCoef(curDstLen + 1 : 2 * curDstLen) = srcCoef(gap + 1 : gap + curDstLen);
        gap = gap + msgLen(i);
    end
end
end

function recData = IDWT(srcCoef, dstLen, Lo_R, Hi_R, filterLen)
recLen = floor((dstLen + filterLen - 1) / 2);
recData = zeros(dstLen, 1);
for i = 0 : dstLen - 1
    recData(i + 1) = 0;
    for j = 0 : recLen - 1
        k = i - 2 * j + filterLen - 2;
        if k >= 0 && k < filterLen
            recData(i + 1) = recData(i + 1) + Lo_R(k + 1) * srcCoef(j + 1) + Hi_R(k + 1) * srcCoef(j + 1 + recLen);
        end
    end
end
end