小波分析实验：二维离散小波变换(Mallat快速算法)

小波分析实验：实验2二维离散小波变换（Mallat快速算法）
实验目的：
在理解离散小波变换原理和Mallat快速算法的基础上，通过编程对图像进行二维离散小波变换，从而加深对二维小波分解和重构的理性和感性认识，并能提高编程能力，为今后的学习和工作奠定基础。

实验工具：
计算机，matlab6.5
分解算法:
重构算法: “"二工必（刃- 2上*［十三g （刃- 2k ）d ［ *
分解算法写成矩阵的形式！ （lb g 的长度为4）
4[0]如]力⑵ h[3] 0 0 0 '
［勺【0】• 记"
h[0] h[\]h[2]山⑶ …
• ••••・ • •
C J
=
勺【1] • •
申[2] h[3] 0
0 0
-.^[0] ^[1]_
.勺［乃-1】_
>[0] g[l] g ⑵ g[3] 0 • • •
e=
• 0 •
g[0] g[l]g ⑵ • • g[3]
■ • •・
■ 0
• D J =
<[i]
■
•
目2] ■
g[3]
0 0
…茎0] 畀］
|g[0] g[l] g[2] g[3] 0 0 0 I
0 0 g[0] g[l]g[2] S [3] - 0
• ••••• • •
•・•・・■ • • g[2] g[3] 0 0
0 ...g[0] g[l]J |_勺4-1[
叨］
I
二
・(2»
于是Mallat分解公式为矩阵变换？丄
Cj- = PC^................. ⑶卩
D j = Q D J-L..... .......... ⑷
重构算法写成矩阵变换：-
C J_I =C$ + Dj------------------------------------ (5) 4
M N
PPq. 一片『峰值信噪比计算公式：P沁沁逻竺皿E卢H耿V 屈E M ｛皿，00分别表示原始图像和重建图像，且
本实验采取的一些小技乐P (I)分SW法…
编程时用如下思想：（h, g 的长度为4）“
今
［1］
勺［刀
-1］
■ V
■
■丐⑼
£［1］ 4刀-1】
将数据。

』和低通〔高通）滤波器进行添零到数据长度诳，再2抽祥，相加得
到
（2） 科瞎
小波 系数丄
后一半 系数门
高通.
重构
滤波 器卩
重构结果
加0] 方[1] 农］ 舛3] 0
0 • •
0 o T
=
纽0] 列1] 風2］疏3］
■
町] ■
列2] ■
扯3]
• •
D 0
•
…0
'sm
g[l] 列2] 或习0
• • • 0
0 11 =
0 烈0] 畀］ 02]烈 3] 0
:.°
■
-g[l]
■
g[2] ■
g[3]
• • 0 0
■ 0
•••0 g[0]J[
%[1]
勺4-1]
怕
】
编程吋用如下思想：2
ii[0]山1]爪)
]
川
］
/
〔0
]
0 00
0 00
•••
•••
••
0 00
•托3] 42]
0・・托
3]
• •
• •
• •
0 ••/(!]
0 0 00…灶3] 夙
0] 0 00 0
•••鼻•
•♦• • •
0 0 00 0
■心[0]
「心[1]
■
•
十
'咕[0「
吆[1]
•
弓4-1]_/j-Jn-l]
5・
1[0] 5 di]
川
]
■
••炷
1]
址
0]
勺
[1]
曲
］「
77+
1
实验结果：“
多尺度分堺E3像
SC1）三级小波分解
图Q
50
100
150
200
250
原始图條逅构图像
50 15C 200 250 50 100 130 200 250
EC2 ）原始凰慷和重构图像
峰值信噪比（psm—252.8923db^
附录：
(1)二维小波分解函数
%二维小波分解函数fun cti on Y=mallatdec2(X,w name,level)
%输入:X 载入的二维图像像数值；
% level 小波分解次(级)数设定值(如果设定值超过最高可分解次数,按最高分解次
数分解)
% wname 小波名字wavelet name
%输岀：丫多极小波分解后的小波系数矩阵
[h,g]=wfilters(wname,'d'); %h,g分别为低通和高通滤波器
X=double(X);
t=1;
hh=size(X,2);
while t<=level
%先进行行小波变换
for row=1:hh
Y(row,1:hh)=mdec1(X(row,1:hh),h,g);
end
%再进行列小波变换
for col=1:hh
temp=mdec1( Y(1:hh,col)',h,g);
Y(1:hh,col)=temp';
end
t=t+1;
hh=hh/2;
X=Y;
end
实现%内部子函数，对一行(row)矢量进行一次小波变换，利用fft
fun cti on y=mdec1(x,h,g)
%输入：x行数组
% h为低通滤波器
% g为高通滤波器
%输岀：y进行一级小波分解后的系数
len x=size(x,2);
len h=size(h,2);
rh=h(e nd:-1:1);
rrh=[zeros(1,(le nx-le nh)),rh];
rrh=circshift(rrh',1)';
rg=g(e nd:-1:1);
rrg=[zeros(1,(le nx-le nh)),rg];
rrg=circshift(rrg',1)';
r1=dyaddow n( ifft(fft(x).*fft(rrh,le nx)),1); %use para 1
r2=dyaddow n(ifft(fft(x).*fft(rrg,le nx)),1);
y=[r1,r2];
(2)二维小波重构函数
%二维小波重构函数
fun cti on Y=mallatrec2(X,w name,level)
%输入:X 载入的小波系数矩阵；
,按最高分解次% level 小波分解次(级)数设定值(如果设定值超过最高可分解次数
数分解)
% wname 小波名字wavelet name
%输岀：Y 重构图像矩阵
[h,g]=wfilters(wname,'d'); %h,g分别为重构低通滤波器和重构高通滤波器
hz=size(X,2);
h1=hz/(2A(level-1));
while h1<=hz
%对列变换
for col=1:h1
temp=mrec1(X(1:h1,col)',h,g)';
X(1:h1,col)=temp;
end
%再对行变换
for row=1:h1
temp=mrec1(X(row,1:h1),h,g);
X(row,1:h1)=temp;
end
h1=h1*2;
end
Y=X;
%内部子函数，对一行小波系数进行重构
fun cti on y=mrec1(x,h,g)
%输入：x行数组
% h为低通滤波器
% g为高通滤波器
%输岀：y进行一级小波重构后值
len x=size(x,2);
r4=dyadup(x(1,(le nx*0.5+1):le nx),O); %use para 0 y=ifft(fft(r3,le nx).*fft(h,le nx))+ ifft(fft(r4,le nx).*fft(g,le nx));
(3)测试函数(主函数)
%测试函数(主函数)
r3=dyadup(x(1,1:le nx*0.5),0); %内插零use para 0
clc;clear;
X=imread('E:\Libin 的文档实验2 要求\exp2\LENA.bmp');% 路径X=double(X);
A = mallatdec2(X,'sym2',3);
image(abs(A));
colormap(gray(255));
ti tle(' 多尺度分解图像’)；
Y= mallatrec2(A,'sym2',3);
Y=real(Y);
figure(2);
subplot(1,2,1);
image(X);
colormap(gray(255));
title(' 原始图像');
subplot(1,2,2);
image(Y);
colormap(gray(255));
title(' 重构图像');
csize=size(X);
sr=csize(1);
sc=csize(2);
mse=sum(sum( (Y-X).A2,1))/(sr*sc);
psn r=10*log(255*255/mse)/log(10)
小波分析实验：实验1连续小波变换
实验目的：
在理解连续小波变换原理的基础上，通过编程实现对一维信
号进行连续小波变换，(实验中采用的是墨西哥帽小波)，从而对连续小波变换增加了理性和感性的认识，并能提高编程能力，为今后的学习和工作奠定基础。

实验工具：
计算机，matlab6.5
ZL
实验原理:
•维连续小波变换公式：°
"；（心）二IG
2
J/（O x y y
（:
—）处
-co
1 Q
t 7
=\^ \2
J /⑺*(一^)刃
—8 "
编程时先离散化,在通过求和代替积分。

卩
3+1）心
咲J
k 炮
A /V /（□）%严'"
） 根据
衫 〒
Cf
一丄
£ — b
"Sr)
当小波函第
/(Ox 1^1 V(
本实验采取的一些小技巧：卩
（1）因为实验中的数据通过lo网mta.ma广）（先将数据拷到work文
件下）命令后，数据存在了dM矩阵中，所以兀财就是datik）•由于数
挺样间隔4为0.03 （常量），所以可以把这个系数忽略，实际上本程序冲
-1 ]c-b
用的是WAa.b）=\a\2匸产（切;w{^L）
J T a
y
.Jr
Q）本实验中墨西哥帽函数为：^x/2（x） = （（1~x2）g /|x|-8…
0. otherwise
〔3）如果按照（◎式计算小波系数矩阵，则算法时间复杂度为。

（沪），
3 通的电脑上要运行5分钟左右才能出正确结果。

通过研究（*♦）,可以後
对于同一尺度—屮d）与必匕二I）有如下关系：0（上匕）从 a a a
二项开始与0（匕）相同，只是第一项不同，（在程序中有说明），所以
a
把算法时间复杂度降为。

（泌），程序执行时间在5秒钟内• q
（4）程序的正确性可由c\vt来验证。

亠
（5）本实验还编了个te前程序来说明存在哪些分量，即对正弦波函鐵
ZL 样后进行连续小波变换，分析图象结果。

卍
实验结果：a
(1)实验数据图像和小波系数图像"
50 100 150 200 250 300 350
400
(2) te前数据图像和小波系数图像〜
orignal dat
5
程序附录：
(1)墨西哥帽小波函数,按照(***)式编程
5
fun cti on Y=mexh(x)
if abs(x)<=8
Y=exp(-x*x/2)*(1-x A 2); else
Y=0;
End
(2) 实验程序，按照
(**)式编程，详细过程请参考“本实验采取的一些小技巧” %
clc;clear;
load('data.mat'); len=len gth(dat);
In a=70;
a=zeros(1. In a); wfab=zeros(l na,le n); mexhab=zeros(1,le
% （尺度a）的长度
%小波系数矩阵
%离散化小波系数矩阵
for s=1: Ina %s表示尺度for k=1:le n
mexhab(k)=mexh(k/s);
end
for t=1:len % t表示位移
wfab(s,t)=(sum(mexhab.*dat))/sqrt(s); 替
mexhab=[mexh(-1*t/s),mexhab(1:le n-
1)];
项并右移%将积分用求和代%mexhab（修改
end
figure(1);
plot(dat);
title(' 原始数据图');
figure(2); %小波系数谱image(wfab);
colorm ap(pi nk(128));
title(' 小波系数图');
%surf(wfab);
%title('小波系数谱网格图');
%pwfab=wfab.*wfab; %%瞬态功率谱
%figure(3);
%subplot(1,2,1);
%surf(pwfab);
%title('瞬态功率谱网格图');
%subplot(1,2,2);
%con tour(pwfab);
%title('瞬态功率谱等值线');
(3)test 函数。

%test函数
clc;clear;
for i=1:200
dat(i)=si n(2*pi*i*0.05); % 正弦波函数end
len=len gth(dat);
In a=40;
wfab=zeros(l na,le n);
mexhab=zeros(1,le n);
for s=1: Ina %s表示尺度
for k=1:le n
mexhab(k)=mexh(k/s);
end
for t=1:le n % t表示位移
wfab(s,t)=(sum(mexhab.*dat))/sqrt(s); %将积分用求和代替
mexhab=[mexh(-1*t/s),mexhab(1:le n-1)]; %mexhab(修改第一项并右移
end
end
figure(1);
plot(dat);
title('orig nal dat');
figure(2); %小波系数谱
image(wfab);
colorm ap(pi nk(128));
title(' 正弦波的小波系数图');
(4)用fft 实现cwt
%按照圆周卷积定理，原周卷积和线性卷积的关系L>=M+N-1
%按照圆周卷积的定义，相关和线性卷积的关系(原始算法和线性卷积的关系) 冰意画图理解
clc;clear;
t仁cputime;
load('data.mat');
len=len gth(dat);
a=zeros(1. In a); % a表示尺度
lna=70; % a (尺度)的长度a=zeros(1. In a); % a表示尺度
wfab=zeros(l na,le n); %小波系数矩阵
mexhab=zeros(1,2*le n-1);
data=[zeros(1,le n-1),dat];
Ydata=fft( data ,4*le n);
for s=1:l na
for k=1:2*le n-1
mexhab(k)=mexh((k-le n) /s);
end
temp=ifft( Ydata.*fft( mexhab,4*le n ) ,4*le n);
wfab(s,:)=real(temp(2*le n-1:3*le n-2))/sqrt(s); %
为什么要取实部而不是取模，我也不是很清楚，可是有种感觉
end
figure(1);
plot(dat);
title(' 原始数据图');
figure(2); %小波系数谱
image(wfab);
colorm ap(pi nk(128));
title(' 小波系数谱');
cputime-tl
4)fft快速计算cwt
%按照圆周卷积的定义,
冰意画图理解
clc;clear;
t仁cputime;
load('data.mat');
len=len gth(dat);
Ina=70; % a (尺度)的长度
a=5;
data=[dat,zeros(1,le n)];
Ydata=fft(dat,2*le n);
for s=1:l na
mexhab=zeros(1,2*le n);
k=[-a*s:1:a*s];
mexhab(k+le n)=mexh2(k./s);
temp=ifft( Ydata.*fft( mexhab,2*len ) ,2*le n);
wfab(s,:)=real(temp(len+1:2*len))/sqrt(s); % 呵呵
要取实部而不是取模, end
figure(1);
plot(dat);
title(' 原始数据图');
figure(2); %小波系数谱
image(wfab);
coIorm ap(pi nk(128));
title(' 小波系数谱');
cputime-tl
5)保存为mexh2.m fun cti on Y=mexh2(x) Y=exp(-x.*x/2).*(1-x.A2);。