基于Matlab的小波分解去噪与重构

1、大作业基于Matlab的小波分解、去噪与重构目录一 作业内容及要求31.1 作业内容31.2 作业要求3二 系统原理32.1 小波变换原理32.2 阈值去噪原理3三 系统分析及设计53.1 图像分解53.2 高频去噪53.3 图像重构6四 程序编写74.1 main函数74.2 分解函数94.2.1 二维分解函数94.2.2 一维分解函数104.3 卷积函数104.4 采样函数114.4.1 下采样函数114.4.2 上采样函数114.5 重构函数124.5.1 二维重构函数124.5.2 一维重构函数13五 结果分析及检验145.1 结果分析145.2 结果检验16六 心得体会18参考文献1

2、9一作业内容及要求1.1 作业内容用小波对图像进行滤波分解、去噪，然后重构。1.2 作业要求用小波对图像进行滤波分解、去噪，然后重构。具体要求：(1) 被处理图像可选择：woman, wbarb, wgatlin, detfingr, tire.；(2) 可以选择db等正交小波、或双正交小波（或用几种小波）；(3) 用选用小波的分解滤波器通过定义的卷积函数conv_my( )对图像二维数组进行小波分解，并进行下采样，获取CA、CV、CD、CH等分解子图；(4) 对高频信号子图进行去噪处理，可以采用软阈值、硬阈值等方法；(5) 用选用小波的综合滤波器对去噪的子图进行图像重构。二系统原理2.1 小

3、波变换原理小波变换的一级分解过程是，先将信号与低通滤波器卷积再下采样可以得到低频部分的小波分解系数再将信号与高通滤波器卷积后下采样得到高频部分的小波分解系数；而多级分解则是对上一级分解得到的低频系数再进行小波分解，是一个递归过程。二维小波分解重构可以用一系列的一维小波分解重构来实现。重构则是分解的逆过程，对低频系数、高频系数分别进行上抽样和低通、高通滤波处理。要注意重构时同一级的低频、高频系数的个数必须相等。2.2 阈值去噪原理图像去噪的方法是：（1）图像的小波分解。选择合适的小波函数以及合适的分解层数对图像进行分解。（2）对分解后的高频系数进行阈值处理。对分解的每一层，选择合适的阈值对该层的

4、水平、垂直和对角三个方向的高频系数进行阈值处理。（3）重构图像。根据小波分解的低频系数和经阈值量化处理后的高频系数进行图像重构。本设计采用软阈值去噪，其原理为：当小波系数的绝对值小于给定的阈值时，令其为0，大于阈值时，令其都减去阈值，即：小波阈值在去噪过程中起到决定性的作用。如果太小，那么施加阈值以后的小波系数中将包含过多的噪声分量，达不到去噪的目的；反之，如果太大，那么将去除一部分信号的分量，从而使由小波系数重建后的信号产生过大的失真。MATLAB中实现阈值获取的函数有ddenmp、thselect、wbmpen和wdcbm，本次设计中采用ddenmp函数进行阈值获取。调用格式为THR,SO

5、RH,KEEPAPP=ddencmp('den','wv',X)，函数ddencmp用于获取信号在消噪或压缩过程中的默认阈值。输入参数X为一维或二维信号；'den'表示进行去噪； 'wv'表示选择小波。返回值THR是返回的阈值；SORH是软阈值或硬阈值选择参数；KEEPAPP表示保存低频信号。三系统分析及设计3.1 图像分解图像分解在程序中分为两部分：一维分解以及二维分解。（1）一维分解先将信号与低通滤波器卷积再下采样可以得到低频部分的小波分解系数，再将信号与高通滤波器卷积后下采样得到高频部分的小波分解系数（2）对输入的信号，也就

6、是前面load图像之后得到的X矩阵，先对每一行进行一维分解，在对分解后得到的矩阵的每一列进行一维分解，最后得到一个矩阵X2=CACHCVCD；可以将CA、CH、 CV、 CD依次输出得到分解后的低频、垂直、水平、对角分量。原始信号分别进行低通、高通滤波，再分别对列进行二元下抽样，就得到低频、高频（也称为平均、细节）两部分系数；再对这两部分系数进行小波分解，步骤与第一级分解一致，最后得到cA：低频分量，cH：水平分量，cV：垂直分量，cD：对角线分量。流程图如下图1所示。列卷积XLo_DHi_D22Lo_DHi_DHi_DLo_D2222CACHCDCV行卷积行卷积列卷积列卷积列卷积图1 分解流

