《现代信号处理》大作业基于Matlab的小波分解、去噪与重构目录一作业内容及要求 (3)1.1 作业内容 (3)1.2 作业要求 (3)二系统原理 (3)2.1 小波变换原理 (3)2.2 阈值去噪原理 (3)三系统分析及设计 (5)3.1 图像分解 (5)3.2 高频去噪 (5)3.3 图像重构 (6)四程序编写 (7)4.1 main函数 (7)4.2 分解函数 (9)4.2.1 二维分解函数 (9)4.2.2 一维分解函数 (10)4.3 卷积函数 (10)4.4 采样函数 (11)4.4.1 下采样函数 (11)4.4.2 上采样函数 (11)4.5 重构函数 (12)4.5.1 二维重构函数 (12)4.5.2 一维重构函数 (13)五结果分析及检验 (14)5.1 结果分析 (14)5.2 结果检验 (16)六心得体会 (18)参考文献 (19)一作业内容及要求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,大于阈值时,令其都减去阈值,即:λλλλ<>-=||,0| |),|)](|([{ww wwwsign小波阈值λ在去噪过程中起到决定性的作用。
如果λ太小,那么施加阈值以后的小波系数中将包含过多的噪声分量,达不到去噪的目的;反之,如果λ太大,那么将去除一部分信号的分量,从而使由小波系数重建后的信号产生过大的失真。
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 矩阵,先对每一行进行一维分解,在对分解后得到的矩阵的每一列进行一维分解,最后得到一个矩阵X2=[CA CH CV CD];可以将CA 、CH 、 CV 、 CD 依次输出得到分解后的低频、垂直、水平、对角分量。
原始信号分别进行低通、高通滤波,再分别对列进行二元下抽样,就得到低频、高频(也称为平均、细节)两部分系数;再对这两部分系数进行小波分解,步骤与第一级分解一致,最后得到cA :低频分量,cH :水平分量,cV :垂直分量,cD :对角线分量。
流程图如下图1所示。
图1 分解流程图 3.2 高频去噪本次设计采用软阈值去噪方法,将对分解出的cH 、cV 、cD 进行阈值去噪,选用ddencmp 函数获取去噪过程中的默认软阈值。
其调用格式为:[THR,SORH,KEEPAPP]=ddencmp(IN1,IN2,X)。
IN1='den',IN2='wv',X 为一维或二维矩阵信号。
在获取默认阈值后,对该阈值进行软阈值处理,调用到函数wthcoef2进行XLo_D Hi_D 2 2 Lo_DHi_D Hi_DLo_D2 2 2 2 CA CH CDCV 行卷积行卷积 列卷积 列卷积 列卷积列卷积处理,其调用过程中关键部分是tmp = (abs(x)-t);tmp = (tmp+abs(tmp))/2;y = sign(x).*tmp;其中t 代表前面获得的默认阈值,x 代表小波系数的大小。
得到的y 即为软阈值处理后的小波系数。
3.3 图像重构重构的过程与分解过程类似,同样分为两部分:一维重构和二维重构。
(1)一维重构:先对平均部分系数进行上采样,再与低通滤波器卷积。
然后对细节部分选取重构所需的细节部分,长度与本层平均部分系数相同,再对其上采样后进行高通滤波器卷积。
将两个得到的结果相加得到新的平均部分系数。
重复以上操作知道细节部分长度小于平均部分长度。
(2)二维重构:根据重构流程图可知,先对cA 和cH 进行重构,再对cV 和cD 重构,在该程序中所以在该程序中X2=[CA CH CV CD],所以先对行左右两半部分进行重构,重构后得到的矩阵在进行列上下两半部分进行重构,最终得到的重构图像。
具体步骤是先将去噪后的子图分别对行进行,再将低频分量和垂直分量进行二元上采样、低通滤波并合成;将水平分量和对角线分量进行二元上采样、高通滤波并合成,就得到低频、高频(也称为平均、细节)两部分系数;再对这两部分系数进行小波合成,步骤与第一级重构一致,最后通过wkeep 选取我们所需的原始图像的大小,得到原始图像。
该流程图如下图2所示。
图2 重构流程图2 Lo_R Hi_RHi_RLo_R 列卷积 列卷积列卷积 CA CH CD CV2 2 2 2 2 Lo_R Lo_R 行卷积行卷积 Wkeep X’列卷积四程序编写4.1 main函数clear all; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % db4小波分解成低通和高通滤波器% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figure;[Lo_D,Hi_D,Lo_R,Hi_R]=wfilters('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('系数');ylabel('Lo_R');subplot(224);stem(Hi_R);grid on;title('综合高通滤波器');xlabel('系数');ylabel('Hi_R'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 导入和显示原始图像% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% load woman;pho to=X; %X含有被装载的信号figure;colormap(map);%设置颜色映射到矩阵图image(pho to);axis off; %关闭原始信号的坐标轴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'分别表示使用硬阈值或软阈值,下同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'); %对新的小波分解结构(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',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维矩阵。