小波作业matlab程序

作业一
读入一个一维非平稳信号，对其进行3层小波分解与重构，提取指定尺度上的低频，高频系数。

保留10%，80%的模极大系数，进行重构，并做分析比较。

在一副图中显示出原始信号，分解后第2,3尺度上的高频，第3尺度上的低频，以及重构信号。

Matlab程序如下：
clear
clc
load noisdopp;
s=noisdopp;
[c,l]=wavedec(s,3,'db3');
ca1=appcoef(c,l,'db3',1);
ca2=appcoef(c,l,'db3',2);
ca3=appcoef(c,l,'db3',3);
%提取低频分量
cd1=detcoef(c,l,1);
cd2=detcoef(c,l,2);
cd3=detcoef(c,l,3);
%提取高频分量
coef1=0.1;
coef2=0.8;
%定义模值系数
[~,id1] = sort(abs(cd1),'descend');%对第一尺度高频信号排序并记录下其索引[m1,n1]=size(cd1);
number=round(n1*coef1);
nid1=id1(1:number);%取前10%的模值个数
ncd1=zeros(m1,n1);
for i=1:n1
if ismember(i,nid1)==1
ncd1(m1,i)=cd1(m1,i);
end
end
%第一尺度高频的10
[~,id2] = sort(abs(cd2),'descend');
[m2,n2]=size(cd2);
number=round(n2*coef1);
nid2=id2(1:number);
ncd2=zeros(m2,n2);
for i=1:n2
if ismember(i,nid2)==1
ncd2(m2,i)=cd2(m2,i);
end
end
%第二尺度高频的10%
[~,id3] = sort(abs(cd3),'descend'); [m3,n3]=size(cd3);
number=round(n3*coef1);
nid3=id3(1:number);
ncd3=zeros(m3,n3);
for i=1:n3
if ismember(i,nid3)==1
ncd3(m3,i)=cd3(m3,i);
end
end
%第三尺度高频的10%
[~,ia3] = sort(abs(ca3),'descend'); [ma3,na3]=size(ca3);
number=round(na3*coef1);
nia3=ia3(1:number);
nca3=zeros(ma3,na3);
for i=1:na3
if ismember(i,nia3)==1
nca3(ma3,i)=ca3(ma3,i);
end
end
%第三尺度低频的10%
c1=[nca3 ncd3 ncd2 ncd1];
x1=waverec(c1,l,'db3');%小波重构
[~,id11] = sort(abs(cd1),'descend'); [m11,n11]=size(cd1);
number2=round(n11*coef2);
nid11=id11(1:number2);
ncd11=zeros(m11,n11);
for i=1:n11
if ismember(i,nid11)==1
ncd11(m11,i)=cd1(m11,i);
end
end
%第一尺度高频的80%
[~,id21] = sort(abs(cd2),'descend'); [m21,n21]=size(cd2);
number2=round(n21*coef2);
nid21=id21(1:number2);
ncd21=zeros(m21,n21);
for i=1:n21
if ismember(i,nid21)==1
ncd21(m21,i)=cd2(m21,i);
end
end
%第二尺度高频的80%
[~,id31] = sort(abs(cd3),'descend'); [m31,n31]=size(cd3);
number2=round(n31*coef2);
nid31=id31(1:number2);
ncd31=zeros(m31,n31);
for i=1:n31
if ismember(i,nid31)==1
ncd31(m31,i)=cd3(m31,i);
end
end
%第三尺度高频的80%
[~,ia31] = sort(abs(ca3),'descend'); [ma31,na31]=size(ca3);
number=round(na31*coef2);
nia31=ia31(1:number);
nca31=zeros(ma31,na31);
for i=1:na31
if ismember(i,nia31)==1
nca31(ma31,i)=ca3(ma31,i);
end
end
%第三尺度低频的80%
cl2=[nca31 ncd31 ncd21 ncd11];
x2=waverec(cl2,l,'db3');%小波重构
subplot(421),plot(s);
subplot(422),plot(ca3);
subplot(423),plot(cd2);
subplot(424),plot(cd3);
subplot(425),plot(x1);
subplot(426),plot(x2);
结果如下图所示：
从图中可以看出，只保留10%模值极大系数，重构的信号与原信号相比偏差较大，保留80%模值极大系数，重构信号与原信号基本一致。

作业四
读入二维信号，加入高斯噪声，对其进行3层小波分解，并用小波阈值进行去噪。

选取两种不同的阈值：一是取Matlab中的阈值，二是自己给定阈值（如根据不同尺度上信号的方差，小波系数相关性等）。

计算均方误差，并比较分析两种不同去噪方法。

显示原始信号，第2,3尺度上的低频，以及硬，软阈值去噪后的图像。

Matlab程序如下：
clear
clc
load woman;
f=X;
x1=X+randn(size(X))*20;
[c,s]=wavedec2(x1,3,'db4');
a2=appcoef2(c,s,'db4',2);
a3=appcoef2(c,s,'db4',3); %提取低频分量
n=[1 2 3];%设置尺度向量
p=[80 23 10];%设置阈值向量
nc=wthcoef2('h',c,s,n,p,'s');
x0=waverec2(nc,s,'db4');%重构信号
[thr, sorh, keepapp] = ddencmp('den','wv',x1);
t = thselect(x1,'sqtwolog');
xd1 = wdencmp('gbl',x1,'db4',3,t,'h','1');%硬阈值去噪
xd2 = wdencmp('gbl',x1,'db4',3,t,'s','1');%软阈值去噪figure;
subplot(421),imshow(f,[]);title('原图')
subplot(422),imshow(x1,[]);title('加噪声后图')
subplot(423),imshow(a2,[]);title('第二尺度低频') subplot(424),imshow(a3,[]);title('第三尺度低频') subplot(425),imshow(xd1,[]);title('硬阈值去噪') subplot(426),imshow(xd2,[]);title('软阈值去噪') subplot(427),imshow(x0,[]);title('自设阈值去噪')
nf=norm(x1-f);%加入噪声的图均方误差
dnf1=norm(xd1-f);%硬阈值去噪后的均方误差
dnf2=norm(xd2-f);软阈值去噪后的均方误差
dnf3=norm(x0-f);%自设阈值去噪后的均方误差
加入噪声的图像与原图相比均方误差为：636.2014
硬阈值去噪后的均方误差为：633.8073
软阈值去噪后的均方误差为：558.2569
自设阈值去噪后的均方差为：589.4099
评价：自设的阈值去噪因为可以自己更改阈值向量的值，使得图像去噪效果向好的方向变换。

从主观视觉上来看，自设阈值的方法和软阈值去噪方法都比硬阈值的方法要好一些。

从客观评价指标，与原图的均方误差值来比较，自设阈值的方法比硬阈值方法要好，但不如matlab自带的软阈值方法。