小波学习笔记（信号分解与重构）——MATLAB

一维信号分解命令：idwt，其实现原理为Mallat算法

%----------Mallat算法和重构算法实现与MATLAB自带函数进行对比---------------%
[Lo_D,Hi_D,Lo_R,Hi_R]=wfilters('db3');
h=Lo_D;
%hn=dbfilter(3);
%A=-1.*ones(1.length(hn));
%B=cumprod(A);%累计积实现了(-1)^n的操作
%gn=(-1).*B.*fliplr(hn);%fliplr命令使低通滤波器反序，即实现h（-n）操作，得到反序滤波器用于信号分解
%h=fliplr(hn);
%g=fliplr(gn);
n=length(h)-1;
X_new=Edge_extend(sumsin,'syms',n);

%-------自编Mallat算法程序----------%
[c,d]=mallat_compostion(X_new,Lo_D,Hi_D);
X_R=mallat_reconstruction(c,d,Lo_R,Hi_R);
dd=max(X_R-sumsin);

%------MATLAB自带小波分解------------%
[c1,d1]=dwt(sumsin,'db3');
X_I=idwt(c1,d1,'db3');

%------两种方法结果对比----------%
max(c-c1);%max(d-d1);max(X_R-X_I)

用到的函数代码

function hn = dbfilter(N)
%求解dbN滤波器
%输入变量N为消失矩，输出变量hn为滤波器


%输入初值dbN的消失矩N，滤波器的长度为2N
a=1;p=1;q=1;
hn=[1,1];
for k=1:N-1
    hn=conv(hn,[1,1]);%卷积得到二项式系数，即描述（1+e^(-iw)）^k的展开系数
    a=-a*0.25*(N-1+k)/k;%计算a（k）
    p=conv(p,[1,-2,1]);%卷积得到三项式系数，即描述（e^(iw)+e^(-iw)-2）^k
    q=[0 q 0]+a*p;%q表示求和运算，计算整个L(w)的系数分布
end
%求出L(w)的多项式的根，并排序
Q_root=sort(roots(q));
%求出L(w)的前N-1个根构成的多项式，real是实部，t的长度为N
t=real(poly(Q_root(1:N-1)));
%完成H(w)=(1+e^(-iw))^N*l(w),hn为前面循环中卷积到N得到的hn,t相当于l(w)
hn=conv(hn,t);
%归一化，保证hn的求和为sqrt(2)
hn=hn*sqrt(2)/sum(hn);



end

function X_new= Edge_extend(X,Edge_X,n)
%子函数：边界延拓函数
%输入参数X为原始信号，Edge_X边界延拓方式，n为延拓的长度
%输出参数X_new为延拓后的信号
%实现零延拓，常数延拓，对称延拓和周期延拓
if strcmp(Edge_X,'zeros')
    B=zeros(1,n);
    X_new=[B,X,B];
elseif strcmp(Edge_X,'const')
    B1=X(1).*ones(1,n);
    B2=X(end).*ones(1,n);
    X_new=[B1,X,B2];
elseif strcmp(Edge_X,'syms')
    B1=fliplr(X(1:n));%反序运算
    B2=fliplr(X(end-n+1:end));
    X_new=[B1,X,B2];
else
    %period extend
    B1=X(end-n+1:end);
    B2=X(1:n);
    X_new=[B1,X,B2];
end
    
end

function [c,d] = mallat_compostion(X_new,h,g)
%Mallat分解算法
%输入参数X_new为边界延拓后的信号，h、g分别为分解用的低通高通滤波器
%输出参数c,d分别对应低通和高通信号
%卷积形式，所用参数valid表示保留没有零延拓信号的部分
c1=conv(X_new,h,'valid');
d1=conv(X_new,g,'valid');
%信号下采样，即抽取原始信号的偶数位置信号
c=dyaddown(c1);
d=dyaddown(d1);
end

function X =mallat_reconstruction(C,D,h,g)
%输入参数C为分解后的低通信号，D为分解后的高通信号，h和g分别为重构用的低通和高通滤波器
%输出参数X为重构信号

C0=dyadup(C);%信号上采样，即将原始信号插零，长度变成一倍
D0=dyadup(D);
X=conv(C0,h,'valid')+conv(D0,g,'valid');%卷积形式，所用参数valid表示保留没有零延拓信号的部分

end