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

1、改進閾值函數(shù)進行語音信號消噪，但是在程序運行過程中頻頻報錯。本人經驗不足調試不出，希望求得各位指導。改進函數(shù)表達式附圖clear all; clc; close all;fs=8000;                  %語音信號采樣頻率為8000xx=wavread('lw1.wav');x1=xx(:,1);%取單聲道t=(0:length(x1)-1)/8000;y1=fft(x1,2048);      

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

3、60;               %做加噪語音信號的時域圖形c,l=wavedec(y,3,'db1');%多尺度一維分解%用db1小波對信號進行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);thr1=thselec

4、t(d1,'rigrsure');%閾值獲取，使用Stein的無偏風險估計原理thr2=thselect(d2,'rigrsure');thr3=thselect(d3,'rigrsure');%利用改進閾值函數(shù)進行去噪處理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');%多尺度重構snr,mse=snrmse(x1,y1')%求得信噪比 均方誤差

5、figure(3);plot(t,y1);function gd=Garrote_gg(a,b)%a為信號分解后的小波系數(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)% 計算信噪比函數(shù)% I :原始信號% In:去噪后信號snr=0;Ps=sum(sum(I-mean(mean(I).2);%signal powerPn=sum(su

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

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

8、length(x1)-1)/8000;y1=fft(x1,2048);           %對信號做2048點FFT變換f=fs*(0:1023)/2048;figure(1)plot(t,x1)                   %做原始語音信號的時域圖形y=awgn(x1',10,'measured');   %加10db的高斯白噪聲snr,ms

9、e=snrmsel(x1',y)     %求得信噪比 均方誤差snr1=SNR_singlech(x1',y)figure(2)plot(t,y)                   %做加噪語音信號的時域圖形c,l=wavedec(y,3,'db1');%多尺度一維分解%用db1小波對信號進行3層分解并提取系數(shù)a3=appcoef(c,l,'db1',3);%a2=appcoef(c,

10、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ù)進行去噪處理gd1=Garrote_gg(d1,thr1);gd2=Garrote_gg(d2,thr2

11、);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=sign(a)*(abs(a)-b/exp(m);else    gd=zeros(size(a);endy1=waverec(c1,l,'db1');%多尺度重構snr,mse=snrmsel(x1',y1) %求得信噪比 均方誤差figure(3);plot(t

12、,y1);小波去噪軟閾值和硬閾值的matlab仿真程序硬閾值、軟閾值這里有一段不知道有用沒%設置信噪比和隨機種子值snr=4;init=2055615866;%產生原始信號sref和高斯白噪聲污染的信號ssref,s=wnoise(1,11,snr,init);%用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;%進行硬閾值處理ythard1=wthresh

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

14、oftd1;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小波除噪，為何硬閾值和軟閾值除躁信噪比一樣了？lo

15、ad leleccum;index=1:1024; f1=leleccum(index); % 產生含噪信號 init=2055615866;randn('seed',init); f2=f1+18*randn(size(x); snr=SNR_singlech(f1,f2) %信噪比subplot(2,2,1);plot(f1);title('含噪信號'); %axis(1,1024,-1,1);subplot(2,2,2);plot(f2);title('含噪信號'); %axis(1,1024,-1,

16、1);  %用db5小波對原始信號進行3層分解并提取系數(shù)c,l=wavedec(f2,3,'db6');a3=appcoef(c,l,'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);%進行硬閾值處理ythard1=wthresh(d1,'h',thr);ythard2=wthresh(d2,'h',thr);ythard3=wthresh(d3

17、,'h',thr);c2=a3 ythard3 ythard2 ythard1;f3=waverec(c2,l,'db6');%進行軟閾值處理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');%對上述信號進行圖示subplot(2,2,3);plot(f3);titl

18、e('硬閾值處理');%axis(1,1024,-1,1);subplot(2,2,4);plot(f4);title('軟閾值處理');%axis(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)% 計算信噪比函數(shù)% I：riginal signal% In:noisy signal(ie. original signal + noise signal)snr=0;Ps=sum(sum

19、(I-mean(mean(I).2);%signal powerPn=sum(sum(I-In).2);           %noise powersnr=10*log10(Ps/Pn);小波去噪程序Matlab小波去噪（默認，強制，給定三種情況）% 利用小波分析對監(jiān)測采集的信號進行去噪處理，恢復原始信號%小波分析進行去噪有3中方法：%1、默認閾值去噪處理。該方法利用函數(shù)ddencmp( )生成信號的默認閾值，然后利用函數(shù)wdencmp( )進行去噪處理；%2、給定閾值去噪處理。在實際的去噪處理過程中，閾值往往可

20、通過經驗公式獲得，且這種閾值比默認閾值的可信度高。在進行閾值量化處理時可利用函數(shù)wthresh( );%3、強制去噪處理。該方法是將小波分解結構中的高頻系數(shù)全部置0，即濾掉所有高頻部分，然后對信號進行小波重構。這種方法比較簡單，且去噪后的信號比較平滑，但是容易丟失信號中的有用成分。% 利用小波分析對監(jiān)測采集的水輪機信號進行去噪處理，恢復原始信號%Program Start% 載入監(jiān)測所得信號load default.txt;  %裝載采集的信號x= default;lx=length(x);t=0:1:length(x)-1' % 繪制監(jiān)測所得信號subp

21、lot(2,2,1);plot(t,x);title('原始信號');grid onset(gcf,'color','w')  set(gca,'fontname','times New Roman')set(gca,'fontsize',14.0)% 用db1小波對原始信號進行3層分解并提取小波系數(shù)c,l=wavedec(x,3,'db1');%sym8ca3=appcoef(c,l,'db1',3);%低頻部分cd3=detcoef(c,l,3

22、);%高頻部分cd2=detcoef(c,l,2);%高頻部分cd1=detcoef(c,l,1);%高頻部分% 對信號進行強制去噪處理并圖示cdd3=zeros(1,length(cd3);cdd2=zeros(1,length(cd2);cdd1=zeros(1,length(cd1);c1=ca3,cdd3,cdd2,cdd1;x1=waverec(c1,1,'db1');subplot(2,2,2);plot(x1);title('強制去噪后信號');grid onset(gcf,'color','w') 

23、60;set(gca,'fontname','times New Roman')set(gca,'fontsize',14.0)% 默認閾值對信號去噪并圖示%用ddencmp( )函數(shù)獲得信號的默認閾值，使用wdencmp( )函數(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('默

24、認閾值去噪后信號');grid onset(gcf,'color','w')  set(gca,'fontname','times New Roman')set(gca,'fontsize',14.0)% 給定的軟閾值進行去噪處理并圖示cd1soft=wthresh(cd1,'x',1.465);%經驗給出軟閾值數(shù)cd2soft=wthresh(cd2,'x',1.823); %經驗給出軟閾值數(shù)cd3soft=wthresh(cd3,'x'

25、,2.768); %經驗給出軟閾值數(shù)c2=ca3,cd3soft,cd2soft,cd1soft;x3=waverec(c2,1,'db1');subplot(2,2,4);plot(x3);title('給定軟閾值去噪后信號');grid onset(gcf,'color','w')  set(gca,'fontname','times New Roman')set(gca,'fontsize',14.0)以上就是三種小波去噪的原程序。但紅色標注的地方，我卻運行

26、不過去。提示分別問：1，? Error using => horzcatCAT arguments dimensions are not consistent.Error in => tt at 33c1=ca3,cdd3,cdd2,cdd1;2，? Error using => wthresh at 30Invalid argument value.Error in => tt at 54cd1soft=wthresh(cd1,'x',1.465);%經驗給出軟閾值數(shù)請問為什么？關于小波新閾值去噪改進方法里，matlab程序編寫的討論 新閾值

27、能有效地克服軟閾值去噪方法中由于估計值與真實值之間的恒定偏差而帶來的去噪誤差，也能有效地抑制硬閾值去噪方法中易產生的信號振蕩現(xiàn)象。但不懂它程序的編寫。    請教高手幫忙解決一下！謝謝了！function y = wthresh(x,sorh,t)%WTHRESH Perform soft or hard thresholding. %   Y = WTHRESH(X,SORH,T) returns soft (if SORH = 's')%   or hard (if SORH = 'h')

28、 T-thresholding  of the input %   vector or matrix X. T is the threshold value.%   Y = WTHRESH(X,'s',T) returns Y = SIGN(X).(|X|-T)+, soft %   thresholding is shrinkage.%   Y = WTHRESH(X,'h',T) returns Y = X.1_(|X|>T), ha

29、rd%   thresholding is cruder.%   See also WDEN, WDENCMP, WPDENCMP.%   M. Misiti, Y. Misiti, G. Oppenheim, J.M. Poggi 12-Mar-96.%   Last Revision: 24-Jul-2007.%   Copyright 1995-2007 The MathWorks, Inc.%        $Revision: .3 $switc

30、h sorh  case 's'    tmp = (abs(x)-t);    tmp = (tmp+abs(tmp)/2;    y   = sign(x).*tmp;  case 'h'    y   = x.*(abs(x)>t);  otherwise    error('Wavelet:FunctionArgVal:Invalid_ArgV

31、al', .        'Invalid argument value.')end這是傳統(tǒng)的閾值去噪法，只要你有新的算法，直接改就行了，但是我對該法有些疑惑，wthresh（x,sorh,thr）,若我用趙瑞珍法算出了每個尺度（假設J=3）的閾值，那么還用wthresh（x,sorh,thr）去噪時，thr是不是應該在之前寫成thr=thr1,thr2,thr3呢？求助高手：Matlab小波去噪中遇見的問題小波程序如下：function y=f1(x);X=xlsread('shoulian.xls

32、'); %讀入表格數(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;%畫出網格%用db4小波對原始信號進行3層分解并提取系數(shù)  還沒有寫好:要增加高頻系數(shù)圖c,l=wavedec(y,3,'db4'); %對y進行3層分解ca3=appcoef(c,l,'db4',3); %提取系數(shù)cd3=detcoef(c,l,3);cd2

33、=detcoef(c,l,2);cd1=detcoef(c,l,1);%對信號進行強制性去噪處理并圖示結果cdd3=zeros(1,length(cd3);cdd2=zeros(1,length(cd2);cdd1=zeros(1,length(cd1);c1=ca3' cdd3 cdd2 cdd1;s1=waverec(c1,l,'db4');subplot(2,2,2);plot(x,s1,'-m*');hold on;title('強制去噪后的信號');grid;%用默認閾值對信號進行去噪處理并圖示結果%用ddencmp()函數(shù)獲得

34、信號的默認閾值,使用wdencmp()命令函數(shù)實現(xiàn)去噪過程thr,sorh,keepapp=ddencmp('den','wv',y);s2=wdencmp('gbl',c,l,'db4',3,thr,sorh,keepapp);subplot(2,2,3);plot(x,s2,'-m*');hold on;title('默認閾值去噪后的信號');grid;%用給定的軟閾值進行去噪處理cdd1soft=wthresh(cd1,'s',1.465);cdd2soft=wthresh(c

35、d2,'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*');title('給定軟閾值去噪后的信號');gridclear;錯誤顯示：? Error using => horzcatAll matrices on a row in the bracketed expression must have th

36、e same number of rows.Error in => xnn at 41c2=ca3 cdd3soft cdd2soft cdd1soft;i以至于最后一個圖形沒法顯示同時處理后的圖形數(shù)據(jù)怎么調出來，以及第一個數(shù)據(jù)應該是【0，0】，但是處理后的不是。請問是我選擇的分解重構函數(shù)部隊還是其他原因請問如何修改 需要說明的是：收斂表格的兩組數(shù)據(jù)為：00.0020.4940.9451.0361.0371.48101.49111.86141.96191.78262.75282.78312.71342.54392.85472.71492.71545.43615.476

37、85.47755.20請教程序問題，心電信號去噪的matlab程序想把這個程序。改為只用軟硬閾值對比的心電信號去噪分析，該怎么辦？ 求高手。程序見下：%信號小波分解%基于Haar小波%ecg=fopen('100.dat','r');N=1201;data=fread(ecg,N,'int16');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('原始

38、心電信號');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',c,l,'Haar',3);a2=wrcoef('a',c,l,'Haar',2);a1=wrcoef('a',c,l,'Haar',1);subplot

39、(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(325);plot(a2);title('a2');axis(0 1201 -2 2);subplot(326);plot(a1);title('a1');axis(0 1201 -2 2);%基于db6小

40、波%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);title('原始心電信號');x=data;wname='db6'level=5;c,l=wavedec(x,level,wname);a5=wrcoef('a',c,l,'db

41、6',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);subplot(322);plot(a5);title('a5');axis(0 1201 -2000 2000);subplot(323);plot(a4);title('a4');axis(0 120

42、1 -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);%基于sym3小波%ecg=fopen('100.dat','r');N=1201;data=fread(ecg,N,'int16

43、9;);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('原始心電信號');x=data;wname='sym3'level=5;c,l=wavedec(x,level,wname);a5=wrcoef('a',c,l,'sym3',5);a4=wrcoef('a',c,l,'sym3',4);a3=wrc

44、oef('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 -2000 2000);subplot(323);plot(a4);title('a4');axis(0 1201 -2000 2000);subplot(324);plot(a3);title('a3');

45、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);%心電信號降噪%Birge-Massart策略閾值降噪%基于小波變換的心電信號的降噪ecg=fopen('100.dat','r');N=1201;data=fread(ecg,N,'int16');data=data/10000;fclo

46、se(ecg);x=data;wavename='db5'level=4;c,l=wavedec(x,level,wavename);alpha=1.5;sorh='h'thr,nkeep=wdcbm(c,l,alpha);%使用硬閾值給信號降噪xc,cxc,lxc,perf0,perfl2=wdencmp('lvd',c,l,wavename,level,thr,sorh);t1=0:0.004:(length(x)-1)*0.004;figure(4);subplot(211);plot(t1,x);title('從人體采集的原始的E

47、CG信號');subplot(212);plot(t1,xc);title('Birge-Massart策略閾值降噪后的ECG信號（wname=db5 level=4）');%最優(yōu)預測軟閾值降噪%基于小波變換的心電信號的壓縮ecg=fopen('100.dat','r');N=1201;data=fread(ecg,N,'int16');data=data/10000;fclose(ecg);x=data;wavename='db3'level=4;xd,cxd,lxd=wden(x,'heursu

48、re','s','mln',level,wavename);t=0:0.004:4.8;figure(5);subplot(311);plot(t,x);title('從人體采集的原始的ECG信號');subplot(312);plot(t,xd);title('heursure規(guī)則閾值降噪后的ECG信號(db3 level=4)');%最小均方差軟閾值降噪%基于小波變換的心電信號的壓縮ecg=fopen('100.dat','r');N=1201;data=fread(ecg,N,'

49、;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(311);plot(t,x);title('從人體采集的原始的ECG信號');subplot(313);plot(t,xd);title('Minimax規(guī)則閾值降噪后的ECG信號(wname=db3 level=

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

51、9;,c,l,wavename,4,thr,sorh);t=0:0.004:4.8;figure(7);subplot(511);plot(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,wavename);alpha=4;sorh='h'thr,nkeep=wdcbm(c,l,alpha);%使用硬閾值壓縮信號xc1,cxc

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

53、wavename,2,thr,sorh);subplot(514);plot(t,xc1);axis(0 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=wdencmp('lvd',c,l,'db3',2,thr,sorh);subplot(515);plot

54、(t,xc1);axis(0 5 -3 3);title('在第2層對高頻系數(shù)量化壓縮后的ECG信號alpha=4 db3');%求助解決信號降噪問題 求高手幫忙，我的一個實驗數(shù)據(jù)要進行信號降噪，需要用到Birge-Massar閾值降噪，penalty閾值降噪和缺省的閾值降噪這3種降噪方法，并通過比較得出最好的處理方法。如果有那位好心高手幫忙，小弟感激不盡啊！請幫忙的高手幫下忙??！按LZ的要求程序和得圖如下xx=load('實驗數(shù)據(jù).txt');x=xx(:,1)'y=xx(:,2)'x=x(end:-1:1);y=y(end:-1:1

55、);% 用sym6小波對信號做5層分解wname ='sym6' lev=5;c,l = wavedec(y,lev,wname);% 通過第1層的細節(jié)系數(shù)估算信號的噪聲強度 sigma = wnoisest(c,l,1);% 使用penalty策略確定降噪的閾值alpha=2;thr1 = wbmpen(c,l,sigma,alpha);% 使用Birge-Massart策略確定降噪的閾值thr2,nkeep = wdcbm(c,l,alpha);xd1 = wdencmp('gbl',c,l,wname,lev,thr1,'s',1

56、);% 用缺省的閾值確定的時候使用硬閾值時系數(shù)進行處理xd2,cxd,lxd,perf0,perf2= wdencmp('lvd',c,l,wname,lev,thr2,'h');% 求得缺省的閾值thr,sorh,keepapp = ddencmp('den','wv',x);% 重建降噪信號xd3 = wdencmp('gbl',c,l,wname,lev,thr,'s',1);figuresubplot(411); plot(y); title('原始信號','font

57、size',10);subplot(412); plot(xd1); title('使用penalty閾值降噪后信號','fontsize',10);subplot(413); plot(xd2); title('使用Birge-Massart閾值降噪后信號','fontsize',10);subplot(414); plot(xd2); title('使用缺省閾值降噪后信號', 'fontsize',10);eb1a.jpg (41.07 KB, 下載次數(shù): 4)小波去噪：為什么

58、閾值去噪后總把有用信息也給濾出了wavename='English.wav'%聲源layer=3;%分解尺度wavelets='db3'%小波選擇SNR_noise=10;%添加噪音的信噪比TPTR='rigrsure'%閾值選取原則sorh='s'%軟閾值 或 硬閾值scal='mln'%閾值每層調整%導入聲源 并 添加白噪聲s=wavread(wavename)'y=awgn(s(1,:),SNR_noise,'measured');%聲源為雙聲道 只提取一個聲道分析subplot(21

59、1),plot(s(1,:);title('Original Sound Wave');subplot(212),plot(y);title('Sound 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:layer    thre(n)=thselect(Dn,TPTR);end%對各系數(shù)進

60、行閾值去噪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)    subplot(layer,2,k+2),plot(D_recn);    title('Detail.cfs level af

61、ter threshing',num2str(n)    k=k+2;end%重構信號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);title('Sound Wave with White-Noise');subplot(313),plot(y_rec);title('

62、;De-Noised Sound Wave');為什么用thselect函數(shù)選取的閾值進行濾噪后 高頻系數(shù)全為零了 我只是想把高頻中的噪聲信號濾出，高手執(zhí)教 本帖最后由 cwjy 于 2010-5-2 07:52 編輯 1.jpg (117.12 KB, 下載次數(shù): 0)(TDC)的語音增強程序，適用于噪聲為白噪聲的情況 %基于Ephraim和Van Trees提出的信號子空間法(TDC)的語音增強程序，適用于噪聲為白噪聲的情況%        magic   

63、60;  SCUclear;%-參數(shù)定義-frame_len=40;              %幀長step_len=0.5*frame_len;    %分幀時的步長，相當于重疊50%N=8;     %計算Toeplitz協(xié)方差矩陣時用到的前后相鄰的幀數(shù),N為偶數(shù)%-讀入帶噪語音文件-filename,pathname=uigetfile('*.wav','請選擇語音文件：');wavin

64、,fs,nbits=wavread(pathname filename);wav_length=length(wavin);%取前3000采樣點用于估計噪聲方差n_var=var(wavin(1:3000),1);xv=n_var;      %定義Rx的特征值判別閾值%-分幀-inframe=(enframe(wavin,frame_len,step_len)'frame_num=size(inframe,2);%-臨時參數(shù)-%記錄dis_SNR以及max_Ax，查看下面的SNR變量時使用，在程序最后，對該變量進行制圖，根據(jù)情況使用max_Ax=zeros(1,frame_num);dis_SNR=zeros(1,frame_num);%-子空間法語音增強-%求噪聲的Toeplitz協(xié)方差矩陣，Rn=n_var*I;Rn=n_var*eye(frame_len);Ry