小波变换语音消噪(改进阈值)

1、精选优质文档-倾情为你奉上改进阈值函数进行语音信号消噪，但是在程序运行过程中频频报错。本人经验不足调试不出，希望求得各位指导。改进函数表达式附图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层分解并提取系数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的无偏风险估计原理thr2=thselect(d2,'rigrsure');thr3=thselect(d3,'rigrsure');%利用改进阈值函数进行去噪处理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为信号分解后的小波系数，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)% 计算信噪比函数% I :原始信号% In:去噪后信号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, 下载次数: 0)改进函数表达式本帖最后由 罗志雄 于 2013-5-16 21:58 编辑function snr,mse=snrmse(I,In)% 计算信噪比函数% I :原始信号% In:去噪后信号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;                  %语音信号采样频率为8000xx=wavread('lw1.wav');x1=xx(:,1);%取单声道x1=x1-m

8、ean(x1);t=(0: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');   %加10d

9、b的高斯白噪声snr,mse=snrmsel(x1',y)     %求得信噪比 均方误差snr1=SNR_singlech(x1',y)figure(2)plot(t,y)                   %做加噪语音信号的时域图形c,l=wavedec(y,3,'db1');%多尺度一维分解%用db1小波对信号进行3层分解并提取系数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的无偏风险估计原理thr2=thselect(d2,'rigrsure');thr3=thselect(d3,'rigrsure');%利用改进阈值函数进行去噪处理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为信号分解后的小波系数，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) %求得信噪比 均方误差fig

