2022年自己动手编写小波信号分解与重构的Matlab程序 .pdf

去年 11 月发布了一系列有关小波变换和图像处理的文章，把学习小波过程中的心得体会和编写的程序放在网上和大家共享交流。半年来，感谢大家的关注和帮助，在相互的讨论交流中， 我不断地从大家提出的问题中拓展自己的知识面，对小波的理论及其应用有了更深入的了解和掌握。根据和大家讨论交流中发现的问题，对博客中的程序进行修正。有关小波图像分解和重构的两篇文章中分享的程序，存在下列问题：（1）程序所用的小波函数只有非标准的Haar 小波，其滤波器组为Lo_D=1/2 1/2, Hi_D=-1/2 1/2，是固化在mydwt2.m 的程序中的，不能选择其他的小波函数；（2）非标准的Haar 小波，其分解出来的系数矩阵中，高频系数的细节内容（轮廓、边缘等特征）不明显；（3）函数mydwt2 中列变换的矩阵对象为输入矩阵，这是错误的，其矩阵对象应该是行变换后的缓存矩阵；（4）函数mydwt2 的输出用 LL,HL,LH,HH表示，不是很规范，应改为cA,cV,cH,cD来表示，即一级小波变换输出的系数矩阵有4 个部分：平均部分、垂直细节部分、水平细节部分和对角线细节部分。（5）函数mywavedec2 的输出y 是与输入矩阵x 相同大小的矩阵，并且已将N 级分解后所有的平均、细节系数组合成一体的。实际上，这种定义只对Haar 小波有效。（6）原程序中要调用modmat 函数对图像矩阵进行修剪，使之能被2 的 N 次方整除，主要是为了生成塔式结构图像而设的，对上述问题修正后，这个modmat 函数已不需使用了。针对上述问题，我对程序作了修正，发布在今天的3 篇文章里，请大家点击查看。新修正的程序更为简洁易懂，功能也有所增强，可以用任意的小波函数进行小波分解，可根据小波分解系数矩阵重构出指定分解级的低频系数和原始图像. Matlab 小波分析工具箱丰富的函数和强大的仿真功能为我们学习小波、用好小波提供了方便、快捷的途径，但是，如果我们要深入掌握小波分析的原理，真正学好、用好小波，就应该尽量用自己编写的程序去实现小波变换和信号分析，尽量在自己的程序中少调用Matlab 提供的函数，多用自己的理解去编写相关的小波函数，这样的过程是一个探索、求知的过程，更能让我们体会到小波的强大和学习的乐趣。下面，我把自己编写的小波一维、二维信号分解和重构Matlab 程序共享出来，也希望有朋友共享自编的程序，共同学习，提高程序的效率和简洁性。首先要说明的一点是，虽然是自己编写Matlab 程序，但并不是说一点也不用Matlab的自带函数。 我们要编写的是实现小波变换的主要功能函数，而绘图等基本功能还是要用到Matlab 函数的。而且，根据小波变换的滤波器组原理，原始信号要通过低通、高通滤波器处理，这里就涉及到卷积这一运算步骤。卷积 FFT 算法的实现，相信很多朋友都能用Matlab 、C 语言等来实现，不过与Matlab 自带的用机器语言编写的FFT 程序相比，运算速度一般会慢几倍、几十倍。所以，我的程序里边涉及卷积的就直接调用Matlab 的 conv() 函数了。名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 1 页，共 6 页 - - - - - - - - - 我们知道，小波变换的一级分解过程是，原始信号分别进行低通、高通滤波，再分别进行二元下抽样，就得到低频、高频（也称为平均、细节）两部分系数；而多级分解则是对上一级分解得到的低频系数再进行小波分解，是一个递归过程。 以下是一维小波分解的程序：function cA,cD = mydwt(x,lpd,hpd,dim); % 函数cA,cD=MYDWT(X,LPD,HPD,DIM) 对输入序列x 进行一维离散小波分解，输出分解序列 cA,cD % 输入参数： x 输入序列；% lpd 低通滤波器；% hpd 高通滤波器；% dim 小波分解级数。% 输出参数： cA 平均部分的小波分解系数；% cD 细节部分的小波分解系数。cA=x; % 初始化 cA，cD cD=; for i=1:dim cvl=conv(cA,lpd); % 低通滤波， 为了提高运行速度，调用 MATLAB 提供的卷积函数 conv() dnl=downspl(cvl); % 通过下抽样求出平均部分的分解系数cvh=conv(cA,hpd); % 高通滤波dnh=downspl(cvh); % 通过下抽样求出本层分解后的细节部分系数cA=dnl; % 下抽样后的平均部分系数进入下一层分解cD=cD,dnh; % 将本层分解所得的细节部分系数存入序列cD end function y=downspl(x); % 函数Y=DOWMSPL(X) 对输入序列进行下抽样，输出序列Y。% 下抽样是对输入序列取其偶数位，舍弃奇数位。例如x=x1,x2,x3,x4,x5，则y=x2,x4. N=length(x); % 读取输入序列长度名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 2 页，共 6 页 - - - - - - - - - M=floor(N/2); % 输出序列的长度是输入序列长度的一半（带小数时取整数部分）i=1:M; y(i)=x(2*i); 而重构则是分解的逆过程，对低频系数、高频系数分别进行上抽样和低通、高通滤波处理。要注意重构时同一级的低频、高频系数的个数必须相等。function y = myidwt(cA,cD,lpr,hpr); % 函数MYIDWT() 对输入的小波分解系数进行逆离散小波变换，重构出信号序列y % 输入参数： cA 平均部分的小波分解系数；% cD 细节部分的小波分解系数；% lpr、hpr 重构所用的低通、高通滤波器。lca=length(cA); % 求出平均、细节部分分解系数的长度lcd=length(cD); while (lcd)=(lca) % 每一层重构中，cA 和 cD 的长度要相等，故每层重构后，% 若 lcd 小于 lca，则重构停止， 这时的cA 即为重构信号序列y 。upl=upspl(cA); % 对平均部分系数进行上抽样cvl=conv(upl,lpr); % 低通卷积cD_up=cD(lcd-lca+1:lcd); % 取出本层重构所需的细节部分系数，长度与本层平均部分系数的长度相等uph=upspl(cD_up); % 对细节部分系数进行上抽样cvh=conv(uph,hpr); % 高通卷积cA=cvl+cvh; % 用本层重构的序列更新cA，以进行下一层重构cD=cD(1:lcd-lca); % 舍弃本层重构用到的细节部分系数，更新cD lca=length(cA); % 求出下一层重构所用的平均、细节部分系数的长度lcd=length(cD); end % lcd lca ，重构完成，结束循环y=cA; % 输出的重构序列y 等于重构完成后的平均部分系数序列cA 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 3 页，共 6 页 - - - - - - - - - function y=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; % 输出序列的长度是输入序列长度的2 倍再减一for i=1:M % 输出序列的偶数位为0，奇数位按次序等于相应位置的输入序列元素if mod(i,2) y(i)=x(i+1)/2); else y(i)=0; end end 我们知道，二维小波分解重构可以用一系列的一维小波分解重构来实现。以下程序是基于 Haar 小波的二维小波分解和重构过程：function LL,HL,LH,HH=mydwt2(x); % 函数MYDWT2() 对输入的r*c 维矩阵x 进行二维小波分解，输出四个分解系数子矩阵 LL,HL,LH,HH % 输入参数： x 输入矩阵，为r*c 维矩阵。% 输出参数： LL,HL,LH,HH 是分解系数矩阵的四个相等大小的子矩阵,大小均为r/2 * c/2 维% LL：低频部分分解系数；HL：垂直方向分解系数；% LH：水平方向分解系数；HH：对角线方向分解系数。lpd=1/2 1/2;hpd=-1/2 1/2; % 默认的低通、高通滤波器row,col=size(x); % 读取输入矩阵的大小for j=1:row % 首先对输入矩阵的每一行序列进行一维离散小波分解tmp1=x(j,:); ca1,cd1=mydwt(tmp1,lpd,hpd,1); 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 4 页，共 6 页 - - - - - - - - - x(j,:)=ca1,cd1; % 将分解系数序列再存入矩阵x 中，得到 L|H end for k=1:col % 再对输入矩阵的每一列序列进行一维离散小波分解tmp2=x(:,k); ca2,cd2=mydwt(tmp2,lpd,hpd,1); x(:,k)=ca2,cd2; % 将分解所得系数存入矩阵x 中，得到 LL,Hl;LH,HH end LL=x(1:row/2,1:col/2); % LL 是矩阵 x 的左上角部分LH=x(row/2+1:row,1:col/2); % LH 是矩阵 x 的左下角部分HL=x(1:row/2,col/2+1:col); % HL 是矩阵 x 的右上角部分HH=x(row/2+1:row,col/2+1:col); % HH 是矩阵 x 的右下角部分function y=myidwt2(LL,HL,LH,HH);% 函数MYIDWT2() 对输入的子矩阵序列进行逆小波变换，重构出矩阵y % 输入参数： LL,HL,LH,HH 是四个大小均为r*c 维的矩阵% 输出参数： y 是一个大小为2r*2c 维的矩阵lpr=1 1;hpr=1 -1; % 默认的低通、高通滤波器tmp_mat=LL,HL;LH,HH; % 将输入的四个矩阵组合为一个矩阵row,col=size(tmp_mat); % 求出组合矩阵的行列数for k=1:col % 首先对组合矩阵tmp_mat 的每一列，分开成上下两半ca1=tmp_mat(1:row/2,k); % 分开的两部分分别作为平均系数序列ca1、细节系数序列 cd1 cd1=tmp_mat(row/2+1:row,k); tmp1=myidwt(ca1,cd1,lpr,hpr); % 重构序列yt(:,k)=tmp1; % 将重构序列存入待输出矩阵yt 的相应列，此时y=L|H end for j=1:row % 将输出矩阵y 的每一行，分开成左右两半ca2=yt(j,1:col/2); % 分开的两部分分别作为平均系数序列ca2、细节系数序列 cd2 cd2=yt(j,col/2+1:col); tmp2=myidwt(ca2,cd2,lpr,hpr); % 重构序列名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 5 页，共 6 页 - - - - - - - - - yt(j,:)=tmp2; % 将重构序列存入待输出矩阵yt 的相应行，得到最终的输出矩阵y=yt end y=yt; 我编写程序时尽量实现模块化，方便不同功能程序之间的调用。以上程序实现了一维、二维小波分解与重构的功能，接下来我们就可以此为基础，探讨小波信号压缩和图像融合等技术了。名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 6 页，共 6 页 - - - - - - - - -