搜索
%设置信噪比和随机种子值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('软阈值处理');%设置信噪比和随机种子值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小波消噪程序（原创）[Yo,FS,NBITS,OPTS]=wavread('normal_s3.wav');fs=FSY=Yo;dt=1/FS;YY=Yo(1800:9500)';N=length(YY);Y=YY+0.14*randn(1,N);time=(0:N-1)*dt;(N-1)*dtfY=fs*(1:N/2)/N;yYY=abs(fft(YY));yY=abs(fft(Y));[c,l]=wavedec(Y,4,'coif4');a4=appcoef(c,l,'coif4',4);a3=appcoef(c,l,'coif4',3);a2=appcoef(c,l,'coif4',2);a1=appcoef(c,l,'coif4',1);d4=detcoef(c,l,4);d3=detcoef(c,l,3);d2=detcoef(c,l,2);d1=detcoef(c,l,1);figure(1)subplot(511);plot(a4);title('a4');subplot(512);plot(d4);title('d4');subplot(513);plot(d3);title('d3');subplot(514);plot(d2);title('d2');subplot(5,1,5);plot(d1);title('d1');SNRin=10*log(norm(YY)^2/(norm(Y-YY))^2)%去噪前信号的信噪比figure(2)subplot(211);plot(time,YY);axis([0time(end)-1.51.5]);xlabel('time/s');title('原音信号');subplot(212);plot(fY,yYY(1:floor(N/2)));xlabel('f/hz');ylabel('频谱');title('原心音信号频谱图');xlim([0500])figure(3)subplot(211);plot(time,Y);axis([0time(end)-1.51.5]);xlabel('time/s');title('加噪心音信号');subplot(212);plot(fY,yY(1:floor(N/2)));xlabel('f/hz');ylabel('频谱');title('加噪心音信号频谱');xlim([0500])%%%%%%%%%%%单级阈值去噪%%%%%%%%%%%%%%%thr=yuzhi(Y,4,'coif4','sqtwolog','sln')dd1=yuzhifunction(d1,thr(1),'soft');dd2=yuzhifunction(d2,thr(2),'soft');dd3=yuzhifunction(d3,thr(3),'soft');dd4=yuzhifunction(d4,thr(4),'soft');cl=[a4,dd4,dd3,dd2,dd1];s=waverec(cl,l,'coif4');N1=length(s);times=(0:N1-1)*dt;fs0=fs*(0:N1/2-1)/N1;ys=abs(fft(s));%%%%%%%%%%%多级阈值去噪%%%%%%%%%%%%%%thr1=yuzhi(Y,4,'coif4','sqtwolog','mln')[s1,cxc,clc,perf2,perf3]=wdencmp('lvd',Y,'coif4',4,thr1,'s');N2=length(s1);times1=(0:N2-1)*dt;fs1=fs*(1:N2/2)/N2;ys1=abs(fft(s1));figure(4);subplot(321);plot(time,YY);axis([0time(end)-1.51.5]);xlabel('time/s');title('原音信号');subplot(322);plot(fY,yYY(1:floor(N/2)));xlabel('f/hz');ylabel('频谱');title('原心音信号频谱图');xlim([0400])subplot(323);plot(times,s);axis([0time(end)-1.51.5]);xlabel('time/s');title('单级阈值去噪心音信号');subplot(324);plot(fs0,ys(1:floor(N1/2)));xlabel('f/hz');ylabel('频谱');title('单级阈值去噪心音信号频谱图');xlim([0400])subplot(325);plot(times1,s1);axis([0time(end)-1.51.5]);xlabel('time/s');title('多级阈值去噪心音信号');subplot(326);plot(fs1,ys1(1:floor(N2/2)));xlabel('f/hz');ylabel('频谱');title('多级阈值去噪心音信号频谱图');xlim([0400])想把这个程序。改为只用软硬阈值对比的心电信号去噪分析，该怎么办？求高手。程序见下：%%%%%%%%%%信号小波分解%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%基于Haar小波%%%%%%%ecg=fopen('100.dat','r');N=1201;data=fread(ecg,N,'int16');data=data/10000;saveECGdatadata;fclose(ecg);x=0:0.004:4.8;figure(1);subplot(321);plot(x,data);axis([04.8-44]);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=wrcoef('a',c,l,'Haar',1);subplot(322);plot(a5);title('a5');axis([01201-22]);subplot(323);plot(a4);title('a4');axis([01201-22]);subplot(324);plot(a3);title('a3');axis([01201-22]);subplot(325);plot(a2);title('a2');axis([01201-22]);subplot(326);plot(a1);title('a1');axis([01201-22]);%%%%%%%%%%%%基于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([04.8-40004000]);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);subplot(322);plot(a5);title('a5');axis([01201-20002000]);subplot(323);plot(a4);title('a4');axis([01201-20002000]);subplot(324);plot(a3);title('a3');axis([01201-20002000]);subplot(325);plot(a2);title('a2');axis([01201-20002000]);subplot(326);plot(a1);title('a1');axis([01201-20002000]);%%%%%%%%%%%%%%%%%基于sym3小波%%%%%%%%%%%%%%%%%%%%%%%ecg=fopen('100.dat','r');N=1201;data