实验目的:通过编程实现离散快速小波变换Mallat 算法,从而加深理解二维
小波变换的分解与合成,同时,提高编程能力和matlab 的应用,为以后的学习打下基础。
实验原理:
1、Mallat 快速算法
本实验使用离散快速小波变换快速算法Mallat 算法,算法原理如下
1(2)j j k n n
c h n k c -=-∑ (1)
1(2)j j k n n
d g n k c -=-∑ (2)
重构算法: 1
(2)(2)j j j n
k k n
n
c h n k c g n k
d -=-+-∑∑ (3)
对于(1)、(2)等效于1
j n c -经过冲击响应为[]h n -和[]g n -的数字滤波器,然后再分别进行“二抽取”,Mallat 分解算法的滤波器表示形式如下图
C j-1
d j (k)
C j (k)
用滤波器表示如下图
d j
C j C j-1(k)
2、 255*255
10lg
PSNR MSE
=
'2
11
()*M
N
ij
ij i j f
f MSE M N
==-=
∑∑
{}ij f '{}ij f 分别表示原始图像和重建后的图像,1,1i M j N ≤≤≤≤。
3、边界延拓方法有零延拓、周期延拓、对称周期延拓、常数连续延拓等,本实验采用以上四种方法进行原图像的1/8延拓,并进行重构,各种延拓方法所对应的函数为yan0(x)、yancir
(x )、yan(x)、yanc(x),在主程序中,需要某种延拓,便调用某种函数。
实验编程思路:
为使程序易于理解,在不考虑算法复杂度的情况下,分解程序采用简洁的循环计算出下一级的分解系数,程序采用的编程思想如下
[][][]11100[0][1][2][3][4][5]
001[1]00[0][1][2][3]00[1][2][3][4][5]00[0][1]12j j j j j j c c h h h h h h c c h h h h n c n h h h h h h c ---??
???????????????
???=???
???????????--?????????????
?
以上矩阵等式左面是进行二抽样的结果,[0][1]2
j j n
c c -是j 分解的低频部分。同理,对
于j 分解的高频部分有如下矩阵形式:
[][][]11
100[0][1][2][3][4][5]0
01[1]00[0][1][2][3]0
0[1][2][3][4][5]00[0][1]12j j j j j d d g g g g g g d d g g g g n d n g g g g g g d ---????
????????????????=?
??
????
???????--?????
?????????
分解程序: lenx=size(x,2);%x 为一维向量
lenh=size(h,2); h=[h,zeros(1,(lenx-lenh))]; g=[g,zeros(1,(lenx-lenh))];
r1(1)=sum(h.*x); r2(1)=sum(g.*x);
for k=1:1:(lenx/2-1) %循环求出下一级低频和高频分量
h=[h(end-1:end),h(1:(end-2))];
r1(k+1)=sum(h.*x); g=[g(end-1:end),g(1:1:(end-2))];
r2(k+1)=sum(g.*x);
end y=[r1,r2];
对于重构算法,其等效形式为
[][][]1(2)(2)j j j n
n
c n h n k c k g n k
d k -=-+-∑∑
上式等号右边部分实质上是对变量k 的数字卷积运算,程序采用频域相乘代替卷积,重建程序为 y=ifft(fft(c3,lenx).*fft(h,lenx))+ ifft(fft(d3,lenx).*fft(g,lenx));
实验结果及分析:
1、多尺度分解与重构图像
二维小波变换采用小波采用db3,其峰值信噪比PNSR=,并对三级分解图像进行归一化,求出0的个数为37626,其所占的百分比为%。
2、延拓重建图像
延拓方法周期延拓对称周期延拓零延拓常数连续延拓PNSR
从PNSR结果可知,在各种延拓中,对称周期延拓的重建图像结果最好,相比之下零延拓图像效果不如其他方法延拓。
3、不同小波下重构图像的性质
用不同小波进行图像重构,所得的重构图像能量分布如下
用各种小波进行重构后的图像的均值方差如下表。 小波 db1 db2 db3 db4 均值 方差 2272
2272
2272
2272
附录:
1、主函数程序
clc;clear;
X=imread('');%路径
X=double(X);
% S=yancir(X);
A=mallatdec2(X,'db3', 3);
image(abs(A));
colormap(gray(255));
title('3级多尺度分解图像');
Y=mallatrec2(A,'db3',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('重建图像');
zerosn=numberzeros(A);
% Y=Y(33:288,33:288); %当调用延拓图像时,从延拓的重建图像进行截取csize=size(X);
sr=csize(1);
sc=csize(2);
mse=sum(sum( (Y-X).^2,1))/(sr*sc);
psnr=10*log(255*255/mse)/log(10)
2、分解程序
function Y=mallatdec2(X,wname,level)
%输入:X 载入的二维图像像数值;
% level 小波分解次(级)数设定值(如果设定值超过最高可分解次数,按最高分解次数分% wname 小波名字wavelet name
%输出:Y 多极小波分解后的小波系数矩阵
[h,g]=wfilters(wname,'d'); %h,g 分别为低通和高通滤波器
X=double(X);
t=1;
[hh,ll]=size(X);
while t<=level
%先进行行小波变换
for row=1:hh
Y(row,1:ll)=mdec1(X(row,1:ll),h,g) ;
end
%再进行列小波变换
for col=1:ll
temp=mdec1( Y(1:hh,col)',h,g);
Y(1:hh,col)=temp';
end
t=t+1;
hh=hh/2;
ll=ll/2;
X=Y;
end
function y=mdec1(x,h,g)
%输入:x 行数组
% h 为低通滤波器
% g 为高通滤波器
%输出: y 进行一级小波分解后的系数
lenx=size(x,2);
lenh=size(h,2);
h=[h,zeros(1,(lenx-lenh))];
g=[g,zeros(1,(lenx-lenh))];
r1(1)=sum(h.*x);
r2(1)=sum(g.*x);
for k=1:1:(lenx/2-1)
h=[h(end-1:end),h(1:(end-2))];
r1(k+1)=sum(h.*x);
g=[g(end-1:end),g(1:1:(end-2))];
r2(k+1)=sum(g.*x);
end
y=[r1,r2];
3、重建程序
function Y=mallatrec2(X,wname,level)
%输入:X 载入的小波系数矩阵;
% level 小波分解次(级)数设定值(如果设定值超过最高可分解次数,按最高分解次数% wname 小波名字wavelet name
%输出:Y 重构图像矩阵
[h,g]=wfilters(wname,'d'); %h,g 分别为重构低通滤波器和重构高通滤波器
hz=size(X,2);
h1=hz/(2^(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;
function y=mrec1(x,h,g)
%输入:x 行数组
% h 为低通滤波器
% g 为高通滤波器
%输出: y 进行一级小波重构后值
lenx=size(x,2);
lenh=size(h,2);
h=[h,zeros(1,(lenx-lenh))];
g=[g,zeros(1,(lenx-lenh))];;
c3=dyadup(x(1,1:lenx*,0,lenx); %内插零use para 0
d3=dyadup(x(1,(lenx*+1):lenx),0,lenx); %use para 0
y=ifft(fft(c3,lenx).*fft(h,lenx))+ ifft(fft(d3,lenx).*fft(g,lenx));
4、延拓程序
%零延拓程序
function T=yan0(x)
[r0,c0]=size(x);
for s=1:r0
a(s,:)=[zeros(1,c0/8),x(s,:),zeros(1,c0/8)];
end
c0=c0+c0/4;
for t=1:c0
T(:,t)=[zeros(r0/8,1)',a(:,t)',zeros(r0/8,1)']';
End
%连续常数延拓程序
function T=yanc(x)
[r0,c0]=size(x);
for s=1:r0
r1(1:c0/8)=x(s,1);
r2(1:c0/8)=x(s,c0);
a(s,:)=[r1,x(s,:),r2];
end
c0=c0+c0/4;
for t=1:c0
c1(1:r0/8)=a(1,t);
c2(1:r0/8)=a(r0,t);
T(:,t)=[c1,a(:,t)',c2]';
end
%对称周期延拓程序
function T=yancir(x)
[r0,c0]=size(x);
for s=1:r0
a(s,:)=[x(s,end-(c0/8-1):1:end),x(s,:),x(s,1:1:c0/8)];
end
c0=c0+c0/4;
for t=1:c0
T(:,t)=[a(end-(r0/8-1):1:end,t)',a(:,t)',a(1:1:r0/8,t)']'; end
%周期延拓程序
function T=yan(x)
[r0,c0]=size(x);
for s=1:r0
a(s,:)=[x(s,c0/8:-1:1),x(s,:),x(s,end:-1:end-(c0/8-1))];
end
c0=c0+c0/4;
for t=1:c0
T(:,t)=[a(r0/8:-1:1,t)',a(:,t)',a(end:-1:end-(r0/8-1),t)']'; end
E RDAS实验报告 图像融合实验 数据来源 采用Erdas中examples文件内的2000年Atlanta多光谱TM数据和高清全色Pan数据。两图为同一地区不同坐标影像,故使用前需预处理从而得到实验区域。 目的 多光谱TM数据分辨率较低但包含多波段色彩,而全色Pan数据只包含一层高清影像,为了得到研究区域的高清彩色影像,我们将TM和Pan数据在Erdas2014中进行融合以达到实验目的。 方法 在遥感领域运用较多的融合方法有主成分变换法、比值变换法、小波变换法和HIS变换法。本实验则运用HIS变换法。IHS属于色度空间变换,从多光谱彩色合成影像上分离出代表信息的明度(I)和代表光谱信息的色调(H)、饱和度(S)等3个分量,并采用相同区域的高分辨率全色波段数据代替明度(I)进行空间信息融合。 步骤 1.几何校正 因原始图像空间坐标不同,需选取控制点进行几何校正。本实验校正方法为多项式法,选取6个控制点进行校正,其校正叠加截图如下:
2.叠加剪切 由校正结果可知两图像只有部分区域重合,所以建立AOI对重合区域进行剪切,以得到研究区域,截图如下: 3.重采样 因多光谱图像分辨率较低,像元点较大,若要与全色图融合出高清影像需进行重采样来调整像元大小,以达到与高清图一致。 4.二次剪切 因图为栅格,统一像元后,边缘区必然会有一定的扩展(如下图),虽说扩展的范围较小,但在科研应用方面不符合要求,故须二次剪切。 5.RGB转HIS
TM图像选取前三层再分别赋予蓝、绿、红三色,转化为HIS格式,如下图: 6.直方图匹配 将高清图像直方图以标准图像的直方图为标准作变换,使全色光图和HIS图中I层两图像的直方图相同和近似,从而使两幅图像具有类似的色调和反差,以便作进一步的运算。 7.图像叠加 运用Layer stack功能将全色光高清图和H、S图层进行叠加即所谓的图像融合。它将多波段图层组合到了一起,从而得到新的包含多个有助于研究者使用的多波段影像。 8.IHS转RGB
小波变换与傅里叶变换的对比、异同 一、基的概念 两者都是基,信号都可以分成无穷多个他们的和(叠加)。而展开系数就是基与信号之间的内积,更通俗的说是投影。展开系数大的,说明信号和基是足够相似的。这也就是相似性检测的思想。但我们必须明确的是,傅里叶是0-2pi 标准正交基,而小波是-inf到inf之间的基。因此,小波在实轴上是紧的。而傅里叶的基(正弦或余弦),与此相反。而小波能不能成为Reisz基,或标准稳定的正交基,还有其它的限制条件。此外,两者相似的还有就是PARSEVAL定理。(时频能量守恒)。 二、离散化的处理 傅里叶变换,是一种数学的精妙描述。但计算机实现,却是一步步把时域和频域离散化而来的。第一步,时域离散化,我们得到离散时间傅里叶变换(DTFT),频谱被周期化;第二步,再将频域离散化,我们得到离散周期傅里叶级数(DFS),时域进一步被周期化。第三步,考虑到周期离散化的时域和频域,我们只取一个周期研究,也就是众所周知的离散傅里叶变换(DFT)。这里说一句,DFT是没有物理意义的,它只是我们研究的需要。借此,计算机的处理才成为可能。所有满足容许性条件(从-INF到+INF积分为零)的函数,都可以成为小波。小波作为尺度膨胀和空间移位的一组函数也就诞生了。但连续取值的尺度因子和平移因子,在时域计算量和频域的混叠来说,都是极为不便的。用更为专业的俗语,叫再生核。也就是,对于任何一个尺度a和平移因子b的小波,和原信号内积,所得到的小波系数,都可以表示成,在a,b附近生成的小波,投影后小波系数的线性组合。这就叫冗余性。这时的连续小波是与正交基毫无关系的东西,它顶多也只能作为一种积分变换或基。但它的显微镜特点和相似性检测能力,已经显现出来了。为了进一步更好的将连续小波变换离散化,以下步骤是一种有效方法。第一步,尺度离散化。一般只将a二进离散化,此时b是任意的。这样小波被称为二进小波。第二步,离散b。怎么离散化呢?b取多少才合适呢?于是,叫小波采样定理的东西,就这样诞生了。也就是小波平移的最小距离(采样间隔),应该大于二倍小波基的最高频率(好像类似,记不清了)。所以b取尺度的整数倍就行了。也就是越胖的小波,对应频谱越窄,平移量应该越大,采样间隔越大。当然,第一二两步的频域理解,即在满足频域窗口中心是3倍的频域窗口半径的前提下,频域就在统计上是完美二分的。(但很多小波满足不了这个条件,而且频域窗口能量不?,所以只是近似二分的).这时的小波变换,称为离散二进小波变换.第三步,引入稳定性条件.也就是经过变换后信号能量和原信号能量有什么不等式关系.满足稳定性条件?后,也就是一个小波框架产生了可能.他是数值稳定性的保证.一个稍弱的稳定条件???,就是? MATLAB小波变换指令及其功能介绍 1 一维小波变换的 Matlab 实现 (1) dwt函数 功能:一维离散小波变换 格式:[cA,cD]=dwt(X,'wname') [cA,cD]=dwt(X,Lo_D,Hi_D)别可以实现一维、二维和 N 维DFT 说明:[cA,cD]=dwt(X,'wname') 使用指定的小波基函数 'wname' 对信号X 进行分解,cA、cD 分别为近似分量和细节分量; [cA,cD]=dwt(X,Lo_D,Hi_D) 使用指定的滤波器组 Lo_D、Hi_D 对信 号进行分解。 (2) idwt 函数 功能:一维离散小波反变换 格式:X=idwt(cA,cD,'wname') X=idwt(cA,cD,Lo_R,Hi_R) X=idwt(cA,cD,'wname',L)函数 fft、fft2 和 fftn 分 X=idwt(cA,cD,Lo_R,Hi_R,L) 说明:X=idwt(cA,cD,'wname') 由近似分量 cA 和细节分量 cD 经 小波反变换重构原始信号 X 。 'wname' 为所选的小波函数 X=idwt(cA,cD,Lo_R,Hi_R) 用指定的重构滤波器 Lo_R 和 Hi_R 经小波反变换重构原始信号 X 。 X=idwt(cA,cD,'wname',L) 和 X=idwt(cA,cD,Lo_R,Hi_R,L) 指定返回信号 X 中心附近的 L 个点。 2 二维小波变换的 Matlab 实现 二维小波变换的函数别可以实现一维、二维和 N 维 DFT 函数名函数功能 --------------------------------------------------- dwt2 二维离散小波变换 wavedec2 二维信号的多层小波分解 idwt2 二维离散小波反变换 waverec2 二维信号的多层小波重构 wrcoef2 由多层小波分解重构某一层的分解信号 upcoef2 由多层小波分解重构近似分量或细节分量 detcoef2 提取二维信号小波分解的细节分量 appcoef2 提取二维信号小波分解的近似分量 upwlev2 二维小波分解的单层重构 dwtpet2 二维周期小波变换 idwtper2 二维周期小波反变换 ----------------------------------------------------------- (1) wcodemat 函数 功能:对数据矩阵进行伪彩色编码函数 fft、fft2 和 fftn 分 格式:Y=wcodemat(X,NB,OPT,ABSOL) Y=wcodemat(X,NB,OPT) Y=wcodemat(X,NB) Y=wcodemat(X) 说明:Y=wcodemat(X,NB,OPT,ABSOL) 返回数据矩阵 X 的编码矩阵 Y ;NB 伪编码的最大值,即编码范围为 0~NB,缺省值 NB=16; OPT 指定了编码的方式(缺省值为 'mat'),即:别可以实现 一维、二维和 N 维 DFT OPT='row' ,按行编码 OPT='col' ,按列编码 南京师范大学物理科学与技术学院 医用电子学论文 论文名称:基于小波变换的心电信号噪声消除 院系:物科院 专业:电路与系统 姓名:聂梦雅 学号: 121002043 指导教师:徐寅林 摘要 以小波变换的多分辨率分析为基础, 通过对体表心电信号(ECG) 及其噪声的分析, 对ECG信号中存在的基线漂移、工频干扰及肌电干扰等几种噪声, 设计了不同的小波消噪算法; 并利用MIT/BIH 国际标准数据库中的ECG 信号和程序模拟所产生的ECG 信号, 分别对算法进行了仿真与实验验证。结果表明, 算法能有效地滤除ECG 信号检测中串入的几类主要噪声, 失真度很小, 可满足临床分析与诊断对ECG 波形的要求。 关键词: ECG 信号, 小波变换, 基线漂移, 工频干扰, 肌电干扰 Abstract We apply the multi-resolution analysis (MRA ) of wavelet transform ( WT ) , which was proposed by Mallat [ 5 ] , to suppress the three main types of noises existing in electrocardiogram ( ECG ) signals : baseline wander, power line interference and electro my ographical interference. We apply Mallat algorithm [ 4 ] to suppress the baseline wander in ECG signals. We apply the sof t-thresholding algorithm, proposed by donohoetal on the basis of MRA of WT , to suppress power line interference in ECG signals. We apply Mallat algorithm and then the algorithm proposed by Donohoetal to suppress the electro my ographical interference in ECG signals ,who sefrequency range varies f rom 5Hz to 2kHz. We performed simulations ,using both ECG signals from MIT/BIH database, and ECG signals generated via computer simulation .The results show that the algorithm can suppress the main no isesexisting in ECG signals efficiently with very little distortion, and can satisfy the requirement s of clinical analysis and diagnosis on ECG waveforms. Key words: ECG (electro cardio gram ) signal, wavelet transform , baseline wander, power line interference , electro my ographical interference MATLAB 小波变换指令及其功能介绍 1 一维小波变换的 Matlab 实现 (1 dwt函数 功能:一维离散小波变换 格式:[cA,cD]=dwt(X,'wname' [cA,cD]=dwt(X,Lo_D,Hi_D别可以实现一维、二维和 N 维 DFT 说明:[cA,cD]=dwt(X,'wname' 使用指定的小波基函数 'wname' 对信号X 进行分解,cA 、cD 分别为近似分量和细节分量; [cA,cD]=dwt(X,Lo_D,Hi_D 使用指定的滤波器组 Lo_D、Hi_D 对信号进行分解。 (2 idwt 函数 功能:一维离散小波反变换 格式:X=idwt(cA,cD,'wname' X=idwt(cA,cD,Lo_R,Hi_R X=idwt(cA,cD,'wname',L函数 fft、fft2 和 fftn 分 X=idwt(cA,cD,Lo_R,Hi_R,L 说明:X=idwt(cA,cD,'wname' 由近似分量 cA 和细节分量 cD 经小波反变换重构原始信号 X 。 'wname' 为所选的小波函数 X=idwt(cA,cD,Lo_R,Hi_R 用指定的重构滤波器 Lo_R 和 Hi_R 经小波反变换重构原始信号 X 。 X=idwt(cA,cD,'wname',L 和 X=idwt(cA,cD,Lo_R,Hi_R,L 指定返回信号 X 中心附近的 L 个点。 2 二维小波变换的 Matlab 实现 二维小波变换的函数别可以实现一维、二维和 N 维 DFT 函数名函数功能 --------------------------------------------------- dwt2 二维离散小波变换 wavedec2 二维信号的多层小波分解 idwt2 二维离散小波反变换 waverec2 二维信号的多层小波重构 wrcoef2 由多层小波分解重构某一层的分解信号 upcoef2 由多层小波分解重构近似分量或细节分量 detcoef2 提取二维信号小波分解的细节分量 appcoef2 提取二维信号小波分解的近似分量 upwlev2 二维小波分解的单层重构 dwtpet2 二维周期小波变换 idwtper2 二维周期小波反变换 ----------------------------------------------------------- (1 wcodemat 函数 功能:对数据矩阵进行伪彩色编码函数 fft、fft2 和 fftn 分格式: Y=wcodemat(X,NB,OPT,ABSOL Y=wcodemat(X,NB,OPT Y=wcodemat(X,NB 《医学图像处理》实验报告 实验十:小波变换 日期: 2014年05月06日 摘要 本次实验的实验目的及主要内容是: 一维小波变换和反变换 二维小波变换和反变换 二维小波细节置零、去噪 一、技术讨论 1.1实验原理 小波变换的原理:是指一组衰减震动的波形,其振幅正负相间变化为零,是具有一定的带宽和中心频率波组。小波变换是用伸缩和平移小波形成的小波基来分解(变换)或重构(反变换)时变信号的过程。不同的小波具有不同带宽和中心频率,同一小波集中的带宽与中心频率的比是不变的,小波变换是一系列的带通滤波响应。它的数学过程与傅立叶分析是相似的,只是在傅立叶分析中的基函数是单频的调和函数,而小波分析中的基函数是小波,是一可变带宽内调和函数的组合。 小波去噪的原理:利用小波变换把含噪信号分解到多尺度中,小波变换多采用二进型,然后在每一尺度下把属于噪声的小波系数去除,保留并增强属于信号的小波系数,最后重构出小波消噪后的信号。其中关键是用什么准则来去除属于噪声的小波系数,增强属于信号的部分。 1.2实验方法 1)dwt函数(实现1-D离散小波变换) [cA,cD]=dwt(X,’wname’)使用指定的小波基函数‘wname’对信号X进行分解,cA和cD分别是近似分量和细节分量; [cA,cD]=dwt(X,Lo_D,Hi_D)用指定的滤波器组Lo_D,Hi_D对信号进行分解 2)idwt函数(实现1-D离散小波反变换) X=idwt(cA,cD,’wname’) X=idwt(cA,cD,Lo_R,Hi_R) X=idwt(cA,cD,’wname’,L) X=idwt(cA,cD,Lo_R,Hi_R,L) 由近似分量cA和细节分量cD经过小波反变换,选择某小波函数或滤波器组,L为信号X中心附近的几个点 3)dwt2函数(实现2-D离散小波变换) [cA,cH,cV,cD]=dwt2(X,’wname’) [cA,cH,cV,cD]=dwt2(X,’wname’) cA近似分量,cH水平细节分量,cV垂直细节分量,cD对角细节分量 4)idwt2函数(实现2-D离散反小波变换) X=idwt2(cA,cH,cV,cD,’wname’) X=idwt2(cA,cH,cV,cD,Lo_R,Hi_R) X=idwt2(cA,cH,cV,cD,’wname’,S) X=idwt2(cA,cH,cV,cD,Lo_R,Hi_R,S)MATLAB小波变换指令及其功能介绍(超级有用)解读
基于小波信号的噪声消除matlab实验报告
MATLAB小波变换指令及其功能介绍(超级有用).
小波变换
第9章小波变换基础