小波学习笔记——MATLAB

最新推荐文章于 2023-11-01 21:45:31 发布

shangyingying_1

最新推荐文章于 2023-11-01 21:45:31 发布

阅读量1.3k

点赞数 5

分类专栏： MATLAB 图像处理

本文链接：https://blog.csdn.net/qq_42261630/article/details/116267091

版权

图像处理同时被 2 个专栏收录

20 篇文章 34 订阅

订阅专栏

MATLAB

11 篇文章 7 订阅

订阅专栏

使用MATLAB进行小波的学习，学习尺度函数、小波函数以及小波函数构造的方法

1.db3小波器的提取

2.对一维信号进行小波滤波

3.对图片进行小波滤波

4.自己构建dbN小波滤波器

5.用Cascade算法求解db4尺度函数和小波函数

详细过程如下代码所示：所调用的function函数见最后

clc,clear all;close all
load sumsin.mat
[LOD,HID,LOR,HIR]=wfilters('db3');%db3小波滤波器的提取
figure,
subplot(2,2,1),plot(LOD,'-'),title('Decomposition Low pass filter(分解低通滤波器)'),xlabel('n'),ylabel('Dh(n)');
subplot(2,2,2),plot(HID,'-'),title('Decomposition High pass filter(分解高通滤波器)'),xlabel('n'),ylabel('Dg(n)');
subplot(2,2,3),plot(LOR,'-'),title('Recomposition Low pass filter(重构低通滤波器)'),xlabel('n'),ylabel('Rh(n)');
subplot(2,2,4),plot(HIR,'-'),title('Decomposition High pass filter(重构高通滤波器)'),xlabel('n'),ylabel('Rg(n)');
[PHI,PSI,XVAL]=wavefun('db3',8);%db3小波对应的尺度函数和小波函数的提取，PHI表示尺度函数，PSI表示小波函数，8表示迭代次数
figure,
subplot(1,2,1),plot(PHI),title('尺度函数PHI'),xlabel('t'),ylabel('\phi(t)');
subplot(1,2,2),plot(PSI,'O'),title('小波函数PSI'),xlabel('t'),ylabel('\psi(t)');
%----------用小波滤波器处理一维信号sumsin---------------%
X_L=filter(LOD,3,sumsin);%其中，LOD为加项滤波器，3为减项滤波器，该值设什么数影响不大
X_H=filter(HID,3,sumsin);
figure,
subplot(3,1,1),plot(sumsin),title('原始一维信号'),xlabel('t'),ylabel('X(t)');
subplot(3,1,2),plot(X_L),title('db3低通滤波器滤波信号'),xlabel('t'),ylabel('X_L(t)');
subplot(3,1,3),plot(X_H),title('db3高通滤波器滤波信号'),xlabel('t'),ylabel('X_H(t)');

[LOD1,HID1,LOR1,HIR1]=wfilters('haar');%haar小波滤波器的提取
figure,
subplot(2,2,1),plot(LOD1,'-'),title('Decomposition Low pass filter(分解低通滤波器)'),xlabel('n'),ylabel('Dh(n)');
subplot(2,2,2),plot(HID1,'-'),title('Decomposition High pass filter(分解高通滤波器)'),xlabel('n'),ylabel('Dg(n)');
subplot(2,2,3),plot(LOR1,'-'),title('Recomposition Low pass filter(重构低通滤波器)'),xlabel('n'),ylabel('Rh(n)');
subplot(2,2,4),plot(HIR1,'-'),title('Decomposition High pass filter(重构高通滤波器)'),xlabel('n'),ylabel('Rg(n)');
[PHI1,PSI1,XVAL1]=wavefun('haar',8);%db3小波对应的尺度函数和小波函数的提取，PHI表示尺度函数，PSI表示小波函数，8表示迭代次数
figure,
subplot(1,2,1),plot(PHI1),title('尺度函数PHI'),xlabel('t'),ylabel('\phi(t)');
subplot(1,2,2),plot(PSI1,'O'),title('小波函数PSI'),xlabel('t'),ylabel('\psi(t)');
X_L1=filter(LOD1,3,sumsin);%其中，LOD为加项滤波器，3为减项滤波器，该值设什么数影响不大
X_H1=filter(HID1,3,sumsin);
figure,
subplot(3,1,1),plot(sumsin),title('原始一维信号'),xlabel('t'),ylabel('X(t)');
subplot(3,1,2),plot(X_L1),title('db3低通滤波器滤波信号'),xlabel('t'),ylabel('X_L(t)');
subplot(3,1,3),plot(X_H1),title('db3高通滤波器滤波信号'),xlabel('t'),ylabel('X_H(t)');


