小波變換語(yǔ)音消噪(改進(jìn)閾值)

1、精選優(yōu)質(zhì)文檔-傾情為你奉上改進(jìn)閾值函數(shù)進(jìn)行語(yǔ)音信號(hào)消噪，但是在程序運(yùn)行過程中頻頻報(bào)錯(cuò)。本人經(jīng)驗(yàn)不足調(diào)試不出，希望求得各位指導(dǎo)。改進(jìn)函數(shù)表達(dá)式附圖clear all; clc; close all;fs=8000;                  %語(yǔ)音信號(hào)采樣頻率為8000xx=wavread('lw1.wav');x1=xx(:,1);%取單聲道t=(0:length(x1)-1)/8000;y1=fft(x1,2048);   

2、60;       %對(duì)信號(hào)做2048點(diǎn)FFT變換f=fs*(0:1023)/2048;figure(1)plot(t,x1)                   %做原始語(yǔ)音信號(hào)的時(shí)域圖形y=awgn(x1',10,'measured');   %加10db的高斯白噪聲snr,mse=snrmse(x1,y')%求得信噪比 均方誤差figure(2)plot(t,y)

3、60;                  %做加噪語(yǔ)音信號(hào)的時(shí)域圖形c,l=wavedec(y,3,'db1');%多尺度一維分解%用db1小波對(duì)信號(hào)進(jìn)行3層分解并提取系數(shù)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)

