一維信號去噪:提升小波閾值去噪



小波變換是現在研究的比較多的時(空)頻域分析理論,離散的小波變換(DWT)的快速算法時最近研究的熱點.Swelden提出的一種不依賴於傅立葉變換
的新的小波構造方案-----lifting scheme,其復雜度只有原來卷積方法的一半左右,因此成為計算離散小波變換的主流方法.
      其實lifting scheme就是為了構造第二代小波,使得不像第一代小波那樣構造,非常依賴Fourier變換.同時已經證明了提升方式可以實現所有的第一代小波變換.
      提升方式的特點:
1. 繼承了第一代小波的多分辨率的特性
2. 不依賴傅立葉變換
3.不占用系統內存
4.反變換很容易從正變換得到,只是改變了數據流的方向和正負號


本文是基於matlab實現的一維信號的提升小波去噪算法實現。

實現步驟:1.用小波提升方案lifting scheme提升小波;

                  2.對該小波進行多層小波分解;

                  3.對每層的高頻系數進行閾值設定,進行消噪處理;

                  4.一層一層地進行小波重構,得到消噪后的信號。


閾值設置:閾值公式:

:每層小波高頻系數絕對值按照從小到大順序排列后的中值/0.6745;

N:每層小波高頻系數長度值;

說明:1.該閾值公式是Donoho發明的,經過證明后發現噪聲的最大限度在很高的概率上是低於Thr()閾值的,所以可以用該閾值公式得到的閾值來區別噪 聲和信號分量,將二者分離,得到去噪后信號。

2.Thr()這個閾值不是最佳的閾值,而是最佳閾值上限,我們在設置時可以自由設置每層閾值(根據信號特點)。


去噪方法:

這里介紹軟閾值和硬閾值去噪算法,一般硬閾值去噪比軟閾值去噪所得信號去噪效果更粗糙。

小波閾值去噪的具體處理過程:選擇一個小波並確定分解層數,然后將含噪信號在各尺度上進行小波分解,保留大尺度低分辨率下的全部小波系數;

對於各尺度高分辨率下的小波系數,可以設定一個閾值,幅值低於該閾值的小波系數置為0,高於該閾值的小波系數或者完整保留,或者作相應的“收縮(shrinkage)”;

最后將處理后獲得的小波系數利用逆小波變換進行重構,恢復出有效的信號。



為閾值,

f(x)為經閾值處理后的高頻系數,

x為高頻小波系數,







MATLAB程序:

%(程序非原創,只做理解)

%軟閾值去噪

%提升小波分解
nx=b1;%一維含噪信號
lsdb4=liftwave('db4');%采用db4小波,得到相應提升方案
els={'p',[-0.125,0.125],0};
lsnew=addlift(lsdb4,els);%添加els到提升方案
xDEC=lwt(nx,lsnew,3);%用lsnew提升小波對信號做3層小波分解
ca1=lwtcoef('ca',xDEC,lsnew,3,1);%ca1、ca2、ca3為小波分解低頻系數
ca2=lwtcoef('ca',xDEC,lsnew,3,2);
ca3=lwtcoef('ca',xDEC,lsnew,3,3);
cd1=lwtcoef('cd',xDEC,lsnew,3,1);%cd1、cd2、cd3為小波分解高頻系數
cd2=lwtcoef('cd',xDEC,lsnew,3,2);
cd3=lwtcoef('cd',xDEC,lsnew,3,3);
%第一層小波閾值設置
len=length(cd1);%得到一層小波高頻系數長度
w=sort(abs(cd1));%對cd1從小至大進行排序

%找cd1系數中值
if rem(len,2)==1%rem函數功能:求余數;判斷len/2余數是否為1,確定cd1長度是奇是偶
    v=w((len+1)/2);%長度為奇,取cd1系數中間值
else
    v=(w(len/2)+w(len/2+1))/2;%長度為偶,取中間兩個值求平均
end
sigma1=abs(v)/0.6745;%求sigma1
Thr1=sigma1*(2*log(len))^(1/2);%求閾值Thr1

%軟閾值處理系數
for ii=1:length(cd1)%for循環
    if(abs(cd1(ii))<=Thr1)%系數絕對值小於閾值,置0
        cd1(ii)=0;
    else if(cd1(ii)>Thr1)%系數絕對值大於閾值
            cd1(ii)=(cd1(ii)-Thr1);%軟閾值公式,對絕對值和sgn函數處理后得
        else
            cd1(ii)=(cd1(ii)+Thr1);
        end
    end
end
%第二層小波閾值設置      
len=length(cd2);
w=sort(abs(cd2));
if rem(len,2)==1
    v=w((len+1)/2);
else
    v=(w(len/2)+w(len/2+1))/2;
end
sigma1=abs(v)/0.6745;
Thr1=1/2*sigma1*(2*log(len))^(1/2);
for ii=1:length(cd2)
    if(abs(cd2(ii))<=Thr1)
        cd2(ii)=0;
    else if(cd2(ii)>Thr1)
            cd2(ii)=cd2(ii)-Thr1;
        else
            cd2(ii)=cd2(ii)+Thr1;
        end
    end
end
%第三層小波閾值設置
len=length(cd3);
w=sort(abs(cd3));
if rem(len,2)==1
    v=w((len+1)/2);
else
    v=(w(len/2)+w(len/2+1))/2;
end
sigma1=abs(v)/0.6745;
Thr1=2/3*sigma1*(2*log(len))^(1/2);
for ii=1:length(cd3)
    if(abs(cd3(ii))<=Thr1)
        cd3(ii)=0;
    else if(cd3(ii)>Thr1)
            cd3(ii)=cd3(ii)-Thr1;
        else
            cd3(ii)=cd3(ii)+Thr1;
        end
    end
end
%小波重建
rec1=ilwt(ca3,cd3,lsnew);%對小波系數進行一層一層重構
rec2=ilwt(rec1,cd2,lsnew);
rec3=ilwt(rec2,cd1,lsnew);%rec3為去噪后信號

figure(6);
subplot(2,1,2);
plot(rec3);
title('提升小波軟閾值消噪后的信號');grid on;





%硬閾值去噪

更改程序片段為:(每層小波閾值設置時都要更改為下面片段)

for ii=1:length(cd1)
    if(abs(cd1(ii))>=Thr1)
        cd1(ii)=cd1(ii);
    else
        cd1(ii)=0;
    end
end



注意!

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



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