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

《小波变换语音消噪(改进阈值)》由会员分享，可在线阅读，更多相关《小波变换语音消噪(改进阈值)（23页珍藏版）》请在装配图网上搜索。

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); %对信号做2048点FFT变换f=fs*(0:1023)/2048;figure(1)plot(t,x1) %做原始语音信号的时域图形y=awgn(x1,10,measured);

2、%加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

3、=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)%求得信噪比 均方误差figure(3);plot(t,y1);function gd=Garrote_gg(a,b)%a为信号分解后的小波系数，b为获得的阈值m=0.2*(a*a)-(b*b)

4、;if (abs(a)=b) gd=sign(a)*(abs(a)-b/exp(m);else (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,y1);小波去噪软阈值和硬阈值的matlab仿真程序硬阈值、软阈值这里有一段不知道有用没%设置信噪比和随机种子值snr=4;init=2055615866;%产生原始信号sref和高斯白噪声污染的信号ssref,s=wn

5、oise(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=a3 ythard3 ythard2 ythard1;s3=waverec(c2,l,db1);%进行软阈值处理ytsoftd1=

6、wthresh(d1,s,thr);ytsoftd2=wthresh(d2,s,thr);ytsoftd3=wthresh(d3,s,thr);c3=a3 ytsoftd3 ytsoftd2 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小波除噪，为何硬阈值和

7、软阈值除躁信噪比一样了？matlab小波除噪，为何硬阈值和软阈值除躁信噪比一样了？load 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,1);%用db5小波对

8、原始信号进行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=a3 ythard3 ythard2 ythard1;f3=waverec(c2,l,db6);%进行软

9、阈值处理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(软阈值处理);%axis(1,1024,-1,1);snr=SNR_singlech(f1,f3)snr=SNR_singlech

10、(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); %noise powersnr=10*log10(Ps/Pn);小波去噪程序Matlab小波去噪（默认，强制，给定三种情况）% 利用小波分析对监测采集的信号进行去噪处理，恢复原始信

11、号%小波分析进行去噪有3中方法：%1、默认阈值去噪处理。该方法利用函数ddencmp( )生成信号的默认阈值，然后利用函数wdencmp( )进行去噪处理；%2、给定阈值去噪处理。在实际的去噪处理过程中，阈值往往可通过经验公式获得，且这种阈值比默认阈值的可信度高。在进行阈值量化处理时可利用函数wthresh( );%3、强制去噪处理。该方法是将小波分解结构中的高频系数全部置0，即滤掉所有高频部分，然后对信号进行小波重构。这种方法比较简单，且去噪后的信号比较平滑，但是容易丢失信号中的有用成分。% 利用小波分析对监测采集的水轮机信号进行去噪处理，恢复原始信号%Program Start% 载入监测

12、所得信号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,fontname,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);%高频部分

13、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)set(gca,fontname,times New Roman)set(gca,fontsize,14.0)%

14、 默认阈值对信号去噪并图示%用ddencmp( )函数获得信号的默认阈值，使用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)set(gca,fontsize,14.0)% 给定的软阈值进行去噪处理并图示cd1soft=wthresh(cd1,x,1.465);

15、%经验给出软阈值数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(给定软阈值去噪后信号);grid onset(gcf,color,w)set(gca,fontname,times New Roman)set(gca,fontsize,14.0)以上就是三种小波去噪的原程序。但红色标注的地方，我却运行不过去。提示分别问：1，?

16、 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);%经验给出软阈值数请问为什么？关于小波新阈值去噪改进方法里，matlab程序编写的讨论新阈值能有效地克服软阈值去噪方法中由于估计值与真实值之间的恒定偏差而带来的去噪误差，也能有效地抑

17、制硬阈值去噪方法中易产生的信号振荡现象。但不懂它程序的编写。 请教高手帮忙解决一下！谢谢了！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) T-thresholdingof the input% vector or matrix X. T is the threshold value.% Y = WTHRESH(X,s,T) returns Y = S

