|
楼主 |
发表于 2012-4-6 20:52
|
显示全部楼层
谢老师您好!我在进一步处理的时候,又碰到问题了。代码如下
本帖最后由 mingmingtree 于 2012-4-6 20:54 编辑
回复 11 # xiezhh 的帖子
- clear all;
- clc;
- a = 5.8;
- b1 = 1528;b2=1535;b3=1560;b4=1583;
- c1 = 0.3;c2=0.2;c3=0.25;c4=0.3;
- g1 = @(t)a*exp(-4*log(2)*(((t-b1)/c1).^2));
- g2 = @(t)a*exp(-4*log(2)*(((t-b2)/c2).^2));
- g3 = @(t)a*exp(-4*log(2)*(((t-b3)/c3).^2));
- g4 = @(t)a*exp(-4*log(2)*(((t-b4)/c4).^2));
- x = 1510:0.005:1592;
- y1 = g1(x);
- y2 = g2(x);
- y3 = g3(x);
- y4 = g4(x);
- s1=y1+y2+y3+y4;
- size(s1);%纯信号
- z=awgn(s1,-5,'measured');%加入噪声
- s=z(1:16002);%去前16002个点作为程序的原始含噪信号
- ls=length(s);
- subplot(2,2,1) ,plot(x,z);title('含噪声的信号');grid;
- %用db1小波对原始信号进行3层分解并提取系数
- [c,l]=wavedec(s,3,'db1');
- ca3=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(cd2));
- cdd1=zeros(1,length(cd1));
- c1=[ca3 cdd3 cdd2 cdd1];
- s1=waverec(c1,l,'db1');
- subplot(2,2,2);
- plot(s1);
- title('强制消噪后的信号');grid;
- %用默认阈值对信号进行消噪处理并图示结果
- %用ddencmp函数获得信号的默认阈值,使用wdencmp命令函数实现消噪过程
- [thr,sorh,keepapp]=ddencmp('den','wv',s);
- s2=wdencmp('gbl',c,l,'db1',3,thr,sorh,keepapp);
- subplot(2,2,3);
- plot(s2);
- title('默认阈值消噪后的信号');grid;
- %用给定的软阈值进行消噪处理
- cd1soft=wthresh(cd1,'s',1.465);
- cd2soft=wthresh(cd2,'s',1.823);
- cd3soft=wthresh(cd3,'s',2.768);
- c2=[ca3 cd3soft cd2soft cd1soft];
- s3=waverec(c2,l,'db1');
- subplot(2,2,4);
- plot(s3);
- title('给定软阈值消噪后的信号');grid
- %function y=snr(s,s1);%s是原始信号,s1,s2,s3分别是三种阈值去噪后的信号
- y1=sum(s1.^2); y2=sum(s2.^2);y3=sum(s3.^2);
- y12=sum((s-s1).^2);y22=sum((s-s2).^2);y32=sum((s-s3).^2);%y2=sum((abs(x1)-abs(x2)).^2)
- y13=10*log10((y1/y12))
- y23=10*log10((y2/y22))
- y33=10*log10((y3/y32))
复制代码
|
|