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

..改进阈值函数进行语音信号消噪，但是在程序运行过程中频频报错。本人经验不足调试不出，希望求得各位指导。改进函数表达式附图clearall;clc;closeall;fs=8000;%语音信号采样频率为8000xx=wavread('lw1.wav');x1=xx(:,1);%取单声道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');%加10db的高斯白噪声[snr,mse]=snrmse(x1,y')%求得信噪比均方误差figure(2)plot(t,y)%做加噪语音信号的时域图形[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);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=[a3gd3gd2gd1];y1=waverec(c2,l,'db1');%多尺度重构[snr,mse]=snrmse(x1,y1')%求得信噪比均方误差figure(3);plot(t,y1);functiongd=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;end..function[snr,mse]=snrmse(I,In)%计算信噪比函数%I:原始信号%In:去噪后信号snr=0;Ps=sum(sum((I-mean(mean(I))).^2));%signalpowerPn=sum(sum((I-In).^2));%noisepowersnr=10*log10(Ps/Pn);mse=Pn/length(I);QQ截图20130516175535.png(11.18KB,下载次数:0)改进函数表达式本帖最后由罗志雄于2013-5-1621:58编辑function[snr,mse]=snrmse(I,In)%计算信噪比函数%I:原始信号%In:去噪后信号snr=0;Ps=sum(sum((I-mean(mean(I))).^2));%signalpowerPn=sum(sum((I-In).^2));%noisepowersnr=10*log10(Ps/Pn);mse=Pn/length(I);修改后程序清单如下：clearall;clc;closeall;fs=8000;%语音信号采样频率为8000xx=wavread('lw1.wav');x1=xx(:,1);%取单声道x1=x1-mean(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');%加10db的高斯白噪声[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);%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=Garrote_gg(d2,thr2);gd3=Garrote_gg(d3,thr3);c1=[a3gd3gd2gd1];functiongd=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));elsegd=zeros(size(a));endy1=waverec(c1,l,'db1');%多尺度重构[snr,mse]=snrmsel(x1',y1)%求得信噪比均方误差figure(3);plot(t,y1);..小波去噪软阈值和硬阈值的matlab仿真程序硬阈值、软阈值这里有一段不知道有用没%设置信噪比和随机种子值snr=4;init=2055615866;%产生原始信号sref和高斯白噪声污染的信号s[sref,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);ythard2=wthresh(d2,'h',thr);ythard3=wthresh(d3,'h',thr);c2=[a3ythard3ythard2ythard1];s3=waverec(c2,l,'db1');%进行软阈值处理ytsoftd1=wthresh(d1,'s',thr);ytsoftd2=wthresh(d2,'s',thr);ytsoftd3=wthresh(d3,'s',thr);c3=[a3ytsoftd3ytsoftd2ytsoftd1];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小波除噪，为何硬阈值和软阈值除躁信噪比一样了？loadleleccum;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,1]);%用db5小波对原始信号进行3层分解并提取系数[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,'h',thr);c2=[a3ythard3ythard2ythard1];f3=waverec(c2,l,'db6');%进行软阈值处理ytsoftd1=wthresh(d1,'s',thr);ytsoftd2=wthresh(d2,'s',thr);ytsoftd3=wthresh(d3,'s',thr);c3=[a3ytsoftd3ytsoftd2ytsoftd1];f4=waverec(c3,l,'db6');%对上述信号进行图示subplot(2,2,3);plot(f3);title('硬阈值处理');%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)信噪比函数SNR_singlech(I,In)functionsnr=SNR_singlech(I,In)%计算信噪比函数%I：riginalsignal%In:noisysignal(ie.originalsignal+noisesignal)..snr=0;Ps=sum(sum((I-mean(mean(I))).^2));%signalpowerPn=sum(sum((I-In).^2));%noisepowersnr=10*log10(Ps/Pn);小波去噪程序Matlab小波去噪（默认，强制，给定三种情况）%%利用小波分析对监测采集的信号进行去噪处理，恢复原始信号%小波分析进行去噪有3中方法：%1、默认阈值去噪处理。该方法利用函数ddencmp()生成信号的默认阈值，然后利用函数wdencmp()进行去噪处理；%2、给定阈值去噪处理。在实际的去噪处理过程中，阈值往往可通过经验公式获得，且这种阈值比默认阈值的可信度高。在进行阈值量化处理时可利用函数wthresh();%3、强制去噪处理。该方法是将小波分解结构中的高频系数全部置0，即滤掉所有高频部分，然后对信号进行小波重构。这种方法比较简单，且去噪后的信号比较平滑，但是容易丢失信号中的有用成分。%%利用小波分析对监测采集的水轮机信号进行去噪处理，恢复原始信号%ProgramStart%%载入监测所得信号loaddefault.txt;%装载采集的信号x=default;lx=length(x);t=[0:1:length(x)-1]';%%绘制监测所得信号subplot(2,2,1);plot(t,x);title('原始信号');gridonset(gcf,'color','w')set(gca,'fontname','timesNewRoman')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(