Matlab小波分析工具箱丰富的函数和强大的仿真功能为我们学习小波、用好小波提供了方便、快捷的途径，但是，如果我们要深入掌握小波分析的原理，真正学好、用好小波，就应该尽量用自己编写的程序去实现小波变换和信号分析，尽量在自己的程序中少调用Matlab提供的函数，多用自己的理解去编写相关的小波函数，这样的过程是一个探索、求知的过程，更能让我们体会到小波的强大和学习的乐趣。下面，我把自己编写的小波一维、二维信号分解和重构Matlab程序共享出来，也希望有朋友共享自编的程序，共同学习，提高程序的效率和简洁性。
首先要说明的一点是，虽然是自己编写Matlab程序，但并不是说一点也不用Matlab的自带函数。我们要编写的是实现小波变换的主要功能函数，而绘图等基本功能还是要用到Matlab函数的。而且，根据小波变换的滤波器组原理，原始信号要通过低通、高通滤波器处理，这里就涉及到卷积这一运算步骤。卷积——FFT算法的实现，相信很多朋友都能用Matlab、C语言等来实现，不过与Matlab自带的用机器语言编写的FFT程序相比，运算速度一般会慢几倍、几十倍。所以，我的程序里边涉及卷积的就直接调用Matlab的conv()函数了。
我们知道，小波变换的一级分解过程是，原始信号分别进行低通、高通滤波，再分别进行二
元下抽样，就得到低频、高频（也称为平均、细节）两部分系数；而多级分解则是对上一级分解得到的低频系数再进行小波分解，是一个递归过程。以下是一维小波分解的程序：
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=[];
fori=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
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);
而重构则是分解的逆过程，对低频系数、高频系数分别进行上抽样和低通、高通滤波处理。要注意重构时同一级的低频、高频系数的个数必须相等。
functiony=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等于重构完成后的平均部分系数序列c
A
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;???????????%输出序列的长度是输入序列长度的2倍再减一
fori=1:M??????????%输出序列的偶数位为0，奇数位按次序等于相应位置的输入序列元素
???ifmod(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/21/2];hpd=[-1/21/2];??????????%默认的低通、高通滤波器
[row,col]=size(x);?????????????????????%读取输入矩阵的大小
forj=1:row????????????????????????????%首先对输入矩阵的每一行序列进行一维离散小波分解
???tmp1=x(j,:);
???[ca1,cd1]=mydwt(tmp1,lpd,hpd,1);
???x(j,:)=[ca1,cd1];??????????????????%将分解系数序列再存入矩阵x中，得到[L|H]
end
fork=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的右下角部分
functiony=myidwt2(LL,HL,LH,HH);
%函数MYIDWT2()对输入的子矩阵序列进行逆小波变换，重构出矩阵y
%输入参数：LL,HL,LH,HH——是四个大小均为r*c维的矩阵
%输出参数：y——是一个大小为2r*2c维的矩阵
lpr=[11];hpr=[1-1];??????????%默认的低通、高通滤波器
tmp_mat=[LL,HL;LH,HH];?????????%将输入的四个矩阵组合为一个矩阵
[row,col]=size(tmp_mat);???????%求出组合矩阵的行列数
matlab学好了有什么用
fork=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
forj=1:row????????????????????%将输出矩阵y的每一行，分开成左右两半
???ca2=yt(j,1:col/2);??????????%分开的两部分分别作为平均系数序列ca2、细节系数序列cd2

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。