Mallat算法小结含MATLAB编程

Mallat算法小结一、发展背景1987年，MallatandMeyer提出了多分辨分析（MultiresolutionAnalysis）,从而成功的统一了在此之前的Stromberg、Meyer、Lemarie和Battle提出的具体小波函数的构造，研究了小波变换的离散化情形。并且Mallat在图像的分解与重构的塔式算法的启发下，根据多分辨率理论，提出小波分解与重构的快速算法，称为Mallat算法。二、算法优点Mallat算法可以避免尺度a值越大，对信息Ψ(t)的采样就得越密的缺点，这一算法在小波分析中地位很重要，相当于快速傅里叶变换(FFT)在经典傅里叶分析中的地位。三、基本原理及步骤1、Mallat算法的信号分解1）根据多分辨率理论，得出结论：(j)(t)=(t)jjnnnPfx其中，Pjf(t):f(t)在Vj中的投影，是f(t)在分辨率j下的平滑逼近。Xn(j):线性组合的权重,也就是f(t)在分辨率j下的平滑逼近的概貌。jn(t):离散后的正交小波基当j=0时，P0f(t)=(0)0nnntx(1)由于(1)1101=p(t),(t)=p(t),(t)kkkffx因为D1f(t)与1(t)k正交，所以01(t),(t)0kDf,所以(1)(0)(0)0101=(t),(t)=(t),(t)nknkkknnnxxx(2)其中010(n-2k)(t),(t)=nkh代入（2）得(1)(0)0(n-2k)=knnhxx（3）注：(1)kx为离散平滑逼近类似还可以得到(1)(0)1(n-2k)=knnhdx（4）注：(1)kd为离散细节信号即小波信号2）从设计滤波器的角度考虑，设'(0)''(0)00=*(k)=(k-n)knnxxhhx（5）经过二抽取后，得(1)'2=kkxx将（5）式代入上式，得(1)(0)0=(n-2k)knnxhx（6）类似还可以得到(1)(0)1=(n-2k)knndhx（7）注：（6）（7）式表现了由V0到V1，W1的分解。仿照前面的步骤推导出如下结论(2)(1)(1)1201=(t),(t)=(t),(t)nknkknnnnxxx（8）(2)(1)(1)1201=(t),(t)=(t),(t)nknkknnnndxx（9）其中分解系数为1n2k0n1k0(n-2k)(t),(t)=(t),(t)=h12011(n-2k)(t),(t)=(t),(t)=nknkh这样我们可以逐级引申，对(1)kx做由V1到V2，W2的分解，得到(2)kx和(2)kd，再对(2)kx做由V2到V3，W3的分解，得到(3)kx和(3)kd，……，对(j)kx由Vj到Vj+1，Wj+1的分解，所需的电路结构不变，且滤波器的系数仍为0(-k)h='0(k)h，1(-k)h='1(k)h。3）如下图所示的网络结构可以重复推演下去。图1网络级联结构2、Mallat算法的信号重建用类似的思路，可以逆推重建过程，由上可知Vj=Vj+1Wj+1所以-1(t)=P(t)+D(t)jjjPfff(j)(j)=(t)+(t)jkjkkkKKxd又(j-1)-1-1,n=P(t),(t)jjnfx因此(j-1)(j)(j)-1,n-1,n=(t),(t)+(t),(t)jkjjkjnkkKKxxd由Mallat算法的信号分解中相同的证明，我们可以得到-1,n10n0(n-2k)(t),(t)=(t),(t)=jkjkh-1,n10n1(n-2k)(t),(t)=(t),(t)=jkjkh所以(j)(j-1)(j)1(n-2k)=+0(n-2k)nkkKKhhxxd(j)(j)1=+(n-2k)0(n-2k)kkKKggxd（10）（10）式反映了相邻两级的反演关系，其中(j)nx是第j级的离散平滑信号，(j)nd是第j级的离散细节信号；(-1)jnx是由(j)nx和(j)nd重建得到的第j-1级离散平滑信号。这里0100k(k)=(t),(t)g1100k(k)=(t),(t)g其中0(k)g、1(k)g与前面的0(k)h、1(k)h一样，为重建系数。下图为信号重建的网络结构。图2信号重建的网络结构四、Mallat算法分解与重建的比较（1）Mallat算法的信号分解系数0h、1h相当于分析滤波器；而Mallat算法的信号重建系数0g、1g相当于综合滤波器；（2）在重建式（10）中，是对k求和，而在分解式（3）和（4）中，是对n求和；（3）在分解算法中信号是先滤波后抽取，而在重建算法中是先插值后滤波。五、Mallat算法的MATLAB算法编程举例：应用Mallat算法实现两个正弦混合波的分解与重构,分别采用db30小波和db4作为小波函数，分解阶数分别取4和6，并将结果比较分析。MATLAB程序：（1）%n是阶数n=6，小波函数为db30clc;clear;%%1信号源：混合正弦波f1=50;%频率1f2=100;%频率2n=6;%n是阶数t=0:1/(2^(n)-1):1;N=length(t);y=sin(2*pi*f1*t)+sin(2*pi*f2*t);%正弦波混合figure(1)plot(y);xlabel('t=0:1')title('两个正弦信号的混合正弦波')figure(2)stem(abs(fft(y)));title('混合信号频谱')%%2.小波滤波器的谱分析h=wfilters('db30','l');%低通g=wfilters('db30','h');%高通h=[h,zeros(1,N-length(h))];%补零（圆周卷积，且增大分辨率变于观察）g=[g,zeros(1,N-length(g))];%补零（圆周卷积，且增大分辨率变于观察）figure(3);%滤波器图stem(abs(fft(h)));title('低通滤波器图')figure(4);stem(abs(fft(g)));title('高通滤波器图')%%3.MALLET分解算法(圆周卷积的快速傅里叶变换实现)sig1=ifft(fft(y).*fft(h));%低通(低频分量)sig2=ifft(fft(y).*fft(g));%高通(高频分量)figure(5);%分解信号图subplot(2,1,1)plot(real(sig1));title('分解信号1')subplot(2,1,2)plot(real(sig2));figure(6);%分解信号频谱图subplot(2,1,1)stem(abs(fft(sig1)));title('分解信号1频谱')subplot(2,1,2)stem(abs(fft(sig2)));title('分解信号2频谱')%%4.MALLET重构算法sig1=dyaddown(sig1);%2抽取sig2=dyaddown(sig2);%2抽取sig1=dyadup(sig1);%2插值sig2=dyadup(sig2);%2插值sig1=sig1(1,[1:N]);%去掉最后一个零sig2=sig2(1,[1:N]);%去掉最后一个零hr=h(end:-1:1);%重构低通gr=g(end:-1:1);%重构高通hr=circshift(hr',1)';%位置调整圆周右移一位gr=circshift(gr',1)';%位置调整圆周右移一位sig1=ifft(fft(hr).*fft(sig1));%低频sig2=ifft(fft(gr).*fft(sig2));%高频sig=sig1+sig2;%源信号figure(7);%信号重构图subplot(2,1,1)plot(real(sig1));title('低频信号重构');subplot(2,1,2)plot(real(sig2));title('高频信号重构');figure(8);subplot(2,1,1)stem(abs(fft(sig1)));title('重构的低频信号频谱');subplot(2,1,2)stem(abs(fft(sig2)));title('重构的高频信号频谱');%%5.重构信号与原信号比较figure(9)%图形比较plot(real(sig),'r','linewidth',2);holdon;plot(y);legend('重构信号','原始信号')title('重构信号与原始信号的比较图')010203040506070-2-1.5-1-0.500.511.52t=0:1两个正弦信号的混合正弦波图3两个正弦信号的混合正弦波010203040506070051015202530混合信号频谱图4混合信号频谱01020304050607000.511.5低通滤波器图图5低通滤波器图01020304050607000.511.5高通滤波器图图6高通滤波器图010203040506070-2-1012分解信号1010203040506070-2-1012图7分解信号图01020304050607001020304050分解信号1频谱010203040506070010203040分解信号2频谱图8分解信号频谱图010203040506070-2-1012低频信号重构010203040506070-2-1012高频信号重构图9信号重构图0102030405060700102030重构的低频信号频谱0102030405060700102030重构的高频信号频谱图10重构信号频谱010203040506070-2-1.5-1-0.500.511.52重构信号与原始信号的比较图重构信号原始信号图11重构信号与原始信号的比较图（2）%n是阶数n=8，小波函数为db30clc;clear;%%1信号源：混合正弦波f1=50;%频率1f2=100;%频率2n=8;%n是阶数t=0:1/(2^(n)-1):1;N=length(t);y=sin(2*pi*f1*t)+sin(2*pi*f2*t);%正弦波混合figure(1)plot(y);xlabel('t=0:1')title('两个正弦信号的混合正弦波')figure(2)stem(abs(fft(y)));title('混合信号频谱')%%2.小波滤波器的谱分析h=wfilters('db30','l');%低通g=wfilters('db30','h');%高通h=[h,zeros(1,N-length(h))];%补零（圆周卷积，且增大分辨率变于观察）g=[g,zeros(1,N-length(g))];%补零（圆周卷积，且增大分辨率变于观察）figure(3);%滤波器图stem(abs(fft(h)));title('低通滤波器图')figure(4);stem(abs(fft(g)));title('高通滤波器图')%%3.MALLET分解算法(圆周卷积的快速傅里叶变换实现)sig1=ifft(fft(y).*fft(h));%低通(低频分量)sig2=ifft(fft(y).*fft(g));%高通(高频分量)figure(5);%分解信号图subplot(2,1,1)plot(real(sig1));title('分解信号1')subplot(2,1,2)plot(real(sig2));figure(6);%分解信号频谱图subplot(2,1,1)stem(abs(fft(sig1)));title('分解信号1频谱')subplot(2,1,2)stem(abs(fft(sig2)));title('分解信号2频谱')%%4.MALLET重构算法sig1=dyaddown(sig1);%2抽取sig2=dyaddown(sig2);%2抽取sig1=dyadup(sig1);%2插值sig2=dy