%---------对detail_1图片的滤波--------%
load detail_1
Y=255./(max(max(X))-min(min(X))).*(X-min(min(X)).*ones(size(X)));
figure,
subplot(2,2,1),imshow(X,[]),title('原始图像');
subplot(2,2,2),imshow(Y,[]),title('灰度拉伸后的图像');
H_average=fspecial('average');%二维均值算子获得
X_average=imfilter(X,H_average);
H_sobel=fspecial('sobel');%sobel横向边缘凸显滤波器获得
X_sobel=imfilter(X,H_sobel);
H_prewitt=fspecial('prewitt');%prewitt横向边缘凸显滤波器获得
X_prewitt=imfilter(X,H_prewitt);
H_haar=[0.5 -0.5;-0.5 0.5];
X_haar=imfilter(X,H_haar);
figure,
subplot(2,2,1),imshow(X_average,[]),title('均值滤波结果');
subplot(2,2,2),imshow(X_sobel,[]),title('sobel滤波结果');
subplot(2,2,3),imshow(X_prewitt,[]),title('prewitt滤波结果');
subplot(2,2,4),imshow(X_haar,[]),title('haar滤波结果');


%--------自己构建dbN滤波器在dbfilter函数中--------------------%
hn=dbfilter(4);%获得db滤波器
X_dbN=imfilter(X,hn);
figure,imshow(X_dbN,[]),title('db4滤波结果');


%-----------用Cascade算法求解db4尺度函数和小波函数-----------------%
%调用系数矩阵的函数，进行初值计算
h=[0.2304 0.7148 0.6309 -0.0280 -0.1870 0.0308 0.0329 -0.0106 ];%输入滤波器
A=Filter_Matrix(h);
[v,d]=eig(A);
phi_n=v(:,1)./sum(v(:,1));%对应特征值的特征向量归一化得到尺度函数整数点的值
[x11,phi11,psi11]=Cascade_f(phi_n,h,3);
figure,plot(x11,phi11),title('Cascade算法—尺度函数'),xlabel('x'),ylabel('\phi(x)')
figure,plot(x11,psi11),title('Cascade算法—小波函数'),xlabel('x'),ylabel('\psi(x)')

主函数中调用的function函数如下：

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 A = Filter_Matrix(h)
%用Cascade算法求解尺度函数和小波函数
%子函数：尺度函数初值的整数点的计算
%输入变量h为滤波器
%输出变量A为整数点的滤波器组成的稀疏矩阵
n=length(h);
c=sqrt(2).*h;
A=zeros(n-2,n-2);
for i=1:n-2
    A_begin=max(1,2*i-n+1);%A矩阵非零元素位置从1开始，或者从2*i-n+1开始
    A_end=min(2*i,n-2);%A矩阵非零元素到2*i结束，或者到n-2结束
    %滤波器填充为倒序
    c_begin=min(2*i,n);
    c_end=max(1,2*i-n+3);
    A(i,A_begin:1:A_end)=c(c_begin:-1:c_end);
end

end

function [x,s_phi,s_psi] = Cascade_f(phi_n,h,M)
%求特征值和特征向量，提取出对应特征值为1的特征向量，并归一化
%子函数：分整数点的计算函数
%输入整数点值，滤波器，分点层数，输出尺度函数s_phi，小波函数s_psi的离散点值
%滤波器长度和尺度函数初始向量长度
N=length(h);
PL=length(phi_n);
%求解phi(x)加细点值
%离散x 值
x=0:1/2^M:PL+1;
nx=length(x);
%为phi(x)赋初值0
s_phi=zeros(1,nx);
s_psi=zeros(1,nx);

%由低通滤波器计算高通滤波器
A=-1.*ones(size(h));
B=cumprod(A);%累计积实现了(-1)^n的操作
g=B.*fliplr(h);%fliplr命令使低通滤波器反序，即实现h(-n)操作

%在整数点赋初值phi(k),k=1,2,...,N-1
s_phi(2^M+1:2^M:nx-1)=phi_n(:);
%计算小波函数整数点的值
for k=1:PL
    nx_1=k*2^M+1;
    for n=1:N
        x_2=2*k-n+1;
        if(x_2>0 & x_2<PL+1)
            s_psi(nx_1)=s_psi(nx_1)+sqrt(2).*g(n)*phi_n(x_2);
        end
    end
end
%迭代求解分整数点值
for m=1:M
    for k=1/2^m:1/2^(m-1):PL+1
        nx_1=find(x==k);
        for n=1:N
            x_2=2*k-n+1;
            if(x_2>0&x_2<PL+1)
                nx_2=find(x==x_2);
                s_phi(nx_1)=s_phi(nx_1)+sqrt(2).*h(n)*s_phi(nx_2);
                s_psi(nx_1)=s_psi(nx_1)+sqrt(2).*g(n)*s_phi(nx_2);
            end
        end
    end
end

            







end