12、ure(3);plot(t,y1);硬阈值、软阈值这里有一段不知道有用没%设置信噪比和随机种子值snr=4;init=;%产生原始信号sref和高斯白噪声污染的信号ssref,s=wnoise(1,11,snr,init);%用db1小波对原始信号进行3层分解并提取系数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(d1,'h',thr

13、);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 ytsoftd1;s4=waverec(c3

14、,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('软阈值处理'); load leleccum;index=1:1024; f1=leleccum(index); % 产生含噪信号 init=;randn(

15、'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,1);  %用db5小波对原始信号进行3层分解并提取系数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);%进行硬阈值处理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');%进行软阈值处理

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');%对上述信号进行图示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)信噪比函数SNR_singlech(I,In)function snr=SNR_singlech(I,In)% 计算信噪比函数% 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);% 利用小波分析对监测采集的信号进行去噪处理，恢复原始信号%小波分析进行去噪有3中方法：%1、默认阈值去噪处理。该方法利用函数ddencmp( )生成信号的默认阈值，然后利用函数wdencmp( )进行去噪处理；%2、给定阈值去噪处理。在实际的去噪处理过程中，阈值往往可通过经验公式获得，且这种阈值比默认阈值的可信度高。在进行阈值量化处理时可利用函数wthresh( );%3、强制去噪处理。该方法是将小波分解结构中的高频系数全部置0，即滤掉所有高频部分，然后对信号进行小波重构。这种方法比较简单，且去噪

20、后的信号比较平滑，但是容易丢失信号中的有用成分。% 利用小波分析对监测采集的水轮机信号进行去噪处理，恢复原始信号%Program Start% 载入监测所得信号load default.txt;  %装载采集的信号x= default;lx=length(x);t=0:1:length(x)-1' % 绘制监测所得信号subplot(2,2,1);plot(t,x);title('原始信号');grid onset(gcf,'color','w')  set(gca,'fontna

21、me','times New Roman')set(gca,'fontsize',14.0)% 用db1小波对原始信号进行3层分解并提取小波系数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);%高频部分% 对信号进行强制去噪处理并图示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('强制去噪后信号');grid onset(gcf,'color','w')  set(gca,'fontname','times New Roman')set(gca,'fontsize',14.0)% 默认阈值对信号去噪并图示%用ddencmp( )函数获

23、得信号的默认阈值，使用wdencmp( )函数实现去噪过程thr,sorh,keepapp=ddencmp('den','wv',x);x2=wdencmp('gbl',c,l,'db1',3,thr,sorh,keepapp);subplot(2,2,3);plot(x2);title('默认阈值去噪后信号');grid onset(gcf,'color','w')  set(gca,'fontname','times New Roman

24、')set(gca,'fontsize',14.0)% 给定的软阈值进行去噪处理并图示cd1soft=wthresh(cd1,'x',1.465);%经验给出软阈值数cd2soft=wthresh(cd2,'x',1.823); %经验给出软阈值数cd3soft=wthresh(cd3,'x',2.768); %经验给出软阈值数c2=ca3,cd3soft,cd2soft,cd1soft;x3=waverec(c2,1,'db1');subplot(2,2,4);plot(x3);title('给定

25、软阈值去噪后信号');grid onset(gcf,'color','w')  set(gca,'fontname','times New Roman')set(gca,'fontsize',14.0)以上就是三种小波去噪的原程序。但红色标注的地方，我却运行不过去。提示分别问：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);%经验给出软阈值数请问为什么？ 新阈值能有效地克服软阈值去噪方法中由于估计值与真实值之间的恒定偏差而带来的去噪误差，也能有效地抑制硬阈值去噪方法中易产生的信号振荡现象。但不懂它程序的编写。    请教高手帮忙解决一下！谢谢了！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这是传统的阈值去噪法，只要你有新的算法，直接改就行了，但是我对该法有些疑惑，wthresh（x,sorh,t

31、hr）,若我用赵瑞珍法算出了每个尺度（假设J=3）的阈值，那么还用wthresh（x,sorh,thr）去噪时，thr是不是应该在之前写成thr=thr1,thr2,thr3呢？小波程序如下：function y=f1(x);X=xlsread('shoulian.xls'); %读入表格数据x=X(:,1); %将表格数据的第一列赋值给x轴y=X(:,2); %将表格数据的第二列赋值给y轴subplot(2,2,1);plot(x,y,'-m*') %画出图形hold on;title('原始图形');grid;%画出网格%用db4小波对原始信

32、号进行3层分解并提取系数  还没有写好:要增加高频系数图c,l=wavedec(y,3,'db4'); %对y进行3层分解ca3=appcoef(c,l,'db4',3); %提取系数cd3=detcoef(c,l,3);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;s1=wa

33、verec(c1,l,'db4');subplot(2,2,2);plot(x,s1,'-m*');hold on;title('强制去噪后的信号');grid;%用默认阈值对信号进行去噪处理并图示结果%用ddencmp()函数获得信号的默认阈值,使用wdencmp()命令函数实现去噪过程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('默认阈值去噪后的信号');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');subplot(2,2,4);plot(x,s3,'-m*

35、');title('给定软阈值去噪后的信号');gridclear;错误显示：? 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以至于最后一个图形没法显示同时处理后的图形数据怎么调出来，以及第一个数据应该是【0，0】，但是处理后的不是。请问是我选择的分解重构函数部队还是其他原

36、因请问如何修改 需要说明的是：收敛表格的两组数据为：00.0020.4940.9451.0361.0371.48101.49111.86141.96191.78262.75282.78312.71342.54392.85472.71492.71545.43615.47685.47755.20想把这个程序。改为只用软硬阈值对比的心电信号去噪分析，该怎么办？ 求高手。程序见下：%信号小波分解%基于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('原始心电信号');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('原始心电信号');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('原始心电信号');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);%心电信号降噪%Birge-Massart策略阈值降

45、噪%基于小波变换的心电信号的降噪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);%使用硬阈值给信号降噪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信号');subplot(212);plot(t1,xc);title('Birge-Massart策略阈值降噪后的ECG信号（wname=db5 level=4）');%最优预测软阈值降噪%基于小波变换的心电信号的压缩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信号');subplot(312);plot(t,xd);title('heursure规则阈值降噪后的EC

48、G信号(db3 level=4)');%最小均方差软阈值降噪%基于小波变换的心电信号的压缩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信号');subplot(313);plot(t,xd);title('Minimax规则阈值降噪后的ECG信号(wname=db3 level=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(

50、x,4,wavename);alpha=2;sorh='h'thr,nkeep=wdcbm(c,l,alpha);%使用硬阈值压缩信号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信号');axis(0 5 -3 3);subplot(512);plot(t,xc);axis(0 5 -3 3);title('在第4层对高频系数量化压缩后的E

51、CG信号（alpha=2 db5）');c,l=wavedec(x,4,wavename);alpha=4;sorh='h'thr,nkeep=wdcbm(c,l,alpha);%使用硬阈值压缩信号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层对高频系数量化压缩后的ECG信号alpha=4 db5');c,l=wavedec(x,2,wavename

52、);alpha=4;sorh='h'thr,nkeep=wdcbm(c,l,alpha);%使用硬阈值压缩信号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层对高频系数量化压缩后的ECG信号alpha=4 db5');c,l=wavedec(x,2,'db3');alpha=4;sorh='h'thr,nkeep=wdcbm(c

