基于Matlab的小波分解去噪与重构.docx
《基于Matlab的小波分解去噪与重构.docx》由会员分享，可在线阅读，更多相关《基于Matlab的小波分解去噪与重构.docx（20页珍藏版）》请在冰豆网上搜索。
基于Matlab的小波分解去噪与重构
《现代信号处理》
大作业
基于Matlab的小波分解、去噪与重构
1作业内容及要求3...
1.1作业内容3...
1.2作业要求3...
2系统原理3...
2.1小波变换原理3...
2.2阈值去噪原理3...
3系统分析及设计5...
3.1图像分解5...
3.2高频去噪5...
3.3图像重构6...
4程序编写7...
4.1main函数7...
4.2分解函数9...
4.2.1二维分解函数9..
4.2.2一维分解函数1..0
4.3卷积函数1..0.
4.4采样函数1..1.
4.4.1下采样函数1..1
4.4.2上采样函数1..1
4.5重构函数1..2.
4.5.1二维重构函数1..2
4.5.2一维重构函数1..3
5结果分析及检验1..4.
5.1结果分析1..4.
5.2结果检验1..6.
6心得体会1..8.
参考文献1..9..
作业内容及要求
1.1作业内容
用小波对图像进行滤波分解、去噪，然后重构。
1.2作业要求
用小波对图像进行滤波分解、去噪，然后重构。
具体要求：
（1）被处理图像可选择：
woman,wbarb,wgatlin,detfingr,tire.；
（2）可以选择db等正交小波、或双正交小波（或用几种小波）；
（3）用选用小波的分解滤波器通过定义的卷积函数conv_my（）对图像二维数组进行小波分解，并进行下采样，获取CA、CV、CD、CH等分解子图；
（4）对高频信号子图进行去噪处理，可以采用软阈值、硬阈值等方法；
（5）用选用小波的综合滤波器对去噪的子图进行图像重构。
系统原理
2.1小波变换原理
小波变换的一级分解过程是，先将信号与低通滤波器卷积再下采样可以得到低频部分的小波分解系数再将信号与高通滤波器卷积后下采样得到高频部分的小波分解系数；而多级分解则是对上一级分解得到的低频系数再进行小波分解，是一个递归过程。
二维小波分解重构可以用一系列的一维小波分解重构来实现。
重构则是分解的逆过程，对低频系数、高频系数分别进行上抽样和低通、高通滤波处理。
要注意重构时同一级的低频、高频系数的个数必须相等。
2.2阈值去噪原理
图像去噪的方法是：
（1）图像的小波分解。
选择合适的小波函数以及合适的分解层数对图像进行分解。
（2）对分解后的高频系数进行阈值处理。
对分解的每一层，选择合适的阈值对该层的水平、垂直和对角三个方向的高频系数进行阈值处理。
（3）重构图像。
根据小波分解的低频系数和经阈值量化处理后的高频系数进行图像重构。
本设计采用软阈值去噪，其原理为：
当小波系数的绝对值小于给定的阈值时，
令其为0，大于阈值时，令其都减去阈值，即：
[sign（w）]（|w|）,|w|
w{
0,|w|
小波阈值入在去噪过程中起到决定性的作用。
如果入太小，那么施加阈值以后的小波系数中将包含过多的噪声分量，达不到去噪的目的；反之，如果入太大，那么将去除一部分信号的分量，从而使由小波系数重建后的信号产生过大的失真。
MATLAB中实现阈值获取的函数有ddenmp、thselect、wbmpen和wdcbm，本次设计中采用ddenmp函数进行阈值获取。
调用格式为[THR,SORH,KEEPAPP]=ddencmp（'den','wv',X），函数ddencmp用于获取信号在消噪或压缩过程中的默认阈值。
输入参数X为一维或二维信号；'den表示进行去噪；'wv'表示选择小波。
返回值THR是返回的阈值；SORH是软阈值或硬阈值选择参数；KEEPAPP表示保存低频信号。
三系统分析及设计
3.1图像分解
图像分解在程序中分为两部分：
一维分解以及二维分解。
（1）一维分解先将信号与低通滤波器卷积再下采样可以得到低频部分的小波分解系数，再将信号与高通滤波器卷积后下采样得到高频部分的小波分解系数
（2）对输入的信号，也就是前面load图像之后得到的X矩阵，先对每一行
进行一维分解，在对分解后得到的矩阵的每一列进行一维分解，最后得到一个矩
QQQQQQQQ
阵X2=[?
?
?
?
?
?
?
?
可以将CA、CH、CV、CD依次输出得到分解后的低频、垂直、水平、对角分量。
原始信号分别进行低通、高通滤波，再分别对列进行二元下抽样，就得到低频、高频（也称为平均、细节）两部分系数；再对这两部分系数进行小波分解，步骤与第一级分解一致，最后得到cA:
低频分量，cH:
水平分量，cV:
垂直分量，cD:
对角线分量。
流程图如下图1所示。
列卷积
图1分解流程图
3.2咼频去噪
本次设计采用软阈值去噪方法，将对分解出的cH、cV、cD进行阈值去噪，
选用ddencmp函数获取去噪过程中的默认软阈值。
其调用格式为：
[THR,SORH,KEEPAPP]=ddencmp（IN1,IN2,X）。
IN仁'den',IN2=W，X为一维或二维矩阵信号。
wthcoef2进行
在获取默认阈值后，对该阈值进行软阈值处理，调用到函数
处理，其调用过程中关键部分是
tmp=（abs（xH）；
tmp=（tmp+abs（tmp））/2;
x代表小波系数的大
y=sign（x）.*tmp;其中t代表前面获得的默认阈值,小。
得到的y即为软阈值处理后的小波系数。
3.3图像重构
重构的过程与分解过程类似，同样分为两部分：
一维重构和二维重构。
（1）一维重构：
先对平均部分系数进行上采样，再与低通滤波器卷积。
然后
对细节部分选取重构所需的细节部分，长度与本层平均部分系数相同，再对其上采样后进行高通滤波器卷积。
将两个得到的结果相加得到新的平均部分系数。
重
复以上操作知道细节部分长度小于平均部分长度。
（2）二维重构：
根据重构流程图可知，先对cA和cH进行重构，再对cV和cD重构，在该程序中所以在该程序中
X2=[?
?
?
?
?
?
?
?
所以先对行左右两半部分进行重构，重构后得到的矩阵在进行列上下两半部分进行重构，最终得到的重构图像。
具体步骤是先将去噪后的子图分别对行进行，再将低频分量和垂直分量进行二元上采样、低通滤波并合成；将水平分量和对角线分量进行二元上采样、高通滤波并合成，就得到低频、高频（也称为平均、细节）两部分系数；再对这两部分系数进行小波合成，步骤与第一级重构一致，最后通过wkeep选取我们所需
的原始图像的大小，得到原始图像。
该流程图如下图2所示。
列卷积
图2重构流程图
X'
四程序编写
4.1main函数
clearall;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%db4小波分解成低通和高通滤波器
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%figure;
'db4'）;
on;
系数'）;ylabel（'Lo_D'）;
on;
系数'）;ylabel（'Hi_D'）;
on;
系数'）;ylabel（'Lo_R'）;
on;
系数'）;ylabel（'Hi_R'）;
[Lo_D,Hi_D,Lo_R,Hi_R]=wfilters（subplot（221）;stem（Lo_D）;gridtitle（'分解低通滤波器'）;xlabel（subplot（222）;stem（Hi_D）;grid
title（'分解高通滤波器'）;xlabel（subplot（223）;stem（Lo_R）;grid
title（'综合低通滤波器'）;xlabel（subplot（224）;stem（Hi_R）;grid
title（'综合高通滤波器'）;xlabel（%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%导入和显示原始图像
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%loadwoman;
photo=X;%X含有被装载的信号
figure;
colormap（map）;%设置颜色映射到矩阵图
image（photo）;
axisoff;%关闭原始信号的坐标轴
title（'原始信号'）;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%用自定义的一维和二维函数分解图像
%得到低频部分（cA）,水平部分（cH）,垂直部分（cV）,对角线部分（cD）子图
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%[cA,cH,cV,cD,x1]=mydwt2（photo,Lo_D,Hi_D,map）;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%显示经过行分解和下采样之后的图像
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[n3,n4]=size（x1）;
L0=x1（1:
n3,1:
n4/2）;H0=x1（1:
n3,n4/2+1:
n4）;
figure;colormap（map）;image（L0）;title（'行分解后的低频分量'）;
figure;colormap（map）;image（H0）;title（'行分解后的高频分量'）;
figure;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%对高频部分（cH,cV,cD）进行去噪处理
%得到新的高频子图（new_cH,new_cV,new_cD）
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%用小波函数sym5对cH进行1层小波分解
[C1,S1]=wavedec2（cH,1,'sym5'）;%二维小波多层分解，这里取1，即为1层，输出
C1为各层分解系数，S1为各层分解系数长度，下同
L=size（x1）;
%求取信号去噪的默认阈值、软阈值、并保留低频部分
[thr1,sorh1,keepapp1]=ddencmp（'den','wv',cH）;
%对三个方向高频系数进行软阈值处理%wthcoef2的调用格式为NC=wthcoef2（'type',C,S,N,T,SORH）
%其中'type'分别取'h','v','d'代表三个方向
%N为尺度向量,
T为阈值向量，SORH取'h'和's'分别表示使用硬阈值或软阈值，下同
new_cH=waverec2（nc1,S1,'sym5'）;%对新的小波分解结构（nc，s）进行重构
new_cHH=wkeep（new_cH丄，'c'）;%提取重构的图像，'c'代表中心部分，但因为L为
nc1=wthcoef2（
'h',C1,S1,1,thr1,
nc1=wthcoef2（
'v',C1,S1,1,thr1,
nc1=wthcoef2（
'd',C1,S1,1,thr1,
's'）;
's'）;
's'）;
完整长度，所以用T或'r'也可以，下同
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%用小波函数sym5对cV进行1层小波分解
[C2，S2]=wavedec2（cV，1，'sym5'）;
%求取信号去噪的默认阈值、软阈值、并保留低频部分
[thr2，sorh2，keepapp2]=ddencmp（
'den'
，'wv'，cV）;
%对三个方向高频系数进行软阈值处理
nc2=wthcoef2（
'h'
，C2，S2，1，thr2，
's'）;
nc2=wthcoef2（
'v'
，C2，S2，1，thr2，
's'）;
nc2=wthcoef2（
'd'
，C2，S2，1，thr2，
's'）;
new_cV=waverec2（nc2，S2，'sym5'）;%对新的小波分解结构（nc，s）进行重构
new_cVV=wkeep（new_cV，L，'c'）;%提取重构的图像
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%用小波函数sym5对cD进行1层小波分解
[C3，S3]=wavedec2（cD，1，'sym5'）;
%求取信号去噪的默认阈值、软阈值、并保留低频部分
[thr3，sorh3，keepapp3]=ddencmp（'den'，'wv'，cD）;
%对三个方向高频系数进行软阈值处理
nc3=wthcoef2（'h'，C3，S3，1，thr3，'s'）;
nc3=wthcoef2（'v'，C3，S3，1，thr3，'s'）;
nc3=wthcoef2（'d',C3,S3,1,thr3,'s'）;new_cD=waverec2（nc3,S3,'sym5'）;%对新的小波分解结构（nc，s）进行重构new_cDD=wkeep（new_cD,L,'c'）;%提取重构的图像%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%subplot（221）;colormap（map）;image（new_cHH）;title（'去噪后的水平分量'）;
subplot（222）;colormap（map）;image（new_cVV）;title（'去噪后的垂直分量'）;
subplot（223）;colormap（map）;image（new_cDD）;title（'去噪后的对角线分量
'）;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%cA和去噪后的子图（new_cH,new_cV,new_cD）进行重构
%reconstruction为重构之后的图形%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
reconstruction=myidwt2（cA,new_cH,new_cV,new_cD,Lo_R,Hi_R,map）;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4.2分解函数
4.2.1二维分解函数
function[LL,HL,LH,HH,x1]=mydwt2（x,Lo_D,Hi_D,map）
%函数MYDWT2（）对输入的r*c维矩阵x进行二维小波分解，输出四个分解系数子矩阵[LL,HL,LH,HH]
%输入参数：
x——输入矩阵，为r*c维矩阵。
%输出参数：
LL,HL,LH,HH——是分解系数矩阵的四个相等大小的子矩阵,大小均为r/2*c/2维
%LL：
低频部分分解系数；HL：
垂直方向分解系数；
%LH：
水平方向分解系数；HH：
对角线方向分解系数。
%首先对输入矩阵的每一行序列进行
%将分解系数序列再存入矩阵x1中，得
%再对输入矩阵的每一列序列进行一维
[row,col]=size（x）;%读取输入矩阵的大小forj=1:
row
维离散小波分解
tmp1=x（j,:
）;
[ca1,cd1,n1,n2]=mydwt（tmp1,Lo_D,Hi_D,1）;
a1=n1+n2;
x1（j,:
）=[ca1,cd1];
到[L|H]
end
L=size（x1）;%取x1的长度
fork=1:
a1离散小波分解tmp2=x1（:
k）;
[ca2,cd2,n1,n2]=mydwt（tmp2,Lo_D,Hi_D,1）;
x2（:
k）=[ca2,cd2];
[LL,HL;LH,HH]
end
[row,col]=size（x2）;
LL=x2（1:
row/2,1:
col/2）;
LH=x2（row/2+1:
row,1:
col/2）;
HL=x2（1:
row/2,col/2+1:
col）;
HH=x2（row/2+1:
row,col/2+1:
col）;
figure;
subplot（221）;colormap（map）;image（LL）;title（subplot（222）;colormap（map）;image（LH）;title（subplot（223）;colormap（map）;image（HL）;title（subplot（224）;colormap（map）;image（HH）;title（
%将分解所得系数存入矩阵x2中，得到
%LL是矩阵x2的左上角部分
%LH是矩阵x2的左下角部分
%HL是矩阵x2的右上角部分
%HH是矩阵x2的右下角部分
'分解后的低频分量'）;'分解后的垂直分量'）;'分解后的水平分量'）;
'分解后的对角分量'）;
4.2.2一维分解函数
对输入序列x进行一维离散小波分
解，
输出分解序列
[cA,cD]
%输入参数：
x—
—输入序列；
%
lpd
——低通滤波器；
%
hpd
——高通滤波器；
%
dim
——小波分解级数。
function[cA,cD,n1,n2]=mydwt（x,Lo_D,Hi_D,dim）
%函数[cA,cD,n1,n2]=MYDWT（X,LPD,HPD,DIM）
%输出参数：
cA——低频部分的小波分解系数；%cD——高频部分的小波分解系数。
cA=x;%初始化cA，cD
cD=[];
fori=1:
dim
cvl=myconv1（cA,Lo_D）;dnl=downspl（cvl）;n1=length（dnl）;
cvh=myconv1（cA,Hi_D）;dnh=downspl（cvh）;n2=length（dnh）;cA=dnl;
%低通滤波
%通过下采样求出平均部分的分解系数
cD=[cD,dnh];
end
%高通滤波
%通过下采样求出本层分解后的细节部分系数
%下采样后的平均部分系数进入下一层分解
%将本层分解所得的细节部分系数存入序列cD
4.3卷积函数
functionc=myconv1（a,b）;
m=length（a）;%获取a矩阵的长度
n=length（b）;%获取b矩阵的长度
fori=1:
m+n-1
c（i）=0;%定义长度为m+n-1的全为0的矩阵c
end
fori=1:
m
forj=1:
n
x=a（i）*b（j）;%将a矩阵中的数乘以b矩阵中的数
c（i+j-1）=c（i+j-1）+x;%将上一行得到的数累加到对应序号的c矩阵数中
end
end
4.4采样函数
4.4.1下采样函数
functiony=downspl（x）
%函数Y=DOWMSPL（X）对输入序列进行下抽样，输出序列Y。
%下抽样是对输入序列取其偶数位，舍弃奇数位。
例如x=[x1,x2,x3,x4,x5]，则
y=[x2,x4].
带小数时取整数部分）
N=length（x）;%读取输入序列长度
M=floor（N/2）;%输出序列的长度是输入序列长度的一半
i=1:
M;
y（i）=x（2*i）;
4.4.2上采样函数
functiony=upspl（x）
%函数Y=UPSPL（X）对输入的一维序列x进行上抽样，即对序列x每个元素之间
%插零，例如x=[x1,x2,x3,x4],上抽样后为y=[x1,0,x2,0,x3,0,x4];
N=length（x）;%读取输入序列长度
M=2*N-1;
fori=1:
M
%输出序列的长度是输入序列长度的2倍再减一
%输出序列的偶数位为0，奇数位按次序等于相应位置的输入序列
ifmod（i,2）y（i）=x（（i+1）/2）;else
y（i）=0;endend
4.5重构函数
4.5.1二维重构函数
functiony=myidwt2（LL,HL,LH,HH,Lo_R,Hi_R,map）%函数MYIDWT2（）对输入的子矩阵序列进行逆小波变换，重构出矩阵y
%输入参数：
LL,HL,LH,HH——是四个大小均为r*c维的矩阵
%输出参数：
y——是一个大小为2r*2c维的矩阵
%将输入的四个矩阵组合为一个矩阵
%求出组合矩阵的行列数
%首先对组合矩阵tmp_mat的每一列，分
%分开的两部分分别作为平均系数序列ca1、
%重构序列
%将重构序列存入待输出矩阵yt的
%将输出矩阵y的每一行，分开成左右两半
%分开的两部分分别作为平均系
tmp_mat=[LL,HL;LH,HH];
flag=0;%循环标志位
fori=1:
2
ifflag==1[row,col]=size（yt）;
fork=1:
col
开成上下两半ca1=yt（1:
row/2,k）;
细节系数序列cd1cd1=yt（row/2+1:
row,k）;tmp1=myidwt（ca1,cd1,Lo_R,Hi_R）;yt1（:
k）=tmp1;
相应列，此时y=[L|H]
end
break;
end
ifflag==0
[row,col]=size（tmp_mat）;
forj=1:
rowca2=tmp_mat（j,1:
col/2）;
数序列ca2、细节系数序列cd2cd2=tmp_mat（j,col/2+1:
col）;tmp2=myidwt（ca2,cd2,Lo_R,Hi_R）;yt（j,:
）=tmp2;
应行，得到最终的输出矩阵y=yt
%重构序列
%将重构序列存入待输出矩阵yt的相
end
flag=1
升级会员 