小波信号分解与重构的Matlab程序
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
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=[];
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); % 读取输入序列长度
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
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);
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的右下角部分