小波變換語音消噪改進閾值

1、改進閾值函數(shù)進行語音信號消噪，但是在程序運行過程中頻頻報錯。本人經(jīng)驗不足調(diào)試不出，希望求得各位指導 改進函數(shù)表達式附圖 clear all; clc; close all;祠音信號采樣頻率為 8000fs=8000;xx=wavread('lw1.wav');x1=xx(:,1);% 取單聲道t=(0:le ngth(x1)-1)/8000;%寸信號做2048點FFT變換y1=fft(x1,2048);f=fs*(0:1023)/2048;% 做原始語音信號的時域圖形%加10db的高斯白噪聲 求得信噪比均方誤差figure(1)plot(t,x1)y=awg n( x1'

2、;,10,'measured');sn r,mse=s nrmse(x1,y')% 做加噪語音信號的時域圖形 多尺度一維分解3層分解并提取系數(shù)figure(2)plot(t,y)c,l=wavedec(y,3,'db1');%用 db1小波對信號進行a3=appcoef(c,l,'db1',3);%a2=appcoef(c,l,'db1',2);%a1=appcoef(c,l,'db1',1);d3=detcoef(c,l,3);d2=detcoef(c,l,2);d1=detcoef(c,l,1);thr

3、仁thselect(d1,'rigrsure');%閾值獲取，使用 Stein的無偏風險估計原理thr2=thselect(d2,'rigrsure');thr3=thselect(d3,'rigrsure');%利用改進閾值函數(shù)進行去噪處理gd 仁 Garrote_gg(d1,thr1);gd2=Garrote_gg(d2,thr2);gd3=Garrote_gg(d3,thr3);c1=a3 gd3 gd2 gd1;y1=waverec(c2,l,'db1');%多尺度重構(gòu)sn r,mse=s nrmse(x1,y1'

4、)%求得信噪比均方誤差figure(3);plot(t,y1);為信號分解后的小波系數(shù)，b為獲得的閾值function gd=Garrote_gg(a,b)%a m=0.2*(a*a)-(b*b);if (abs(a)>=b)gd=sig n( a)*(abs(a)-b/exp(m); else (abs(a)<b)gd=0;end function sn r,mse=s nrmse(l,l n)%計算信噪比函數(shù)% I :原始信號% In:去噪后信號snr=0;Ps=sum(sum(l-mea n( mea n(l)92);%sig nal powerPn=sum(sum(l-l

5、n)42);%no ise powersn r=10*log10(Ps/P n);mse=P n/len gth(l);QC截圖 20130516175535.png (11.18 KB, 下載次數(shù)：0)改進函數(shù)表達式IIWx Jsign(W) (|W| _ 麗為盯 |W| > TA ( 0|W| < T唄為閾值處理后的小波系數(shù)；w信號分解后的小波系數(shù)；T為獲得的閾值；M值已在函數(shù)中賦值為02本帖最后由 羅志雄 于2013-5-16 21:58 編輯function sn r,mse=s nrmse(l,l n)%計算信噪比函數(shù)%l :原始信號% ln:去噪后信號snr=0;Ps=

6、sum(sum(l-mea n( mea n(l)92);%sig nal power Pn=sum(sum(l-l n)42);%no ise powersn r=10*log10(Ps/P n);mse=P n/len gth(l);修改后程序清單如下：clear all; clc; close all;fs=8000;祠音信號米樣頻率為 8000xx=wavread('lw1.wav');x1=xx(:,1);%取單聲道x1=x1-mea n( x1);t=(0:le ngth(x1)-1)/8000;y1=fft(x1,2048);%寸信號做2048點FFT變換f=fs

7、*(0:1023)/2048;figure(1)plot(t,x1)% 做原始語音信號的時域圖形y=awg n( x1',10,'measured');%加10db的高斯白噪聲sn r,mse=s nrmsel(x1',y)%求得信噪比均方誤差sn r1=SNR si nglech(x1',y)figure(2)plot(t,y)% 做加噪語音信號的時域圖形c,l=wavedec(y,3,'db1');%多尺度一維分解%用 db1小波對信號進行3層分解并提取系數(shù)a3=appcoef(c,l,'db1',3);%a2=app

8、coef(c,l,'db1',2);%a1=appcoef(c,l,'db1',1);d3=detcoef(c,l,3);d2=detcoef(c,l,2);d1=detcoef(c,l,1);thr1=thselect(d1,'rigrsure');%閾值獲取，使用 Stein的無偏風險估計原理thr2=thselect(d2,'rigrsure');thr3=thselect(d3,'rigrsure');%利用改進閾值函數(shù)進行去噪處理gd 仁 Garrote gg(d1,thr1);gd2=Garrote_g

9、g(d2,thr2);gd3=Garrote gg(d3,thr3);c1=a3 gd3 gd2 gd1;function gd=Garrote_gg(a,b)%a為信號分解后的小波系數(shù)，b為獲得的閾值m=0.2*(a.*a)-(b*b);if (abs(a)>=b)gd=sig n( a)*(abs(a)_b/exp(m);elsegd=zeros(size(a);多尺度重構(gòu)endy1=waverec(c1,l,'db1');%sn r,mse=s nrmsel(x1',y1) %求得信噪比均方誤差figure(3); plot(t,y1);小波去噪軟閾值和硬閾

10、值的matlab仿真程序硬閾值、軟閾值這里有一段不知道有用沒 %設置信噪比和隨機種子值 snr=4;in it=2055615866;%產(chǎn)生原始信號sref和高斯白噪聲污染的信號ssref,s=w no ise(1,11,s nr,i nit);%用db1小波對原始信號進行3層分解并提取系數(shù) c,l=wavedec(s,3,'db1'); a3=appcoef(c,l,'db1',3); d3=detcoef(c,l,3); d2=detcoef(c,l,2); d1=detcoef(c,l,1);thr=1;%進行硬閾值處理ythard 1= wthresh(

11、d1,'h',thr); ythard2=wthresh(d2,'h',thr); ythard3=wthresh(d3,'h',thr); c2=a3 ythard3 ythard2 ythard1; s3=waverec(c2,l,'db1');%進行軟閾值處理ytsoftd仁 wthresh(d1,'s',thr); ytsoftd2=wthresh(d2,'s',thr); ytsoftd3=wthresh(d3,'s',thr); c3=a3 ytsoftd3 ytsoft

12、d2 ytsoftd1; s4=waverec(c3,l,'db1');參考信號'); 染噪信號');硬閾值處理'); 軟閾值處理');%對上述信號進行圖示 subplot(5,1,1);plot(sref);title(' subplot(5,1,2);plot(s);title(' subplot(5,1,3);plot(s3);title(' subplot(5,1,4);plot(s4);title('matlab小波除噪，為何硬閾值和軟閾值除躁信噪比一樣了?matlab小波除噪，為何硬閾值和軟閾值除躁信噪

13、比一樣了？ load leleccum;in dex=1:1024;f1=leleccu m(in dex); %產(chǎn)生含噪信號in it=2055615866;randn ('seed',i nit); f2=f1+18*ra ndn (size(x);sn r=SNR_si nglech(f1,f2) % subplot(2,2,1);plot(f1);title(' subplot(2,2,2);plot(f2);title('%用 db5小波對原始信號進行c,l=wavedec(f2,3,'db6'); a3=appcoef(c,l,

14、9;db6',3); d3=detcoef(c,l,3);d2=detcoef(c,l,2);d1=detcoef(c,l,1);sigma=w no isest(c,l,1); thr=wbmpe n( c,l,sigma,2);%進行硬閾值處理ythard 1= wthresh(d1,'h',thr); ythard2=wthresh(d2,'h',thr); ythard3=wthresh(d3,'h',thr); c2=a3 ythard3 ythard2 ythardl; f3=waverec(c2,l,'db6'

15、;);%進行軟閾值處理ytsoftd仁 wthresh(d1,'s',thr); ytsoftd2=wthresh(d2,'s',thr); ytsoftd3=wthresh(d3,'s',thr); c3=a3 ytsoftd3 ytsoftd2 ytsoftdl; f4=waverec(c3,l,'db6');%對上述信號進行圖示subplot(2,2,3);plot(f3);title(' subplot(2,2,4);plot(f4);title(' sn r=SNR_si nglech(f1,f3) sn

16、 r=SNR_si nglech(f1,f4)信噪比含噪信號');%axis(1,1024,-1,1);含噪信號');%axis(1,1024,-1,1);3層分解并提取系數(shù)硬閾值處理');%axis(1,1024,-1,1);軟閾值處理');%axis(1,1024,-1,1);信噪比函數(shù) SNR_singlech(l,ln) function sn r=SNR_s in glech(l,l n)%計算信噪比函數(shù)% I : riginal signal% ln:no isy sig nal(ie. origi nal sig nal + no ise sig

17、nal)snr=O;Ps=sum(sum(l-mea n( mea n(l)92);%sig nal power Pn=sum(sum(l-l n)42);%no ise powersn r=10*log10(Ps/P n);小波去噪程序Matlab小波去噪(默認，強制，給定三種情況)%利用小波分析對監(jiān)測采集的信號進行去噪處理，恢復原始信號%、波分析進行去噪有 3中方法：%1默認閾值去噪處理。該方法利用函數(shù)ddencmp()生成信號的默認閾值，然后利用函數(shù)wdencmp()進行去噪%2給定閾值去噪處理。在實際的去噪處理過程中，閾值往往可通過經(jīng)驗公式獲得，且這種閾值比默認閾值的可 進行閾值量化處

18、理時可利用函數(shù)wthresh();%3強制去噪處理。該方法是將小波分解結(jié)構(gòu)中的高頻系數(shù)全部置0，即濾掉所有高頻部分，然后對信號進行小種方法比較簡單，且去噪后的信號比較平滑，但是容易丟失信號中的有用成分。%利用小波分析對監(jiān)測采集的水輪機信號進行去噪處理，恢復原始信號%Program Start%載入監(jiān)測所得信號load default.txt;%裝載采集的信號x= default;lx=le ngth(x);t=0:1:le ngth(x)-1'%繪制監(jiān)測所得信號subplot(2,2,1);plot(t,x);title(' 原始信號');grid onset(gcf,

19、'color','w')set(gca,'fo ntn ame','times New Roma n')set(gca,'fo ntsize',14.0)%用db1小波對原始信號進行 3層分解并提取小波系數(shù)c,l=wavedec(x,3,'db1');%sym8 ca3=appcoef(c,l,'db1',3);%低頻部分cd3=detcoef(c,l,3);% 高頻部分 cd2=detcoef(c,l,2);% 高頻部分 cd1=detcoef(c,l,1);% 高頻部分%對信號進

20、行強制去噪處理并圖示cdd3=zeros(1,le ngth(cd3);cdd2=zeros(1,le ngth(cd2);cdd1=zeros(1,le ngth(cd1);c1=ca3,cdd3,cdd2,cdd1;x1=waverec(c1,1,'db1');subplot(2,2,2);plot(xl);title('強制去噪后信號');grid onset(gcf,'color','w')set(gca,'fo ntn ame','times New Roma n')set(gca,

21、9;fo ntsize',14.O)%默認閾值對信號去噪并圖示%用 ddencmp()函數(shù)獲得信號的默認閾值，使用wdencmp()函數(shù)實現(xiàn)去噪過程thr,sorh,keepapp=dde ncmp('de n','wv',x);x2=wde ncmp('gbl',c,l,'db1',3,thr,sorh,keepapp);subplot(2,2,3);plot(x2);title('默認閾值去噪后信號)；grid onset(gcf,'color','w')set(gca,'

22、;fo ntn ame','times New Roma n')set(gca,'fo ntsize',14.0)%給定的軟閾值進行去噪處理并圖示cd1soft=wthresh(cd1,'x',1.465);%經(jīng)驗給出軟閾值數(shù)cd2soft=wthresh(cd2,'x',1.823); %經(jīng)驗給出軟閾值數(shù)cd3soft=wthresh(cd3,'x',2.768); %經(jīng)驗給出軟閾值數(shù)c2=ca3,cd3soft,cd2soft,cd1soft;x3=waverec(c2,1,'db1')

23、;subplot(2,2,4);plot(x3);title('給定軟閾值去噪后信號');grid onset(gcf,'color','w')set(gca,'fo ntn ame','times New Roma n')set(gca,'fo ntsize',14.0)以上就是三種小波去噪的原程序。但紅色標注的地方，我卻運行不過去。提示分別問：1， ? Error using => horzcatCAT argume nts dime nsions are not con siste nt

24、.Error in => tt at 33c1=ca3,cdd3,cdd2,cdd1;2， ? Error using => wthresh at 30In valid argume nt value.Error in => tt at 54cd1soft=wthresh(cd1,'x',1.465);%經(jīng)驗給出軟閾值數(shù)請問為什么？ ？時間域關于小波新閾值去噪改進方法里，matlab程序編寫的討論新閾值能有效地克服軟閾值去噪方法中由于估計值與真實值之間的恒定偏差而 帶來的去噪誤差，也能有效地抑制硬閾值去噪方法中易產(chǎn)生的信號振蕩現(xiàn)象。 但不懂它程序的編寫。請教高

25、手幫忙解決一下！謝謝了！function y = wthresh(x,sorh,t)%WTHRESH Perform soft or hard thresholdi ng.% Y = WTHRESH(X,SORH,T) returns soft (if SORH = 's')% or hard (if SORH = 'h') T-thresholdi ngof the in put% vector or matrix X. T is the threshold value.% Y = WTHRESH(X,'s',T) returns Y = SI

26、GN(X).(|X|-T)+, soft% thresholding is shrinkage.% Y = WTHRESH(X,'h',T) returns Y = X.1_(|X|>T), hard% thresholding is cruder.% See also WDEN, WDENCMP, WPDENCMP.% M. Misiti, Y. Misiti, G. Oppenheim, J.M. Poggi 12-Mar-96.% Last Revisi on: 24-Jul-2007.% Copyright 1995-2007 The MathWorks, Inc

27、.%$Revisi on: 1.11.4.3 $ switch sorhcase 's'tmp = (abs(x)-t);tmp = (tmp+abs(tmp)/2;y = sig n( x).*tmp;case 'h'y = x.*(abs(x)>t);otherwiseerror('Wavelet:Fu nctio nArgVal:l nvalid_ArgVal', .'In valid argume nt value.')end這是傳統(tǒng)的閾值去噪法，只要你有新的算法，直接改就行了，但是我對該法有些疑惑， 珍法算出了每個

28、尺度(假設J=3)的閾值，那么還用wthresh ( x,sorh,thrthr=thr1,thr2,thr3 呢？wthresh ( x,sorh,thr )去噪時，thr 是不是應該求助高手：Matlab小波去噪中遇見的問題小波程序如下:fun ctio n y=f1(x);X=xlsread('shoulia n.xls'); %讀入表格數(shù)據(jù)x=X(:,1); % 將表格數(shù)據(jù)的第一列賦值給 x軸 y=X(:,2); % 將表格數(shù)據(jù)的第二列賦值給 y軸subplot(2,2,1);plot(x,y,'-m*') % 畫出圖形 hold on;title(&#

29、39;原始圖形');grid;%畫出網(wǎng)格%用db4小波對原始信號進行 c,l=wavedec(y,3,'db4'); % ca3=appcoef(c,l,'db4',3); % cd3=detcoef(c,l,3);cd2=detcoef(c,l,2);3層分解并提取系數(shù) 對y進行3層分解 提取系數(shù)還沒有寫好：要增加高頻系數(shù)圖cd1= detcoef(c,l,1);%對信號進行強制性去噪處理并圖示結(jié)果cdd3=zeros(1,le ngth(cd3); cdd2=zeros(1,le ngth(cd2); cdd1=zeros(1,le ngth(cdl

30、); c1=ca3' cdd3 cdd2 cddl; s1=waverec(c1,l,'db4');subplot(2,2,2);plot(x,s1,'-m*');hold on;title('強制去噪后的信號');grid;%用默認閾值對信號進行去噪處理并圖示結(jié)果%用 ddencmp()函數(shù)獲得信號的默認閾值,使用wdencmp()命令函數(shù)實現(xiàn)去噪過程thr,sorh,keepapp=dde ncmp('de n','wv',y);s2=wde ncmp('gbl',c,l,'db

31、4',3,thr,sorh,keepapp);subplot(2,2,3);plot(x,s2,'-m*');hold on;title('默認閾值去噪后的信號');grid;%用給定的軟閾值進行去噪處理cdd1soft=wthresh(cd1,'s',1.465);cdd2soft=wthresh(cd2,'s',1.823);cdd3soft=wthresh(cd3,'s',2.768);c2=ca3 cdd3soft cdd2soft cdd1soft;s3=waverec(c2,l,'db4

32、');subplot(2,2,4);plot(x,s3,'-m*');title('給定軟閾值去噪后的信號');gridclear;錯誤顯示：? Error using => horzcatAll matrices on a row in the bracketed expression must have thesame nu mber of rows.Error in => xnn at 41c2=ca3 cdd3soft cdd2soft cdd1soft;i以至于最后一個圖形沒法顯示同時處理后的圖形數(shù)據(jù)怎么調(diào)出來，以及第一個數(shù)據(jù)應該是

33、【0, 0】 但是處理后的不是。請問是我選擇的分解請問如何修改原貽圖形D 20406080強制去噪后的信號6420默認閾值去噪后的信號MatlabtfifciLoveMatlab.cai需要說明的是：收斂表格的兩組數(shù)據(jù)為:00.0020.4940.9451.0361.0371.48101.49111.86141.96191.78262.75282.78312.71342.54392.85472.71492.71545.43615.47685.47755.20請教程序問題，心電信號去噪的matlab程序想把這個程序。改為只用軟硬閾值對比的心電信號去噪分析，該怎么辦？求高手。程序見下:%(信%波分

34、解%基)于%Haar 小波 %ecg=fope n( '100.dat',丫')；N=1201;data=fread(ecg,N,'i nt16');data=data/10000;save ECGdata data;fclose(ecg);x=0:0.004:4.8;figure(1);subplot(321);plot(x,data);axis(0 4.8 -4 4);title(' 原始心電信號');x=data;wn ame='Haar'level=5;c,l=wavedec(x,level,w name);a5=

35、wrcoef('a',c,l,'Haar',5);a4=wrcoef('a',c,l,'Haar',4);a3=wrcoef('a',c,l,'Haar',3);a2=wrcoef('a',c,l,'Haar',2);a1=wrcoef('a',c,l,'Haar',1);subplot(322);plot(a5);title('a5');axis(0 1201 -2 2);subplot(323);plot(a4);t

36、itle('a4');axis(0 1201 -2 2);subplot(324);plot(a3);title('a3');axis(0 1201 -2 2);subplot(325);plot(a2);title('a2');axis(0 1201 -2 2);subplot(326);plot(a1);title('a1');axis(0 1201 -2 2);%6 小波 %ecg=fope n( '100.dat',丫')；N=1201;data=fread(ecg,N,'i nt16

37、9;);data=data/10;fclose(ecg);x=0:0.004:4.8;figure(2);subplot(321);plot(x,data);axis(0 4.8 -4000 4000);title(' 原始心電信號');x=data;wn ame='db6'level=5;c,l=wavedec(x,level,w name);a5=wrcoef('a',c,l,'db6',5);a4=wrcoef('a',c,l,'db6',4);a3=wrcoef('a',c,

38、l,'db6',3);a2=wrcoef('a',c,l,'db6',2);a1=wrcoef('a',c,l,'db6',1);subplot(322);plot(a5);title('a5');axis(0 1201 -2000 2000);subplot(323);plot(a4);title('a4');axis(0 1201 -2000 2000);subplot(324);plot(a3);title('a3');axis(0 1201 -2000 200

39、0);subplot(325);plot(a2);title('a2');axis(0 1201 -2000 2000);subplot(326);plot(a1);title('a1');axis(0 1201 -2000 2000);%基于%小波 %ecg=fope n( '100.dat',丫')；N=1201;data=fread(ecg,N,'i nt16');data=data/10;save ECGdata data;fclose(ecg);x=0:0.004:4.8;figure(3);subplot(32

40、1);plot(x,data);axis(0 4.8 -4000 4000);title(' 原始心電信號');x=data;wn ame='sym3'level=5;c,l=wavedec(x,level,w name);a5=wrcoef('a',c,l,'sym3',5);a4=wrcoef('a',c,l,'sym3',4);a3=wrcoef('a',c,l,'sym3',3);a2=wrcoef('a',c,l,'sym3'

41、,2);a1=wrcoef('a',c,l,'sym3',1);subplot(322);plot(a5);title('a5');axis(0 1201 -2000 2000);subplot(323);plot(a4);title('a4');axis(0 1201 -2000 2000);subplot(324);plot(a3);title('a3');axis(0 1201 -2000 2000);subplot(325);plot(a2);title('a2');axis(0 1201

42、-2000 2000);subplot(326);plot(a1);title('a1');axis(0 1201 -2000 2000);%心電信號噪%Birge-M策略閾值降噪%基于小波變換的心電信號的降噪ecg=fope n( '100.dat',丫')；N=1201;data=fread(ecg,N,'i nt16');data=data/10000;fclose(ecg);x=data;wave name='db5'level=4;c,l=wavedec(x,level,wave name);alpha=1.5;

43、sorh='h'thr, nkeep=wdcbm(c,l,alpha);%使用硬閾值給信號降噪xc,cxc,lxc,perf0,perfl2=wde ncmp('lvd',c,l,wave name,level,thr,sorh);t1=0:0.004:(le ngth(x)-1)*0.004;figure (4);subplot(211);plot(t1,x);title('從人體采集的原始的ECG言號');subplot(212);plot(t1,xc);title('Birge-Massart策略閾值降噪后的ECG言號(wname=

44、db5 level=4 )');%最優(yōu)%測軟閾值降噪%基于小波變換的心電信號的壓縮ecg=fope n( '100.dat','r');N=1201;data=fread(ecg,N,'i nt16');data=data/10000;fclose(ecg);x=data;wave name='db3'level=4;xd,cxd,lxd=wde n(x,'heursure','s','ml n',level,wave name);t=0:0.004:4.8;figure(5

45、);subplot(311);plot(t,x);title('從人體采集的原始的ECG言號');subplot(312);plot(t,xd);title('heursure規(guī)則閾值降噪后的ECG言號(db3 level=4)');%最小均方差軟閾值降噪何于小波變換的心電信號的壓縮ecg=fope n( '100.dat',丫')；N=1201;data=fread(ecg,N,'i nt16');data=data/10000;fclose(ecg);x=data;wave name='db3'leve

46、l=4;xd=wde n(x,'mi ni maxi','s','ml n',level,wave name);t=0:0.004:4.8;figure(6);subplot(311);plot(t,x);title('從人體采集的原始的ECG言號');subplot(313);plot(t,xd);title('Minimax規(guī)則閾值降噪后的ECG言號(wname=db3 level=4)');%心電信號壓縮%基于小波變換的心電信號的壓縮ecg=fope n( '1OO.dat',丫')；

47、N=1201;data=fread(ecg,N,'i nt16');data=data/10000;fclose(ecg);x=data;wave name='db5'c,l=wavedec(x,4,wave name);alpha=2;sorh='h'thr, nkeep=wdcbm(c,l,alpha);%使用硬閾值壓縮信號xc,cxc,lxc,perf0,perfl2=wde ncmp('lvd',c,l,wave name,4,thr,sorh);t=0:0.004:4.8;figure(7);subplot(511);p

48、lot(t,x);title(' 原始的ECG言號');axis(0 5 -3 3);subplot(512);plot(t,xc);axis(0 5 -3 3);title(' 在第4層對高頻系數(shù)量化壓縮后的ECG言號(alpha=2 db5 )');c,l=wavedec(x,4,wave name);alpha=4;sorh='h'thr, nkeep=wdcbm(c,l,alpha);%使用硬閾值壓縮信號xc1,cxc1,lxc1,perf01,perfl21=wde ncmp('lvd',c,l,wave name,4,

49、thr,sorh); subplot(513);plot(t,xc1);axis(0 5 -3 3);title(' 在第4層對高頻系數(shù)量化壓縮后的ECG言號alpha=4 db5');c,l=wavedec(x,2,wave name);alpha=4;sorh='h'thr, nkeep=wdcbm(c,l,alpha);%使用硬閾值壓縮信號xc1,cxc1,lxc1,perf01,perfl21=wde ncmp('lvd',c,l,wave name,2,thr,sorh); subplot(514);plot(t,xc1);axis(0

50、 5 -3 3);title(' 在第2層對高頻系數(shù)量化壓縮后的ECG言號alpha=4 db5');c,l=wavedec(x,2,'db3');alpha=4;sorh='h'thr, nkeep=wdcbm(c,l,alpha);% 使用硬閾值壓縮信號xc1,cxc1,lxc1,perf01,perfl21=wde ncmp('lvd',c,l,'db3',2,thr,sorh);subplot(515);plot(t,xc1);axis(0 5 -3 3);title(' 在第2層對高頻系數(shù)量化壓縮

51、后的ECG言號alpha=4 db3');%求助解決信號降噪問題求高手幫忙，我的一個實驗數(shù)據(jù)要進行信號降噪，需要用到Birge-Massar閾值降噪，pen alty閾值降噪和缺省的3種降噪方法，并通過比較得出最好的處理方法。如果有那位好心高手幫忙，小弟感激不盡??！請幫忙的高手幫 按LZ的要求程序和得圖如下xx=load('實驗數(shù)據(jù).txt');x=xx(:,1)：y=xx(:,2)：x=x(e nd:-1:1);y=y(e nd:-1:1);%用sym6小波對信號做5層分解 wn ame ='sym6： lev=5;c,l = wavedec(y,lev,w

52、name);%通過第1層的細節(jié)系數(shù)估算信號的噪聲強度 sigma = wno isest(c,l,1);%使用penalty策略確定降噪的閾值alpha=2;thr1 = wbmpe n( c,l,sigma,alpha);%使用Birge-Massart 策略確定降噪的閾值thr2 ,n keep = wdcbm(c,l,alpha);xd1 = wde ncmp('gbl',c,l,w name,lev,thr1,'s',1);%用缺省的閾值確定的時候使用硬閾值時系數(shù)進行處理xd2,cxd,lxd,perf0,perf2= wde ncmp('lvd

53、',c,l,w name,lev,thr2,'h');%求得缺省的閾值 thr,sorh,keepapp = dde ncmp('de n','wv',x);%重建降噪信號figuresubplot(411); plot(y); title('subplot(412); plot(xd1); title('subplot(413); plot(xd2); title('subplot(414); plot(xd2); title('xd3 = wde ncmp('gbl',c,l,w nam

54、e,lev,thr,'s',1);原始信號','fontsize',10);使用penalty閾值降噪后信號','fontsize',10);使用 Birge-Massart閾值降噪后信號','fontsize',10);使用缺省閾值降噪后信號','fon tsize',10);eb1a.jpg (41.07 KB,下載次數(shù)：4)原始信號300250100Mat lab論紜0200400600800200小波去噪：為什么閾值去噪后總把有用信息也給濾出了wave name='E

55、 nglish.wav'% 聲源layer=3;%分解尺度wavelets='db3'%小波選擇SNR_ noise=10;%添加噪音的信噪比TPTR='rigrsure'%閾值選取原則sorh='s'%軟閾值或硬閾值scal='ml n'%閾值每層調(diào)整%導入聲源并添加白噪聲s=wavread(wave name)：y=awg n(s(1,:),SNR_ noise,'measured');%聲源為雙聲道只提取一個聲道分析subplot(211),plot(s(1,:);title('Origi n

56、al Sou nd Wave');subplot(212),plot(y);title('Sou nd Wave with White-Noise');%對含噪聲源進行離散小波分解，并提取高低頻系數(shù)C,L=wavedec(y,layer,wavelets);A=appcoef(C,L,wavelets,layer);D=detcoef(C,L,wavelets,layer);%各層閾值選擇for n=1:layerthre( n)=thselect(D n,TPTR);end%對各系數(shù)進行閾值去噪for n=1:layerD_rec n=wthresh(D n,sorh

57、,thre( n);endfigure;k=0;for n=1:layersubplot(layer,2,k+1),plot(D n);title('Detail.cfs level before thresh in g', nu m2str( n)subplot(layer,2,k+2),plot(D_rec n);title('Detail.cfs level after threshi ng',n um2str( n)k=k+2;end%重構(gòu)信號6= A;for n=1:layerC1=C1,D n;endy_rec=waverec(C1,L,wavele

58、ts);figure;subplot(311),plot(s(1,:);title('Origi nal Sou nd Wave');subplot(312),plot(y);title('Sou nd Wave with White-Noise');subplot(313),plot(y_rec);title('De-Noised Sou nd Wave');為什么用thselect函數(shù)選取的閾值進行濾噪后高頻系數(shù)全為零了我只是想把高頻中的噪聲信號濾出，高手執(zhí)本帖最后由cwjy 于2010-5-2 07:52 編輯1.jpg (117.12 KB,下載次數(shù)：0)Original Sound