4、;thr1=thselect(d1,'rigrsure');%閾值獲取，使用Stein的無(wú)偏風(fēng)險(xiǎn)估計(jì)原理thr2=thselect(d2,'rigrsure');thr3=thselect(d3,'rigrsure');%利用改進(jìn)閾值函數(shù)進(jìn)行去噪處理gd1=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)snr,mse=snrmse(x1,y1'

5、;)%求得信噪比 均方誤差figure(3);plot(t,y1);function gd=Garrote_gg(a,b)%a為信號(hào)分解后的小波系數(shù)，b為獲得的閾值m=0.2*(a*a)-(b*b);if (abs(a)>=b)    gd=sign(a)*(abs(a)-b/exp(m);else (abs(a)<b)    gd=0;endfunction snr,mse=snrmse(I,In)% 計(jì)算信噪比函數(shù)% I :原始信號(hào)% In:去噪后信號(hào)snr=0;Ps=sum(sum(I-mean(mean(I).2);%signal p

6、owerPn=sum(sum(I-In).2);           %noise powersnr=10*log10(Ps/Pn);mse=Pn/length(I); (11.18 KB, 下載次數(shù): 0)改進(jìn)函數(shù)表達(dá)式本帖最后由 羅志雄 于 2013-5-16 21:58 編輯function snr,mse=snrmse(I,In)% 計(jì)算信噪比函數(shù)% I :原始信號(hào)% In:去噪后信號(hào)snr=0;Ps=sum(sum(I-mean(mean(I).2);%signal powerPn=sum(su

7、m(I-In).2);           %noise powersnr=10*log10(Ps/Pn);mse=Pn/length(I);修改后程序清單如下：clear all; clc; close all;fs=8000;                  %語(yǔ)音信號(hào)采樣頻率為8000xx=wavread('lw1.wav');x1=xx(:,1);%取單聲道x1=x1-m

8、ean(x1);t=(0:length(x1)-1)/8000;y1=fft(x1,2048);           %對(duì)信號(hào)做2048點(diǎn)FFT變換f=fs*(0:1023)/2048;figure(1)plot(t,x1)                   %做原始語(yǔ)音信號(hào)的時(shí)域圖形y=awgn(x1',10,'measured');   %加10d

9、b的高斯白噪聲snr,mse=snrmsel(x1',y)     %求得信噪比 均方誤差snr1=SNR_singlech(x1',y)figure(2)plot(t,y)                   %做加噪語(yǔ)音信號(hào)的時(shí)域圖形c,l=wavedec(y,3,'db1');%多尺度一維分解%用db1小波對(duì)信號(hào)進(jìn)行3層分解并提取系數(shù)a3=appcoef(c,l,'db1',3);%

10、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);thr1=thselect(d1,'rigrsure');%閾值獲取，使用Stein的無(wú)偏風(fēng)險(xiǎn)估計(jì)原理thr2=thselect(d2,'rigrsure');thr3=thselect(d3,'rigrsure');%利用改進(jìn)閾值函數(shù)進(jìn)行去噪處理gd1=Garrote_gg(d1,thr1);gd2=Garro

11、te_gg(d2,thr2);gd3=Garrote_gg(d3,thr3);c1=a3 gd3 gd2 gd1;function gd=Garrote_gg(a,b)%a為信號(hào)分解后的小波系數(shù)，b為獲得的閾值m=0.2*(a.*a)-(b*b);if (abs(a)>=b)    gd=sign(a)*(abs(a)-b/exp(m);else    gd=zeros(size(a);endy1=waverec(c1,l,'db1');%多尺度重構(gòu)snr,mse=snrmsel(x1',y1) %求得信噪比 均方誤差fig

12、ure(3);plot(t,y1);硬閾值、軟閾值這里有一段不知道有用沒%設(shè)置信噪比和隨機(jī)種子值snr=4;init=;%產(chǎn)生原始信號(hào)sref和高斯白噪聲污染的信號(hào)ssref,s=wnoise(1,11,snr,init);%用db1小波對(duì)原始信號(hào)進(jìn)行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;%進(jìn)行硬閾值處理ythard1=wthresh(d1,'h',thr

13、);ythard2=wthresh(d2,'h',thr);ythard3=wthresh(d3,'h',thr);c2=a3 ythard3 ythard2 ythard1;s3=waverec(c2,l,'db1');%進(jìn)行軟閾值處理ytsoftd1=wthresh(d1,'s',thr);ytsoftd2=wthresh(d2,'s',thr);ytsoftd3=wthresh(d3,'s',thr);c3=a3 ytsoftd3 ytsoftd2 ytsoftd1;s4=waverec(c3

14、,l,'db1');%對(duì)上述信號(hào)進(jìn)行圖示subplot(5,1,1);plot(sref);title('參考信號(hào)');subplot(5,1,2);plot(s);title('染噪信號(hào)');subplot(5,1,3);plot(s3);title('硬閾值處理');subplot(5,1,4);plot(s4);title('軟閾值處理'); load leleccum;index=1:1024; f1=leleccum(index); % 產(chǎn)生含噪信號(hào) init=;randn(

15、'seed',init); f2=f1+18*randn(size(x); snr=SNR_singlech(f1,f2) %信噪比subplot(2,2,1);plot(f1);title('含噪信號(hào)'); %axis(1,1024,-1,1);subplot(2,2,2);plot(f2);title('含噪信號(hào)'); %axis(1,1024,-1,1);  %用db5小波對(duì)原始信號(hào)進(jìn)行3層分解并提取系數(shù)c,l=wavedec(f2,3,'db6');a3=appcoef(c,l,&#

16、39;db6',3);d3=detcoef(c,l,3);d2=detcoef(c,l,2);d1=detcoef(c,l,1);sigma=wnoisest(c,l,1);thr=wbmpen(c,l,sigma,2);%進(jìn)行硬閾值處理ythard1=wthresh(d1,'h',thr);ythard2=wthresh(d2,'h',thr);ythard3=wthresh(d3,'h',thr);c2=a3 ythard3 ythard2 ythard1;f3=waverec(c2,l,'db6');%進(jìn)行軟閾值處理

17、ytsoftd1=wthresh(d1,'s',thr);ytsoftd2=wthresh(d2,'s',thr);ytsoftd3=wthresh(d3,'s',thr);c3=a3 ytsoftd3 ytsoftd2 ytsoftd1;f4=waverec(c3,l,'db6');%對(duì)上述信號(hào)進(jìn)行圖示subplot(2,2,3);plot(f3);title('硬閾值處理');%axis(1,1024,-1,1);subplot(2,2,4);plot(f4);title('軟閾值處理');%a

18、xis(1,1024,-1,1);snr=SNR_singlech(f1,f3)snr=SNR_singlech(f1,f4)信噪比函數(shù)SNR_singlech(I,In)function snr=SNR_singlech(I,In)% 計(jì)算信噪比函數(shù)% I：riginal signal% In:noisy signal(ie. original signal + noise signal)snr=0;Ps=sum(sum(I-mean(mean(I).2);%signal powerPn=sum(sum(I-In).2);        &#

19、160;  %noise powersnr=10*log10(Ps/Pn);% 利用小波分析對(duì)監(jiān)測(cè)采集的信號(hào)進(jìn)行去噪處理，恢復(fù)原始信號(hào)%小波分析進(jìn)行去噪有3中方法：%1、默認(rèn)閾值去噪處理。該方法利用函數(shù)ddencmp( )生成信號(hào)的默認(rèn)閾值，然后利用函數(shù)wdencmp( )進(jìn)行去噪處理；%2、給定閾值去噪處理。在實(shí)際的去噪處理過程中，閾值往往可通過經(jīng)驗(yàn)公式獲得，且這種閾值比默認(rèn)閾值的可信度高。在進(jìn)行閾值量化處理時(shí)可利用函數(shù)wthresh( );%3、強(qiáng)制去噪處理。該方法是將小波分解結(jié)構(gòu)中的高頻系數(shù)全部置0，即濾掉所有高頻部分，然后對(duì)信號(hào)進(jìn)行小波重構(gòu)。這種方法比較簡(jiǎn)單，且去噪

20、后的信號(hào)比較平滑，但是容易丟失信號(hào)中的有用成分。% 利用小波分析對(duì)監(jiān)測(cè)采集的水輪機(jī)信號(hào)進(jìn)行去噪處理，恢復(fù)原始信號(hào)%Program Start% 載入監(jiān)測(cè)所得信號(hào)load default.txt;  %裝載采集的信號(hào)x= default;lx=length(x);t=0:1:length(x)-1' % 繪制監(jiān)測(cè)所得信號(hào)subplot(2,2,1);plot(t,x);title('原始信號(hào)');grid onset(gcf,'color','w')  set(gca,'fontna

21、me','times New Roman')set(gca,'fontsize',14.0)% 用db1小波對(duì)原始信號(hào)進(jìn)行3層分解并提取小波系數(shù)c,l=wavedec(x,3,'db1');%sym8ca3=appcoef(c,l,'db1',3);%低頻部分cd3=detcoef(c,l,3);%高頻部分cd2=detcoef(c,l,2);%高頻部分cd1=detcoef(c,l,1);%高頻部分% 對(duì)信號(hào)進(jìn)行強(qiáng)制去噪處理并圖示cdd3=zeros(1,length(cd3);cdd2=zeros(1,length(c

22、d2);cdd1=zeros(1,length(cd1);c1=ca3,cdd3,cdd2,cdd1;x1=waverec(c1,1,'db1');subplot(2,2,2);plot(x1);title('強(qiáng)制去噪后信號(hào)');grid onset(gcf,'color','w')  set(gca,'fontname','times New Roman')set(gca,'fontsize',14.0)% 默認(rèn)閾值對(duì)信號(hào)去噪并圖示%用ddencmp( )函數(shù)獲

23、得信號(hào)的默認(rèn)閾值，使用wdencmp( )函數(shù)實(shí)現(xiàn)去噪過程thr,sorh,keepapp=ddencmp('den','wv',x);x2=wdencmp('gbl',c,l,'db1',3,thr,sorh,keepapp);subplot(2,2,3);plot(x2);title('默認(rèn)閾值去噪后信號(hào)');grid onset(gcf,'color','w')  set(gca,'fontname','times New Roman

24、')set(gca,'fontsize',14.0)% 給定的軟閾值進(jìn)行去噪處理并圖示cd1soft=wthresh(cd1,'x',1.465);%經(jīng)驗(yàn)給出軟閾值數(shù)cd2soft=wthresh(cd2,'x',1.823); %經(jīng)驗(yàn)給出軟閾值數(shù)cd3soft=wthresh(cd3,'x',2.768); %經(jīng)驗(yàn)給出軟閾值數(shù)c2=ca3,cd3soft,cd2soft,cd1soft;x3=waverec(c2,1,'db1');subplot(2,2,4);plot(x3);title('給定

25、軟閾值去噪后信號(hào)');grid onset(gcf,'color','w')  set(gca,'fontname','times New Roman')set(gca,'fontsize',14.0)以上就是三種小波去噪的原程序。但紅色標(biāo)注的地方，我卻運(yùn)行不過去。提示分別問：1，? Error using => horzcatCAT arguments dimensions are not consistent.Error in => tt at 33c1=ca3,cdd3

26、,cdd2,cdd1;2，? Error using => wthresh at 30Invalid argument value.Error in => tt at 54cd1soft=wthresh(cd1,'x',1.465);%經(jīng)驗(yàn)給出軟閾值數(shù)請(qǐng)問為什么？ 新閾值能有效地克服軟閾值去噪方法中由于估計(jì)值與真實(shí)值之間的恒定偏差而帶來的去噪誤差，也能有效地抑制硬閾值去噪方法中易產(chǎn)生的信號(hào)振蕩現(xiàn)象。但不懂它程序的編寫。    請(qǐng)教高手幫忙解決一下！謝謝了！function y = wthresh(x,sorh,t)%WTHRESH Pe

27、rform soft or hard thresholding. %   Y = WTHRESH(X,SORH,T) returns soft (if SORH = 's')%   or hard (if SORH = 'h') T-thresholding  of the input %   vector or matrix X. T is the threshold value.%   Y = WTHRESH(X,'s',T

28、) returns Y = SIGN(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. Pogg

29、i 12-Mar-96.%   Last Revision: 24-Jul-2007.%   Copyright 1995-2007 The MathWorks, Inc.%        $Revision: 1.11.4.3 $switch sorh  case 's'    tmp = (abs(x)-t);    tmp = (tmp+abs(tmp)/2;    y   = sign(x).

30、*tmp;  case 'h'    y   = x.*(abs(x)>t);  otherwise    error('Wavelet:FunctionArgVal:Invalid_ArgVal', .        'Invalid argument value.')end這是傳統(tǒng)的閾值去噪法，只要你有新的算法，直接改就行了，但是我對(duì)該法有些疑惑，wthresh（x,sorh,t

31、hr）,若我用趙瑞珍法算出了每個(gè)尺度（假設(shè)J=3）的閾值，那么還用wthresh（x,sorh,thr）去噪時(shí)，thr是不是應(yīng)該在之前寫成thr=thr1,thr2,thr3呢？小波程序如下：function y=f1(x);X=xlsread('shoulian.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('原始圖形');grid;%畫出網(wǎng)格%用db4小波對(duì)原始信

32、號(hào)進(jìn)行3層分解并提取系數(shù)  還沒有寫好:要增加高頻系數(shù)圖c,l=wavedec(y,3,'db4'); %對(duì)y進(jìn)行3層分解ca3=appcoef(c,l,'db4',3); %提取系數(shù)cd3=detcoef(c,l,3);cd2=detcoef(c,l,2);cd1=detcoef(c,l,1);%對(duì)信號(hào)進(jìn)行強(qiáng)制性去噪處理并圖示結(jié)果cdd3=zeros(1,length(cd3);cdd2=zeros(1,length(cd2);cdd1=zeros(1,length(cd1);c1=ca3' cdd3 cdd2 cdd1;s1=wa

33、verec(c1,l,'db4');subplot(2,2,2);plot(x,s1,'-m*');hold on;title('強(qiáng)制去噪后的信號(hào)');grid;%用默認(rèn)閾值對(duì)信號(hào)進(jìn)行去噪處理并圖示結(jié)果%用ddencmp()函數(shù)獲得信號(hào)的默認(rèn)閾值,使用wdencmp()命令函數(shù)實(shí)現(xiàn)去噪過程thr,sorh,keepapp=ddencmp('den','wv',y);s2=wdencmp('gbl',c,l,'db4',3,thr,sorh,keepapp);subplot(2,2,3

34、);plot(x,s2,'-m*');hold on;title('默認(rèn)閾值去噪后的信號(hào)');grid;%用給定的軟閾值進(jìn)行去噪處理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');subplot(2,2,4);plot(x,s3,'-m*

35、');title('給定軟閾值去噪后的信號(hào)');gridclear;錯(cuò)誤顯示：? Error using => horzcatAll matrices on a row in the bracketed expression must have the same number of rows.Error in => xnn at 41c2=ca3 cdd3soft cdd2soft cdd1soft;i以至于最后一個(gè)圖形沒法顯示同時(shí)處理后的圖形數(shù)據(jù)怎么調(diào)出來，以及第一個(gè)數(shù)據(jù)應(yīng)該是【0，0】，但是處理后的不是。請(qǐng)問是我選擇的分解重構(gòu)函數(shù)部隊(duì)還是其他原

36、因請(qǐng)問如何修改 需要說明的是：收斂表格的兩組數(shù)據(jù)為：00.0020.4940.9451.0361.0371.48101.49111.86141.96191.78262.75282.78312.71342.54392.85472.71492.71545.43615.47685.47755.20想把這個(gè)程序。改為只用軟硬閾值對(duì)比的心電信號(hào)去噪分析，該怎么辦？ 求高手。程序見下：%信號(hào)小波分解%基于Haar小波%ecg=fopen('100.dat','r');N=1201;data=fread(ecg,N,'int16');data=dat

37、a/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('原始心電信號(hào)');x=data;wname='Haar'level=5;c,l=wavedec(x,level,wname);a5=wrcoef('a',c,l,'Haar',5);a4=wrcoef('a',c,l,'Haar',4);a3=wrcoef('a'

38、,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);title('a4');axis(0 1201 -2 2);subplot(324);plot(a3);title('a3');axis(0 1201 -2 2);subplot(3

39、25);plot(a2);title('a2');axis(0 1201 -2 2);subplot(326);plot(a1);title('a1');axis(0 1201 -2 2);%基于db6小波%ecg=fopen('100.dat','r');N=1201;data=fread(ecg,N,'int16');data=data/10;fclose(ecg);x=0:0.004:4.8;figure(2);subplot(321);plot(x,data);axis(0 4.8 -4000 4000);

40、title('原始心電信號(hào)');x=data;wname='db6'level=5;c,l=wavedec(x,level,wname);a5=wrcoef('a',c,l,'db6',5);a4=wrcoef('a',c,l,'db6',4);a3=wrcoef('a',c,l,'db6',3);a2=wrcoef('a',c,l,'db6',2);a1=wrcoef('a',c,l,'db6',1);

41、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 -2000 2000);subplot(326);plot(a1);title('a1&

42、#39;);axis(0 1201 -2000 2000);%基于sym3小波%ecg=fopen('100.dat','r');N=1201;data=fread(ecg,N,'int16');data=data/10;save ECGdata data;fclose(ecg);x=0:0.004:4.8;figure(3);subplot(321);plot(x,data);axis(0 4.8 -4000 4000);title('原始心電信號(hào)');x=data;wname='sym3'level=5;c,l

43、=wavedec(x,level,wname);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',2);a1=wrcoef('a',c,l,'sym3',1);subplot(322);plot(a5);title('a5');axis(0 1201 -2

44、000 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 -2000 2000);subplot(326);plot(a1);title('a1');axis(0 1201 -2000 2000);%心電信號(hào)降噪%Birge-Massart策略閾值降

45、噪%基于小波變換的心電信號(hào)的降噪ecg=fopen('100.dat','r');N=1201;data=fread(ecg,N,'int16');data=data/10000;fclose(ecg);x=data;wavename='db5'level=4;c,l=wavedec(x,level,wavename);alpha=1.5;sorh='h'thr,nkeep=wdcbm(c,l,alpha);%使用硬閾值給信號(hào)降噪xc,cxc,lxc,perf0,perfl2=wdencmp('lvd

46、9;,c,l,wavename,level,thr,sorh);t1=0:0.004:(length(x)-1)*0.004;figure(4);subplot(211);plot(t1,x);title('從人體采集的原始的ECG信號(hào)');subplot(212);plot(t1,xc);title('Birge-Massart策略閾值降噪后的ECG信號(hào)（wname=db5 level=4）');%最優(yōu)預(yù)測(cè)軟閾值降噪%基于小波變換的心電信號(hào)的壓縮ecg=fopen('100.dat','r');N=1201;data=fread

47、(ecg,N,'int16');data=data/10000;fclose(ecg);x=data;wavename='db3'level=4;xd,cxd,lxd=wden(x,'heursure','s','mln',level,wavename);t=0:0.004:4.8;figure(5);subplot(311);plot(t,x);title('從人體采集的原始的ECG信號(hào)');subplot(312);plot(t,xd);title('heursure規(guī)則閾值降噪后的EC

48、G信號(hào)(db3 level=4)');%最小均方差軟閾值降噪%基于小波變換的心電信號(hào)的壓縮ecg=fopen('100.dat','r');N=1201;data=fread(ecg,N,'int16');data=data/10000;fclose(ecg);x=data;wavename='db3'level=4;xd=wden(x,'minimaxi','s','mln',level,wavename);t=0:0.004:4.8;figure(6);subplot(3

49、11);plot(t,x);title('從人體采集的原始的ECG信號(hào)');subplot(313);plot(t,xd);title('Minimax規(guī)則閾值降噪后的ECG信號(hào)(wname=db3 level=4)');%心電信號(hào)壓縮%基于小波變換的心電信號(hào)的壓縮ecg=fopen('100.dat','r');N=1201;data=fread(ecg,N,'int16');data=data/10000;fclose(ecg);x=data;wavename='db5'c,l=wavedec(

50、x,4,wavename);alpha=2;sorh='h'thr,nkeep=wdcbm(c,l,alpha);%使用硬閾值壓縮信號(hào)xc,cxc,lxc,perf0,perfl2=wdencmp('lvd',c,l,wavename,4,thr,sorh);t=0:0.004:4.8;figure(7);subplot(511);plot(t,x);title('原始的ECG信號(hào)');axis(0 5 -3 3);subplot(512);plot(t,xc);axis(0 5 -3 3);title('在第4層對(duì)高頻系數(shù)量化壓縮后的E

51、CG信號(hào)（alpha=2 db5）');c,l=wavedec(x,4,wavename);alpha=4;sorh='h'thr,nkeep=wdcbm(c,l,alpha);%使用硬閾值壓縮信號(hào)xc1,cxc1,lxc1,perf01,perfl21=wdencmp('lvd',c,l,wavename,4,thr,sorh);subplot(513);plot(t,xc1);axis(0 5 -3 3);title('在第4層對(duì)高頻系數(shù)量化壓縮后的ECG信號(hào)alpha=4 db5');c,l=wavedec(x,2,wavename

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

53、,l,alpha);%使用硬閾值壓縮信號(hào)xc1,cxc1,lxc1,perf01,perfl21=wdencmp('lvd',c,l,'db3',2,thr,sorh);subplot(515);plot(t,xc1);axis(0 5 -3 3);title('在第2層對(duì)高頻系數(shù)量化壓縮后的ECG信號(hào)alpha=4 db3');% 求高手幫忙，我的一個(gè)實(shí)驗(yàn)數(shù)據(jù)要進(jìn)行信號(hào)降噪，需要用到Birge-Massar閾值降噪，penalty閾值降噪和缺省的閾值降噪這3種降噪方法，并通過比較得出最好的處理方法。如果有那位好心高手幫忙，小弟感激不盡

54、??！請(qǐng)幫忙的高手幫下忙?。“碙Z的要求程序和得圖如下xx=load('實(shí)驗(yàn)數(shù)據(jù).txt');x=xx(:,1)'y=xx(:,2)'x=x(end:-1:1);y=y(end:-1:1);% 用sym6小波對(duì)信號(hào)做5層分解wname ='sym6' lev=5;c,l = wavedec(y,lev,wname);% 通過第1層的細(xì)節(jié)系數(shù)估算信號(hào)的噪聲強(qiáng)度 sigma = wnoisest(c,l,1);% 使用penalty策略確定降噪的閾值alpha=2;thr1 = wbmpen(c,l,sigma,alpha);% 使用Birg

55、e-Massart策略確定降噪的閾值thr2,nkeep = wdcbm(c,l,alpha);xd1 = wdencmp('gbl',c,l,wname,lev,thr1,'s',1);% 用缺省的閾值確定的時(shí)候使用硬閾值時(shí)系數(shù)進(jìn)行處理xd2,cxd,lxd,perf0,perf2= wdencmp('lvd',c,l,wname,lev,thr2,'h');% 求得缺省的閾值thr,sorh,keepapp = ddencmp('den','wv',x);% 重建降噪信號(hào)xd3 = wdencm

56、p('gbl',c,l,wname,lev,thr,'s',1);figuresubplot(411); plot(y); title('原始信號(hào)','fontsize',10);subplot(412); plot(xd1); title('使用penalty閾值降噪后信號(hào)','fontsize',10);subplot(413); plot(xd2); title('使用Birge-Massart閾值降噪后信號(hào)','fontsize',10);subplot(41

57、4); plot(xd2); title('使用缺省閾值降噪后信號(hào)', 'fontsize',10); (41.07 KB, 下載次數(shù): 4)wavename='English.wav'%聲源layer=3;%分解尺度wavelets='db3'%小波選擇SNR_noise=10;%添加噪音的信噪比TPTR='rigrsure'%閾值選取原則sorh='s'%軟閾值 或 硬閾值scal='mln'%閾值每層調(diào)整%導(dǎo)入聲源 并 添加白噪聲s=wavread(wavename)

58、'y=awgn(s(1,:),SNR_noise,'measured');%聲源為雙聲道 只提取一個(gè)聲道分析subplot(211),plot(s(1,:);title('Original Sound Wave');subplot(212),plot(y);title('Sound Wave with White-Noise');%對(duì)含噪聲源進(jìn)行離散小波分解，并提取高低頻系數(shù)C,L=wavedec(y,layer,wavelets);A=appcoef(C,L,wavelets,layer);D=detcoef(C,L,wavelets,

59、layer);%各層閾值選擇for n=1:layer    thre(n)=thselect(Dn,TPTR);end%對(duì)各系數(shù)進(jìn)行閾值去噪for n=1:layer    D_recn=wthresh(Dn,sorh,thre(n);endfigure;k=0;for n=1:layer    subplot(layer,2,k+1),plot(Dn);    title('Detail.cfs level before threshing',num2str(n)    su

60、bplot(layer,2,k+2),plot(D_recn);    title('Detail.cfs level after threshing',num2str(n)    k=k+2;end%重構(gòu)信號(hào)C1=A;for n=1:layer    C1=C1,Dn;endy_rec=waverec(C1,L,wavelets);figure;subplot(311),plot(s(1,:);title('Original Sound Wave');subplot(312),plot(y);titl

61、e('Sound Wave with White-Noise');subplot(313),plot(y_rec);title('De-Noised Sound Wave');為什么用thselect函數(shù)選取的閾值進(jìn)行濾噪后 高頻系數(shù)全為零了 我只是想把高頻中的噪聲信號(hào)濾出，高手執(zhí)教 本帖最后由 cwjy 于 2010-5-2 07:52 編輯  (117.12 KB, 下載次數(shù): 0) %基于Ephraim和Van Trees提出的信號(hào)子空間法(TDC)的語(yǔ)音增強(qiáng)程序，適用于噪聲為白噪聲的情況%   

62、08.3.22     magic      SCUclear;%-參數(shù)定義-frame_len=40;              %幀長(zhǎng)step_len=0.5*frame_len;    %分幀時(shí)的步長(zhǎng)，相當(dāng)于重疊50%N=8;     %計(jì)算Toeplitz協(xié)方差矩陣時(shí)用到的前后相鄰的幀數(shù),N為偶數(shù)%-讀入帶噪語(yǔ)音文件-filename,pathnam

63、e=uigetfile('*.wav','請(qǐng)選擇語(yǔ)音文件：');wavin,fs,nbits=wavread(pathname filename);wav_length=length(wavin);%取前3000采樣點(diǎn)用于估計(jì)噪聲方差n_var=var(wavin(1:3000),1);xv=n_var;      %定義Rx的特征值判別閾值%-分幀-inframe=(enframe(wavin,frame_len,step_len)'frame_num=size(inframe,2);%-臨時(shí)參數(shù)-%記錄dis_

64、SNR以及max_Ax，查看下面的SNR變量時(shí)使用，在程序最后，對(duì)該變量進(jìn)行制圖，根據(jù)情況使用max_Ax=zeros(1,frame_num);dis_SNR=zeros(1,frame_num);%-子空間法語(yǔ)音增強(qiáng)-%求噪聲的Toeplitz協(xié)方差矩陣，Rn=n_var*I;Rn=n_var*eye(frame_len);Ry=zeros(frame_len,frame_len);          %定義幀信號(hào)的Toeplitz協(xié)方差矩陣seq=zeros(1,frame_len);     &#

65、160;           %定義相鄰6幀的自相關(guān)序列，長(zhǎng)度等于幀長(zhǎng)L=N*frame_len;                          %定義相鄰N幀的長(zhǎng)度outframe=zeros(frame_len,frame_num);    %定義增強(qiáng)后的矩陣%對(duì)每一幀進(jìn)行處理，得到增強(qiáng)后的信號(hào),因?yàn)橛?jì)算Toe

66、plitz協(xié)方差矩陣要用到N幀，其中當(dāng)前%幀處在第(N/2+1)幀（注意：計(jì)算Toeplitz協(xié)方差矩陣的幀是不重疊的，而inframe相鄰幀%之間重疊50%），因此，inframe的前N幀和最后(N/2-1)*2幀無(wú)法估計(jì)Toeplitz協(xié)方差矩陣，%且這(N-1)*2幀一般是噪聲信號(hào)，故不必對(duì)其處理，直接將out_frame的這些幀置0for i=(N+1):(frame_num-N-2)%1.構(gòu)造當(dāng)前幀的Toeplitz協(xié)方差矩陣    for j=0:(frame_len-1)        bgn_point=(i-N-1)*step_len+1;       %相鄰6幀的起始點(diǎn)        end_point=(i+N-1)*step_len;         %相鄰6幀的終點(diǎn)       &