18、IGN(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 Revision: 24-Jul-2007.% Copyright 1995-2007 The MathWorks, Inc.% $Revision:

19、1.11.4.3 $switch sorhcase 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_ArgVal, . Invalid argument value.)end这是传统的阈值去噪法，只要你有新的算法，直接改就行了，但是我对该法有些疑惑，wthresh（x,sorh,thr）,若我用赵瑞珍法算出了每个尺度（假设J=3）的阈值，那么还用wthresh（x,sorh,thr）

20、去噪时，thr是不是应该在之前写成thr=thr1,thr2,thr3呢？求助高手：Matlab小波去噪中遇见的问题小波程序如下：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小波对原始信号进行3层分解并提取系数还没有写好:要增加高频系数图c,l=wavedec(y,3,db4); %对y进行3层分解ca3=

21、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=waverec(c1,l,db4);subplot(2,2,2);plot(x,s1,-m*);hold on;title(强制去噪后的信号);grid;%用默认阈值对信号进行去噪处理并图示结果%用d

22、dencmp()函数获得信号的默认阈值,使用wdencmp()命令函数实现去噪过程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(cd2,s,1.823);cdd3soft=wthresh(cd3,s,2.768);c2=ca3 cdd3soft

23、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 thesame number of rows.Error in = xnn at 41c2=ca3 cdd3soft cdd2soft cdd1soft;i以至于最后一个图形没法显示同时处理后的图形数据怎么调出来，以及第一个数据

24、应该是【0，0】，但是处理后的不是。请问是我选择的分解重构函数部队还是其他原因请问如何修改需要说明的是：收敛表格的两组数据为：00.0020.4940.9451.0361.0371.48101.49111.86141.96191.78262.75282.78312.71342.54392.85472.71492.71545.43615.47685.47755.20请教程序问题，心电信号去噪的matlab程序想把这个程序。改为只用软硬阈值对比的心电信号去噪分析，该怎么办？ 求高手。程序见下：%信号小波分解%基于Haar小波%ecg=fopen(100.dat,r);N=1201;data=fre

25、ad(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(原始心电信号);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=wrcoe

26、f(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(325);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=12

27、01;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,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,d

28、b6,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 -2000 2000);subplot(326);plot(a1);title(a1);axis(0 1201 -2000 2000);%基于sym3小波%ecg=

29、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=wavedec(x,level,wname);a5=wrcoef(a,c,l,sym3,5);a4=wrcoef(a,c,l,sym3,4);a3=wrcoef(a,c,l,sym3,

30、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);axis(0 1201 -2000 2000);subplot(325);plot(a2);title(a2);axis(0 1201 -2000 2000);subplot(326);plot(a1);tit

31、le(a1);axis(0 1201 -2000 2000);%心电信号降噪%Birge-Massart策略阈值降噪%基于小波变换的心电信号的降噪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(lv

32、d,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(ecg,N,int16);data=data/10000;fclose(ecg)

33、;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规则阈值降噪后的ECG信号(db3 level=4);%最小均方差软阈值降噪%基于小波变换的心电信号的压缩ecg=fopen(100.dat,r);N=1201;data=fread(ecg,N,int16);data=data/

34、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规则阈值降噪后的ECG信号(wname=db3 level=4);%心电信号压缩%基于小波变换的心电信号的压缩ecg=fopen(100.dat,r);N=1201;data=fread(ecg,N,int16);

35、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,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

36、 3);title(在第4层对高频系数量化压缩后的ECG信号（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);alpha=

37、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,l,alpha);%使用硬阈值压缩信号xc1,cxc1,lxc1,perf01,perfl21=wdencmp