7、程图3.2 高频去噪本次设计采用软阈值去噪方法，将对分解出的cH、cV、cD进行阈值去噪，选用ddencmp函数获取去噪过程中的默认软阈值。其调用格式为：THR,SORH,KEEPAPP=ddencmp(IN1,IN2,X)。IN1='den'，IN2='wv'，X为一维或二维矩阵信号。在获取默认阈值后，对该阈值进行软阈值处理，调用到函数wthcoef2进行处理，其调用过程中关键部分是tmp = (abs(x)-t);tmp = (tmp+abs(tmp)/2; y = sign(x).*tmp;其中t代表前面获得的默认阈值，x代表小波系数的大小。得到的y即为软

8、阈值处理后的小波系数。3.3 图像重构重构的过程与分解过程类似，同样分为两部分：一维重构和二维重构。(1)一维重构：先对平均部分系数进行上采样，再与低通滤波器卷积。然后对细节部分选取重构所需的细节部分，长度与本层平均部分系数相同，再对其上采样后进行高通滤波器卷积。将两个得到的结果相加得到新的平均部分系数。重复以上操作知道细节部分长度小于平均部分长度。(2)二维重构：根据重构流程图可知，先对cA和cH进行重构，再对cV和cD重构，在该程序中所以在该程序中X2=CACHCVCD，所以先对行左右两半部分进行重构，重构后得到的矩阵在进行列上下两半部分进行重构，最终得到的重构图像。列卷积具体步骤是先将去

9、噪后的子图分别对行进行，再将低频分量和垂直分量进行二元上采样、低通滤波并合成；将水平分量和对角线分量进行二元上采样、高通滤波并合成，就得到低频、高频（也称为平均、细节）两部分系数；再对这两部分系数进行小波合成，步骤与第一级重构一致，最后通过wkeep选取我们所需的原始图像的大小，得到原始图像。该流程图如下图2所示。2Lo_RHi_RHi_RLo_R列卷积列卷积列卷积CACHCDCV22222Lo_RLo_R行卷积行卷积WkeepX图2 重构流程图四程序编写4.1 main函数clear all;% db4小波分解成低通和高通滤波器%figure;Lo_D,Hi_D,Lo_R,Hi_R=wfil

10、ters('db4');subplot(221);stem(Lo_D);grid on;title('分解低通滤波器');xlabel('系数');ylabel('Lo_D');subplot(222);stem(Hi_D);grid on;title('分解高通滤波器');xlabel('系数');ylabel('Hi_D');subplot(223);stem(Lo_R);grid on;title('综合低通滤波器');xlabel('系数');

11、ylabel('Lo_R');subplot(224);stem(Hi_R);grid on;title('综合高通滤波器');xlabel('系数');ylabel('Hi_R');% 导入和显示原始图像%load woman;photo=X;%X含有被装载的信号figure;colormap(map);%设置颜色映射到矩阵图image(photo);axis off; %关闭原始信号的坐标轴title('原始信号');% 用自定义的一维和二维函数分解图像% 得到低频部分（cA），水平部分（cH），垂直部分（cV

12、），对角线部分（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

13、_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'分别

14、取'h','v','d'代表三个方向%N为尺度向量，T为阈值向量，SORH取'h'和's'分别表示使用硬阈值或软阈值，下同nc1=wthcoef2('h',C1,S1,1,thr1,'s');nc1=wthcoef2('v',C1,S1,1,thr1,'s');nc1=wthcoef2('d',C1,S1,1,thr1,'s');new_cH=waverec2(nc1,S1,'sym5'); %对新的小波

15、分解结构(nc，s)进行重构new_cHH=wkeep(new_cH,L,'c');%提取重构的图像，'c'代表中心部分，但因为L为完整长度，所以用'l'或'r'也可以，下同%用小波函数sym5对cV进行1层小波分解C2,S2=wavedec2(cV,1,'sym5');%求取信号去噪的默认阈值、软阈值、并保留低频部分thr2,sorh2,keepapp2=ddencmp('den','wv',cV);%对三个方向高频系数进行软阈值处理nc2=wthcoef2('h'

16、,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');%求取信号去噪的默认

17、阈值、软阈值、并保留低频部分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)进

18、行重构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）进行重构% reconst

19、ruction为重构之后的图形%reconstruction = myidwt2(cA,new_cH,new_cV,new_cD,Lo_R,Hi_R,map);%4.2 分解函数二维分解函数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：垂直方向分解

20、系数；% LH：水平方向分解系数； HH：对角线方向分解系数。row,col=size(x); % 读取输入矩阵的大小for j=1:row % 首先对输入矩阵的每一行序列进行一维离散小波分解 tmp1=x(j,:); ca1,cd1,n1,n2=mydwt(tmp1,Lo_D,Hi_D,1); a1=n1+n2; x1(j,:)=ca1,cd1; % 将分解系数序列再存入矩阵x1中，得到L|HendL=size(x1); %取x1的长度for k=1:a1 % 再对输入矩阵的每一列序列进行一维离散小波分解 tmp2=x1(:,k); ca2,cd2,n1,n2=mydwt(tmp2,Lo_D

21、,Hi_D,1); x2(:,k)=ca2,cd2; % 将分解所得系数存入矩阵x2中，得到LL,HL;LH,HHendrow,col=size(x2); LL=x2(1:row/2,1:col/2); % LL是矩阵x2的左上角部分LH=x2(row/2+1:row,1:col/2); % LH是矩阵x2的左下角部分HL=x2(1:row/2,col/2+1:col); % HL是矩阵x2的右上角部分HH=x2(row/2+1:row,col/2+1:col); % HH是矩阵x2的右下角部分figure;subplot(221);colormap(map);image(LL);title(

22、'分解后的低频分量');subplot(222);colormap(map);image(LH);title('分解后的垂直分量');subplot(223);colormap(map);image(HL);title('分解后的水平分量');subplot(224);colormap(map);image(HH);title('分解后的对角分量');一维分解函数function cA,cD,n1,n2 = mydwt(x,Lo_D,Hi_D,dim)% 函数 cA,cD,n1,n2=MYDWT(X,LPD,HPD,DIM) 对输

23、入序列x进行一维离散小波分解，输出分解序列cA,cD% 输入参数：x输入序列；% lpd低通滤波器；% hpd高通滤波器；% dim小波分解级数。% 输出参数：cA低频部分的小波分解系数；% cD高频部分的小波分解系数。cA=x; % 初始化cA，cDcD=;for i=1:dimcvl=myconv1(cA,Lo_D); % 低通滤波dnl=downspl(cvl); % 通过下采样求出平均部分的分解系数 n1=length(dnl);cvh=myconv1(cA,Hi_D); % 高通滤波dnh=downspl(cvh); % 通过下采样求出本层分解后的细节部分系数 n2=length(d

24、nh);cA=dnl; % 下采样后的平均部分系数进入下一层分解cD=cD,dnh; % 将本层分解所得的细节部分系数存入序列cDend4.3 卷积函数function c=myconv1(a,b);m=length(a);%获取a矩阵的长度n=length(b);%获取b矩阵的长度for i=1:m+n-1c(i)=0;%定义长度为m+n-1的全为0的矩阵cendfor i=1:m for j=1:n x=a(i)*b(j);%将a矩阵中的数乘以b矩阵中的数c(i+j-1)=c(i+j-1)+x;%将上一行得到的数累加到对应序号的c矩阵数中endend4.4 采样函数4.下采样函数funct

25、ion y=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.上采样函数function y=upspl(x)% 函数 Y = UPSPL(X) 对输入的一维序列x进行上抽样，即对序列x每个元素之间% 插零，例如 x=x1,x2,x3,x4,上抽样后为 y=x1,0,x2,0,

26、x3,0,x4;N=length(x); % 读取输入序列长度M=2*N-1; % 输出序列的长度是输入序列长度的2倍再减一for i=1:M % 输出序列的偶数位为0，奇数位按次序等于相应位置的输入序列元素if mod(i,2) y(i)=x(i+1)/2);else y(i)=0;endend4.5 重构函数二维重构函数function y=myidwt2(LL,HL,LH,HH,Lo_R,Hi_R,map)% 函数 MYIDWT2() 对输入的子矩阵序列进行逆小波变换，重构出矩阵 y% 输入参数：LL,HL,LH,HH 是四个大小均为 r*c 维的矩阵% 输出参数：y 是一个大小为 2r