53、,l,alpha);%使用硬阈值压缩信号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层对高频系数量化压缩后的ECG信号alpha=4 db3');% 求高手帮忙，我的一个实验数据要进行信号降噪，需要用到Birge-Massar阈值降噪，penalty阈值降噪和缺省的阈值降噪这3种降噪方法，并通过比较得出最好的处理方法。如果有那位好心高手帮忙，小弟感激不尽

54、啊！请帮忙的高手帮下忙啊！按LZ的要求程序和得图如下xx=load('实验数据.txt');x=xx(:,1)'y=xx(:,2)'x=x(end:-1:1);y=y(end:-1:1);% 用sym6小波对信号做5层分解wname ='sym6' lev=5;c,l = wavedec(y,lev,wname);% 通过第1层的细节系数估算信号的噪声强度 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);% 用缺省的阈值确定的时候使用硬阈值时系数进行处理xd2,cxd,lxd,perf0,perf2= wdencmp('lvd',c,l,wname,lev,thr2,'h');% 求得缺省的阈值thr,sorh,keepapp = ddencmp('den','wv',x);% 重建降噪信号xd3 = wdencm

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

57、4); plot(xd2); title('使用缺省阈值降噪后信号', 'fontsize',10); (41.07 KB, 下载次数: 4)wavename='English.wav'%声源layer=3;%分解尺度wavelets='db3'%小波选择SNR_noise=10;%添加噪音的信噪比TPTR='rigrsure'%阈值选取原则sorh='s'%软阈值 或 硬阈值scal='mln'%阈值每层调整%导入声源 并 添加白噪声s=wavread(wavename)

58、'y=awgn(s(1,:),SNR_noise,'measured');%声源为双声道 只提取一个声道分析subplot(211),plot(s(1,:);title('Original Sound Wave');subplot(212),plot(y);title('Sound Wave with White-Noise');%对含噪声源进行离散小波分解，并提取高低频系数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%对各系数进行阈值去噪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%重构信号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函数选取的阈值进行滤噪后 高频系数全为零了 我只是想把高频中的噪声信号滤出，高手执教 本帖最后由 cwjy 于 2010-5-2 07:52 编辑  (117.12 KB, 下载次数: 0) %基于Ephraim和Van Trees提出的信号子空间法(TDC)的语音增强程序，适用于噪声为白噪声的情况%   

62、08.3.22     magic      SCUclear;%-参数定义-frame_len=40;              %帧长step_len=0.5*frame_len;    %分帧时的步长，相当于重叠50%N=8;     %计算Toeplitz协方差矩阵时用到的前后相邻的帧数,N为偶数%-读入带噪语音文件-filename,pathnam

63、e=uigetfile('*.wav','请选择语音文件：');wavin,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);%-临时参数-%记录dis_

64、SNR以及max_Ax，查看下面的SNR变量时使用，在程序最后，对该变量进行制图，根据情况使用max_Ax=zeros(1,frame_num);dis_SNR=zeros(1,frame_num);%-子空间法语音增强-%求噪声的Toeplitz协方差矩阵，Rn=n_var*I;Rn=n_var*eye(frame_len);Ry=zeros(frame_len,frame_len);          %定义帧信号的Toeplitz协方差矩阵seq=zeros(1,frame_len);     &#

65、160;           %定义相邻6帧的自相关序列，长度等于帧长L=N*frame_len;                          %定义相邻N帧的长度outframe=zeros(frame_len,frame_num);    %定义增强后的矩阵%对每一帧进行处理，得到增强后的信号,因为计算Toe

66、plitz协方差矩阵要用到N帧，其中当前%帧处在第(N/2+1)帧（注意：计算Toeplitz协方差矩阵的帧是不重叠的，而inframe相邻帧%之间重叠50%），因此，inframe的前N帧和最后(N/2-1)*2帧无法估计Toeplitz协方差矩阵，%且这(N-1)*2帧一般是噪声信号，故不必对其处理，直接将out_frame的这些帧置0for i=(N+1):(frame_num-N-2)%1.构造当前帧的Toeplitz协方差矩阵    for j=0:(frame_len-1)        bgn_point=(i-N-1)*step_len+1;       %相邻6帧的起始点        end_point=(i+N-1)*step_len;         %相邻6帧的终点       &