小波分析工具箱丰富的功能和强大的仿真功能为我们学习和使用好小波提供了方便快捷的方法。 但是,如果我们想深入掌握小波分析的原理,真正学好小波、用好小波,还得尽力而为。 使用您自己的程序来实现小波变换和信号分析。 尽量在自己的程序中少调用提供的函数,用自己的理解来编写相关的小波函数。 这个过程是一个探索、求知的过程,可以让人更加舒服。 我们意识到小波的力量和学习的乐趣。 下面分享一下我写的小波一维和二维信号分解重构程序。 也希望朋友们能够分享自己写的程序,共同学习,提高程序的效率和简洁性。
首先要注意的是,虽然你自己写了程序,但并不意味着你根本不用它。 我们要写的是实现小波变换的主要功能函数,绘图等基本功能还是需要使用函数。 而且,根据小波变换的滤波器组原理,原始信号必须经过低通和高通滤波器的处理,这涉及到卷积步骤。 卷积 - FFT 算法的实现。 相信很多朋友都会用,C语言等,但是与机器语言编写的内置FFT程序相比,运算速度一般会慢几倍、几十倍。 所以我的程序中涉及到卷积的时候就直接调用了conv()函数。
我们知道,小波变换的一级分解过程是对原始信号分别进行低通和高通滤波,然后分别进行二值下采样,得到低频和高频(也称为平均和细节)系数; 而多级分解就是对上一级分解得到的低频系数进行小波分解,这是一个递归过程。 以下是一维小波分解的过程:
[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=[];
对于 i=1:暗淡
cvl=conv(ca,lpd); % 低通滤波,为了提高运行速度,调用提供的卷积函数conv()
dnl=(cvl); % 通过下采样求出平均部分的分解系数
CVH = CV(CA,HPD); %高通滤波
dnh=(cvh); % 通过下采样求分解后本层的细节系数
cA=dnl; % 下采样后的平均偏系数进入下一级分解
cD=[cD,dnh]; % 将本层分解得到的细节系数存储到序列cD中
结尾
y=(x);
% 函数 Y=(X) 对输入序列进行下采样并输出序列 Y。
% 下采样就是取输入序列的偶数位,丢弃奇数位。 例如,x=[x1,x2,x3,x4,x5],则 y=[x2,x4]。
N=(x); % 读取输入序列长度
M=下限(N/2); % 输出序列的长度为输入序列长度的一半(有小数时取整数部分)
我=1:M;
y(i)=x(2*i);
重构是分解的逆过程,分别对低频系数和高频系数进行上采样、低通和高通滤波。 需要注意的是,重建时同一层的低频和高频系数的数量必须相等。
y = (ca,cD,lpr,hpr);
% ()对输入的小波分解系数进行逆离散小波变换,重构信号序列y
% 输入参数:cA——平均部分的小波分解系数;
%cD——细节部分的小波分解系数;
% lpr、hpr - 重建中使用的低通和高通滤波器。
lca=(ca); % 求平均和细节分解系数的长度
液晶屏=(CD);
while (lcd)>=(lca) % 在每一层重构中,cA和cD的长度必须相等,所以在每一层重构之后,
% 如果lcd小于lca,则重构停止,此时的cA就是重构后的信号序列y。
upl=upspl(cA); % 对平均部分系数进行上采样
cvl=转换(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=(ca); % 求下一层重建时使用的均值和细节系数的长度
液晶屏=(CD);
end %lcd < lca,重建完成,结束循环
y = cA; % 输出重构序列y等于重构完成后的平均部分系数序列cA
y=uppl(x);
% 函数 Y = UPSPL(X) 对输入一维序列 x 进行上采样,即序列 x 的各个元素之间
% 零插值,例如x=[x1,x2,x3,x4],上采样后为y=[x1,0,x2,0,x3,0,x4];
N=(x); % 读取输入序列长度
M=2*N-1; % 输出序列的长度是输入序列长度的2倍减1。
for i=1:M % 输出序列的偶数位为0,奇数位等于输入序列按顺序对应位置的元素。
如果 mod(i,2)
y(i)=x((i+1)/2);
别的
y(i)=0;
结尾
结尾
我们知道,二维小波分解和重构可以通过一系列一维小波分解和重构来实现。 下面的程序是基于Haar小波的二维小波分解和重构过程:
[LL,HL,LH,HH]=(x);
% ()对输入的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]; % 默认低通和高通滤波器
[行,列]=大小(x); % 读取输入矩阵的大小
for j=1:row % 首先对输入矩阵的每一行序列进行一维离散小波分解
tmp1=x(j,:);
[ca1,cd1]=mydwt(tmp1,lpd,hpd,1);
x(j,:)=[ca1,cd1]; % 将分解系数序列存入矩阵x,得到[L|H]
结尾
for k=1:col % 然后对输入矩阵的各列序列进行一维离散小波分解
tmp2=x(:,k);
[ca2,cd2]=mydwt(tmp2,lpd,hpd,1);
x(:,k)=[ca2,cd2]; % 将分解得到的系数存入矩阵x,得到[LL,Hl;LH,HH]
结尾
LL=x(1:行/2,1:列/2); % LL是矩阵x的左上角
LH=x(行/2+1:行,1:列/2); %LH为矩阵x的左下角
HL=x(1:行/2,列/2+1:列); % HL为矩阵x的右上角
HH=x(行/2+1:行,列/2+1:列); %HH是矩阵x的右下角
y=(LL,HL,LH,HH);
% ()对输入的子矩阵序列进行逆小波变换,重构矩阵y
% 输入参数:LL、HL、LH、HH - 是四个 r*c 维度的矩阵
% 输出参数:y——是2r*2c维度的矩阵
lpr=[1 1];hpr=[1 -1]; % 默认低通和高通滤波器
=[LL,HL;LH,HH]; % 将四个输入矩阵合并为一个矩阵
[行,列]=大小(); % 求组合矩阵的行数和列数
for k=1:col % 首先将组合矩阵的每一列分成上半部和下半部
ca1=(1:行/2,k); % 两个独立的部分分别是平均系数序列ca1和细节系数序列cd1
cd1=(行/2+1:行,k);
tmp1=(ca1,cd1,lpr,hpr); %重建序列
yt(:,k)=tmp1; % 将重构后的序列存放在要输出的矩阵yt的对应列中,此时y=[L|H]
结尾
for j=1:row % 将输出矩阵 y 的每一行分为左右两半
ca2=yt(j,1:列/2); % 分开的两个部分分别作为平均系数序列ca2和细节系数序列cd2
cd2=yt(j,列/2+1:列);
tmp2=(ca2,cd2,lpr,hpr); %重建序列
yt(j,:)=tmp2; % 将重构后的序列存放在要输出的矩阵yt的对应行中,得到最终的输出矩阵y=yt
结尾
y=yt;