38、(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种降噪方法，并通过比较得出最好的处理方法。如果有那位好心高手帮忙，小弟感激不尽啊！请帮忙的高手帮下忙啊！按LZ的要求程序和得图如下xx=load(实验数据.txt);x=xx(:,1);y=xx(:,2);x=x(end:-1:1);y=y(

39、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);% 使用Birge-Massart策略确定降噪的阈值thr2,nkeep = wdcbm(c,l,alpha);xd1 = wdencmp(gbl,c,l,wname,lev,thr1,s,1);% 用缺省的阈值确定的时候使用硬阈值时系数进行处理

40、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(原始信号,fontsize,10);subplot(412); plot(xd1); title(使用penalty阈值降噪后信号,fontsize,10);subplot(413); plot(xd2); title(使

41、用Birge-Massart阈值降噪后信号,fontsize,10);subplot(414); plot(xd2); title(使用缺省阈值降噪后信号, fontsize,10);eb1a.jpg(41.07 KB, 下载次数: 4)小波去噪：为什么阈值去噪后总把有用信息也给滤出了wavename=English.wav;%声源layer=3;%分解尺度wavelets=db3;%小波选择SNR_noise=10;%添加噪音的信噪比TPTR=rigrsure;%阈值选取原则sorh=s;%软阈值 或 硬阈值scal=mln;%阈值每层调整%导入声源 并 添加白噪声s=wavread(wav

42、ename);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,layer);%各层阈值选择for n=1:layer

43、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) subplot(layer,2,k+2),plot(D_recn); title(Detail.cfs level after threshing,num2str(n) k=k+2;end%重构信号C1=A;f

44、or 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(De-Noised Sound Wave);为什么用thselect函数选取的阈值进行滤噪后 高频系数全为零了 我只是想把高频中的噪声信号滤出，高手执教本帖最后由 cwjy 于 2010-5-2 07

45、:52 编辑1.jpg(117.12 KB, 下载次数: 0)(TDC)的语音增强程序，适用于噪声为白噪声的情况%基于Ephraim和Van Trees提出的信号子空间法(TDC)的语音增强程序，适用于噪声为白噪声的情况% 08.3.22 magic SCUclear;%-参数定义-frame_len=40; %帧长step_len=0.5*frame_len; %分帧时的步长，相当于重叠50%N=8; %计算Toeplitz协方差矩阵时用到的前后相邻的帧数,N为偶数%-读入带噪语音文件-filename,pathname=uigetfile(*.wav,请选择语音文件：);wavin,fs,

46、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_SNR以及max_Ax，查看下面的SNR变量时使用，在程序最后，对该变量进行制图，根据情况使用max_Ax=zeros(1,frame_num);dis_SNR=zeros

47、(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); %定义相邻6帧的自相关序列，长度等于帧长L=N*frame_len; %定义相邻N帧的长度outframe=zeros(frame_len,frame_num); %定义增强后的矩阵%对每一帧进行处理，得到增强后的信号,因为计算Toeplitz协方差矩阵要用到N帧，其中当前%帧处在第(N/2

48、+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帧的终点 seq(j+1)=

49、wavin(bgn_point:(end_point-j)*. wavin(bgn_point+j):end_point)/L; end; Ry=Toeplitz(seq);%2.对Ry进行特征分解，求它的特征向量矩阵Uy和特征值矩阵Ay% Uy,Ay=eig(Ry,nobalance); Uy,Ay=eig(Ry);%3.求出Ay中大于阈值xv的个数M，即为信号子空间维数，并由此求出语音子空间的特征值矩%阵Ax=Ay-Rn(维数:M*M),以及对应的特征向量矩阵Ux(维数：frame_len*M) I,J=find(Ayxv); M=length(I); Ax=zeros(M,M); A=A

50、y; Ux=zeros(frame_len,M); Ay_seq=zeros(1,frame_len); A=sort(Ay); seq1=A(frame_len,:); Ay_seq,IX=sort(seq1); for k=1:M num=frame_len-k+1; Ax(k,k)=Ay_seq(num)-xv; Ux(:,k)=Uy(:,IX(num); end;%4.求当前帧的增强后的信号 SNR=10*log10(sum(Ay_seq(frame_len-M+1:frame_len)/(frame_len*n_var);%SNR=-5时，u=5；-5SNR20时，u=u0-SNR/s；20=SNR时，u=1；其中u0=4.2,s=4.25 u0=4.2; s=6.25; if SNR=20 u=1; else u=u0-SNR/s; end;% u=3; %固定u值