27、*2c维的矩阵tmp_mat=LL,HL;LH,HH; % 将输入的四个矩阵组合为一个矩阵flag=0; %循环标志位for i=1:2if flag=1 row,col=size(yt); % 求出组合矩阵的行列数for k=1:col % 首先对组合矩阵tmp_mat的每一列，分开成上下两半 ca1=yt(1:row/2,k); % 分开的两部分分别作为平均系数序列ca1、细节系数序列cd1 cd1=yt(row/2+1:row,k); tmp1=myidwt(ca1,cd1,Lo_R,Hi_R);% 重构序列 yt1(:,k)=tmp1; % 将重构序列存入待输出矩阵yt的相应列，此时

28、y=L|Hendbreak;endif flag=0 row,col=size(tmp_mat); for j=1:row % 将输出矩阵 y 的每一行，分开成左右两半 ca2=tmp_mat(j,1:col/2); % 分开的两部分分别作为平均系数序列ca2、细节系数序列cd2 cd2=tmp_mat(j,col/2+1:col); tmp2=myidwt(ca2,cd2,Lo_R,Hi_R); % 重构序列yt(j,:)=tmp2; % 将重构序列存入待输出矩阵yt的相应行，得到最终的输出矩阵 y=ytend flag=1;endendy=yt1;L=256,256; %取数据原始大小fi

29、gure;y=wkeep(yt1,L,'c') colormap(map);image(y);title('重构图像');一维重构函数function y=myidwt(cA,cD,Lo_R,Hi_R)% 函数 MYIDWT() 对输入的小波分解系数进行逆离散小波变换，重构出信号序列 y% 输入参数：cA 平均部分的小波分解系数；% cD 细节部分的小波分解系数；% Lo_R、Hi_R 重构所用的低通、高通滤波器。lca=length(cA); % 求出平均、细节部分分解系数的长度lcd=length(cD);while (lcd)>=(lca) % 每一层重构中，cA和cD的长度要相等，故每层重构后，% 若lcd小于lc