第7章-频域处理--数字图像处理-教学课件.ppt

第七章 频域处理 第七章 频域处理 7.1 频域世界与频域变换频域世界与频域变换 7.2 傅立叶变换傅立叶变换 7.3 频域变换的一般表达式频域变换的一般表达式 7.4 离散余弦变换离散余弦变换 7.5 离散沃尔什哈达玛变换离散沃尔什哈达玛变换 7.6 用用MatrixC+库实现图像变换的库实现图像变换的VC+编程编程7.7 小波变换简介小波变换简介 第七章 频域处理 7.1 频域世界与频域变换频域世界与频域变换图7-1任意波形可分解为正弦波的加权和第七章 频域处理 图7-2正弦波的振幅A和相位第七章 频域处理 图7-3图7-1（a）波形的频域表示(a)幅频特性；(b)相频特性第七章 频域处理 第七章 频域处理 7.2 傅傅 立立 叶叶 变变 换换 7.2.1 连续函数的傅立叶变换连续函数的傅立叶变换若把一个一维输入信号作一维傅立叶变换，该信号就被变换到频域上的一个信号，即得到了构成该输入信号的频谱，频谱反映了该输入信号由哪些频率构成。这是一种分析与处理一维信号的重要手段。当一个一维信号f(x)满足狄里赫莱条件，即f(x)（1）具有有限个间断点；（2）具有有限个极值点；（3）绝对可积。第七章 频域处理 则其傅立叶变换对（傅立叶变换和逆变换）一定存在。在实际应用中，这些条件一般总是可以满足的。一维傅立叶变换对的定义为(7-3)(7-4)式中：，x称为时域变量，u称为频域变量。第七章 频域处理 以上一维傅立叶变换可以很容易地推广到二维，如果二维函数f(x,y)满足狄里赫莱条件，则它的二维傅立叶变换对为(7-5)(7-6)式中：x,y为时域变量；u,v为频域变量。第七章 频域处理 第七章 频域处理(7-7)(7-8)式中：x，u=0，1，2，,N1。注：式（7-8）中的系数1/N也可以放在式（7-7）中，有时也可在傅立叶正变换和逆变换前分别乘以，这是无关紧要的，只要正变换和逆变换前系数乘积等于1/N即可。第七章 频域处理 由欧拉公式可知(7-9)将式（7-9）代入式（7-7），并利用cos()=cos()，可得(7-10)可见，离散序列的傅立叶变换仍是一个离散的序列，每一个u对应的傅立叶变换结果是所有输入序列f(x)的加权和（每一个f(x)都乘以不同频率的正弦和余弦值），u决定了每个傅立叶变换结果的频率。第七章 频域处理 第七章 频域处理 通常称|F(u)|为f(x)的频谱或傅立叶幅度谱，(u)为f(x)的相位谱。频谱的平方称为能量谱或功率谱，它表示为(7-15)考虑到两个变量，就很容易将一维离散傅立叶变换推广到二维。二维离散傅立叶变换对定义为(7-16)(7-17)第七章 频域处理 第七章 频域处理 7.2.3 离散傅立叶变换的性质离散傅立叶变换的性质 表表7-1 二维离散傅立叶变换的性质二维离散傅立叶变换的性质 第七章 频域处理 第七章 频域处理 第七章 频域处理 2.平移性质平移性质平移性质表明，只要将f(x,y)乘以因子(1)x+y，再进行离散傅立叶变换，即可将图像的频谱原点（0，0）移动到图像中心（M2,N2）处。图7-5是简单方块图像平移的结果。图7-5傅立叶频谱平移示意图(a)原图像；（b）无平移的傅立叶频谱；（c）平移后的傅立叶频谱(a)(b)(c)第七章 频域处理 3.旋转不变性旋转不变性由旋转不变性可知，如果时域中离散函数旋转0角度，则在变换域中该离散傅立叶变换函数也将旋转同样的角度。离散傅立叶变换的旋转不变性如图7-6所示。图7-6离散傅立叶变换的旋转不变性（a）原始图像；（b）原始图像的傅立叶频谱；（c）旋转45后的图像；（d）图像旋转后的傅立叶频谱(a)(b)(d)(c)第七章 频域处理 7.2.4 快速离散傅立叶变换快速离散傅立叶变换离散傅立叶变换计算量非常大，运算时间长。可以证明其运算次数正比于N2，特别是当N较大时，其运算时间将迅速增长，以至于无法容忍。为此，研究离散傅立叶变换的快速算法（FastFourierTransform，FFT）是非常有必要的。下面介绍一种称为逐次加倍法的快速傅立叶变换算法（FFT），它是1965年Cooley和Tukey首先提出的。采用该FFT算法，其运算次数正比于NlbN，当N很大时计算量可以大大减少。例如，FFT的运算次数和DFT的运算次数之比，当N=1024时，比值为1102.4；当N=4096时，比值可达1341.3。第七章 频域处理 由于二维离散傅立叶变换具有可分离性，即它可由两次一维离散傅立叶变换计算得到，因此，仅研究一维离散傅立叶变换的快速算法即可。先将式（7-7）写成(7-21)式中，W=e-j2N，称为旋转因子。第七章 频域处理 这样，可将式(7-21)所示的一维离散傅立叶变换（DFT）用矩阵的形式表示为式中，由Wux构成的矩阵称为W阵或系数矩阵。(7-22)第七章 频域处理 第七章 频域处理 由W的周期性得：W4W0，W6W2，W9W1；再由W的对称性可得：W3W1，W2W0。于是式(7-23)可变为(7-24)第七章 频域处理 可见N=4的W阵中只需计算W0和W1两个系数即可。这说明W阵的系数有许多计算工作是重复的，如果把一个离散序列分解成若干短序列，并充分利用旋转因子W的周期性和对称性来计算离散傅立叶变换，便可以简化运算过程，这就是FFT的基本思想。设N为2的正整数次幂，即如令M为正整数，且N=2M(7-25)(7-26)第七章 频域处理 将式（7-26）代入式（7-21），离散傅立叶变换可改写成如下形式：由旋转因子W的定义可知，因此式(7-27)变为现定义(7-27)(7-28)(7-29)(7-30)第七章 频域处理 第七章 频域处理 在此，以计算N=8的DFT为例，此时n=3，M=4。由式(7-31)和式(7-32)可得(7-33)第七章 频域处理 式（7-33）中，u取07时的F(u)、Fe(u)和Fo(u)的关系可用图7-7描述。左方的两个节点为输入节点，代表输入数值；右方两个节点为输出节点，表示输入数值的叠加，运算由左向右进行。线旁的W18和W18为加权系数，定义由F(1)、F(5)、Fe(1)和Fo(1)所构成的结构为蝶形运算单元,其表示的运算为(7-34)第七章 频域处理 第七章 频域处理 由于Fe(u)和Fo(u)都是4点的DFT，因此，如果对它们再按照奇偶进行分组，则有（7-35a）（7-35b）第七章 频域处理 图7-84点DFT分解为2点DFT的蝶形流程图第七章 频域处理 图7-98点DFT的蝶形流程图第七章 频域处理 图7-108点DFT逐级分解框图第七章 频域处理 表表7-2 自然顺序与码位倒序（自然顺序与码位倒序（N=8）第七章 频域处理 上述FFT是将f(x)序列按x的奇偶进行分组计算的，称之为时间抽选FFT。如果将频域序列的F(u)按u的奇偶进行分组计算，也可实现快速傅立叶计算，这称为频率抽选FFT。至此，读者应该对傅立叶变换的理论基础及其实现方式有所了解。对于计算机专业的学生而言，每个人都应该尝试编写快速傅立叶变换的程序。当然，有关傅立叶变换的算法还有很多，网上的FFT算法源代码也非常多，但不建议大家拿来就用。当你得到类似的代码后，一定要认真分析其实现过程和思路，只有这样才能不断地提高编程水平。第七章 频域处理 7.3 频域变换的一般表达式频域变换的一般表达式 7.3.1 可分离变换可分离变换二维傅立叶变换可用通用的关系式来表示：(7-36)(7-37)式中：x,u=0,1,2,M1；y,v=0,1,2,N1；g(x,y,u,v)和h(x,y,u,v)分别称为正向变换核和反向变换核。第七章 频域处理 如果g(x,y,u,v)=g1(x,u)g2（y,v）(7-38)h(x,y,u,v)=h1(x,u)h2(y,v)(7-39)则称正、反变换核是可分离的。进一步，如果g1和g2，h1和h2在函数形式上一样，则称该变换核是对称的。二维傅立叶变换对是式（7-36）和式（7-37）的一个特殊情况，它们的核为第七章 频域处理 可见，它们都是可分离的和对称的。如前所述，二维傅立叶变换可以利用变换核的可分离性，用两次一维变换来实现，即可先对f(x,y)的每一行进行一维变换得到F(x,v)，再沿F(x,v)每一列取一维变换得到变换结果F(u,v)。对于其他的图像变换，只要其变换核是可分离的，同样也可用两次一维变换来实现。如果先对f(x,y)的每一列进行一维变换得到F(y,u)，再沿F(y,u)每一行取一维变换得到F(u,v)，其最终结果是一样的。该结论对反变换核也适用。第七章 频域处理 7.3.2 图像变换的矩阵表示图像变换的矩阵表示数字图像都是实数矩阵，设f(x,y)为MN的图像灰度矩阵，通常为了分析、推导方便，可将可分离变换写成矩阵的形式：F=PfQF=P-1FQ-1其中，F、f是二维MN的矩阵；P是MM矩阵；Q是NN矩阵。(7-44)(7-43)(7-42)式中，u=0,1,2,M1，v=0,1,2,N1。第七章 频域处理 对二维离散傅立叶变换，则有(7-45)(7-46)实践中，除了DFT变换之外，还采用许多其他的正交变换。例如：离散余弦变换、沃尔什-哈达玛变换、K-L变换等，下面将对常用的变换作一简要介绍。第七章 频域处理 7.4 离散余弦变换（离散余弦变换（DCT）离散余弦变换（DiscreteCosineTransform，DCT）的变换核为余弦函数。DCT除了具有一般的正交变换性质外，它的变换阵的基向量能很好地描述人类语音信号和图像信号的相关特征。因此，在对语音信号、图像信号的变换中，DCT变换被认为是一种准最佳变换。近年颁布的一系列视频压缩编码的国际标准建议中，都把DCT作为其中的一个基本处理模块。除此之外，DCT还是一种可分离的变换。第七章 频域处理 7.4.1 一维离散余弦变换一维离散余弦变换一维DCT的变换核定义为式中，x,u=0,1,2,N1；(7-47)(7-48)一维DCT定义如下：设f(x)|x=0,1,N-1为离散的信号列。(7-49)式中，u,x=0,1,2,N1。第七章 频域处理 将变换式展开整理后，可以写成矩阵的形式，即F=Gf（7-50）其中(7-51)第七章 频域处理 一维DCT的逆变换IDCT定义为(7-52)式中，x,u=0,1,2,N1。可见一维DCT的逆变换核与正变换核是相同的。第七章 频域处理 7.4.2 二维离散余弦变换二维离散余弦变换考虑到两个变量，很容易将一维DCT的定义推广到二维DCT。其正变换核为(7-53)式中，C(u)和C(v)的定义同式（7-48）；x,u=0,1,2,M1；y,v=0,1,2,N1。二维DCT定义如下：设f(x,y)为MN的数字图像矩阵，则(7-54)第七章 频域处理 式中：x,u=0,1,2,M1；y,v=0,1,2,N1。二维DCT逆变换定义如下：(7-55)式中：x,u=0,1,2,M1；y,v=0,1,2,N1。类似一维矩阵形式的DCT，可以写出二维DCT的矩阵形式如下：F=GfGT(7-56)第七章 频域处理 同时，由式(7-55)和式(7-54)可知二维DCT的逆变换核与正变换核相同，且是可分离的，即(7-57)式中：C(u)和C(v)的定义同式（7-48）；x,u=0,1,2,M1；y,v=0,1,2,N1。第七章 频域处理 通常根据可分离性，二维DCT可用两次一维DCT来完成，其算法流程与DFT类似，即(7-58)第七章 频域处理 7.4.3 快速离散余弦变换快速离散余弦变换离散余弦变换的计算量相当大，在实用中非常不方便，也需要研究相应的快速算法。目前已有多种快速DCT（FCT），在此介绍一种由FFT的思路发展起来的FCT。首先，将f(x)延拓为x=0,1,2,N-1x=N，N+1,2N-1(7-59)按照一维DCT的定义，fe(x)的DCT为(7-60)第七章 频域处理 式中，Re表示取复数的实部。第七章 频域处理 由于为fe(x)的2N点DFT。因此，在作DCT时，可把长度为N的f(x)的长度延拓为2N点的序列fe(x)，然后对fe(x)作DFT，最后取DFT的实部便可得到DCT的结果。同理对于离散余弦逆变换IDCT，可首先将F(u)延拓为u=0,1,2,N-1u=N,N+1,2N-1(7-62)第七章 频域处理 由式(7-52)可得，DCT的IDCT为（7-63）由式（7-63）可见，IDCT可由的2N点的IDFT来实现。第七章 频域处理 最后要注意的是二维DCT的频谱分布，其谱域分布与DFT相差一倍，如图7-11所示。从图中可以看出，对于DCT而言，(0,0)点对应于频谱的低频成分，(N-1,N-1)点对应于高频成分，而同阶的DFT中，(N2,N2)点对应于高频成分（注：此频谱图中未作频谱中心平移）。由于DFT和IDFT已有快速算法FFT和IFFT，因此可用它们实现快速DCT和IDCT算法FCT及IFCT。不过，由于FFT及IFFT中要涉及到复数运算，因此这种FCT及IFCT算法并不是最佳的。第七章 频域处理 图7-11DFT和DCT的频谱分布（a）DFT频谱分布；（b）DCT频谱分布第七章 频域处理 7.5 离散沃尔什离散沃尔什-哈达玛变换（哈达玛变换（WHT）7.5.1 一维离散沃尔什一维离散沃尔什-哈达玛变换哈达玛变换 1.沃尔什函数沃尔什函数沃尔什函数是1923年由美国数学家沃尔什提出的。沃尔什函数系是一个完备正交函数系，其值只能取1和1。从排列次序上可将沃尔什函数分为三种定义方法：一种是按照沃尔什排列来定义（按列率排序）；另一种是按照佩利排列来定义（按自然排序）；第三种是按照哈达玛排列来定义。由于哈达玛排序的沃尔什函数是由2n（n=0，1，2，）阶哈达玛矩阵（HadamardMatrix）得到的，而哈达玛矩阵的最大优点在于它具有简单的递推关系，即高阶矩阵可用两个低阶矩阵的克罗内克积求得，因此在此只介绍哈达玛排列定义的沃尔什变换。第七章 频域处理 N=2n(n=0,1,2,)阶哈达玛矩阵每一行的符号变化规律对应于某一个沃尔什函数的符号变化规律，即N=2n(n=0,1,2,)阶哈达玛矩阵的每一行对应于一个离散沃尔什函数，哈达玛矩阵与沃尔什函数系不同之处仅仅是行的次序不同。2n阶哈达玛矩阵有如下形式：(7-64)(7-65)(7-66)第七章 频域处理 哈达玛矩阵的阶数是按N2n(n0,1,2,)规律排列的，阶数较高的哈达玛矩阵，可以利用矩阵的克罗内克积运算，由低阶哈达玛矩阵递推得到，即（7-67）第七章 频域处理 矩阵的克罗内克积(KroneckerProduct)运算用符号记作AB，其运算规律如下：设第七章 频域处理 则(7-68)(7-69)第七章 频域处理 2.离散沃尔什离散沃尔什-哈达玛变换哈达玛变换一维离散沃尔什变换定义为(7-70)一维离散沃尔什逆变换定义为(7-71)式中，Walsh(u,x)为沃尔什函数。若将Walsh(u,x)用哈达玛矩阵表示，并将变换表达式写成矩阵形式，则式（7-70）和式（7-71）分别为第七章 频域处理 和(7-72)(7-73)式中，HN为N阶哈达玛矩阵。第七章 频域处理 由哈达玛矩阵的特点可知，沃尔什-哈达玛变换的本质上是将离散序列f(x)的各项值的符号按一定规律改变后，进行加减运算，因此，它比采用复数运算的DFT和采用余弦运算的DCT要简单得多。第七章 频域处理 7.5.2 二维离散沃尔什变换二维离散沃尔什变换很容易将一维WHT的定义推广到二维WHT。二维WHT的正变换核和逆变换核分别为(7-74)和(7-75)式中：x,u=0,1,2,M1；y,v=0,1,2,N1。第七章 频域处理 二维离散沃尔什变换的矩阵形式表达式为和求这两个信号的二维WHT。第七章 频域处理 根据题意，式（7-76）中的M=N=4，其二维WHT变换核为第七章 频域处理 所以第七章 频域处理 第七章 频域处理 再如，图7-12是一幅数字图像及对其进行二维WHT变换的结果。图7-12二维WHT结果（a）原图像；（b）WHT结果第七章 频域处理 注：图7-12中的结果是按哈达玛变换后再用沃尔什排序的结果。从以上例子可看出，二维WHT具有能量集中的特性，而且原始数据中数字越是均匀分布，经变换后的数据越集中于矩阵的边角上。因此，二维WHT可用于压缩图像信息。第七章 频域处理 7.5.3 快速沃尔什变换（快速沃尔什变换（FWHT）类似于FFT，WHT也有快速算法FWHT，也可将输入序列f(x)按奇偶进行分组，分别进行WHT。FWHT的基本关系为(7-78)WHT的变换核是可分离和对称的，因此二维WHT也可分为两个一维的WHT分别用FWHT进行变换而得到最终结果，由此便可实现二维的FWHT。第七章 频域处理 综上所述，WHT是将一个函数变换成取值为1或1的基本函数构成的级数，用它来逼近数字脉冲信号时要比FFT有利。同时，WHT只需要进行实数运算，存储量比FFT要少得多，运算速度也快得多。因此，WHT在图像传输、通信技术和数据压缩中被广泛使用。第七章 频域处理 7.6 用用MatrixC+库实现图像库实现图像变换的变换的VC+编程编程 7.6.1 Matrix简介及其与简介及其与VC+工程的集成工程的集成1.Matrix简介简介MatrixC+数学库是MathTools公司利用Matcom技术开发的一个面向专业从事工程技术和科学计算人员的矩阵运算动态链接库。这个C+库提供了绝大多数的关于矩阵类、矩阵操作、数值计算等函数的定义，它提供了一个双精度Matrix类型Mm，它可以定义的矩阵是复数矩阵、实数矩阵、稀疏矩阵甚至n维矩阵。同时还提供了近千个与矩阵运算相关的函数，并对如、*、/、()等操作符进行了重载。库中所提供的函数涉及线性代数、多项式数学、信号处理、文件输入输出、绘图等各个方面。充分利用这些库函数，便可完成数字图像变换的各种操作。第七章 频域处理 2.在在C+工程中集成工程中集成MatrixC+数学库数学库为了将MatrixC+库集成到C+的工程文件中，需要做如下的操作：（1）将库文件v4501v.lib库文件添加到C+的工程中；（2）将库声明头文件matlib.h包含到工程中；（3）初始化库；（4）创建矩阵变量；（5）访问矩阵单元；（6）调用MatrixC+库函数完成矩阵操作。第七章 频域处理 为了使读者掌握如何在VC+中使用MatrixC+数学库的方法，下面详细介绍从创建工程到最后调用MatrixC+库函数的全过程。1）创建一个新的VC+工程HTSS在VC+中，采用默认选项创建一个MDIMFC工程。当然MatrixC+库可以用在任何其他类型的工程中，如Console、OWL、ActiveX等。第七章 频域处理 2）将库文件添加到工程中为使C+编译器能够将MatrixC+库链接到最后的执行文件中去，必须将MatrixC+库加入到C+的工程文件中。用 菜 单 命 令 Project/Add to Project/Files将 C:matcom45libv4501v.lib（如果v4501v.lib存在不同的文件夹中，请选择正确的路径）加入到C+工程中。第七章 频域处理 3）包含matlib.h头文件在要使用MatrixC+库的.h或.cpp文件的头部，将matlib.h头文件包含到工程中，即在使用MatrixC+库的.h或.cpp文件头部添加如下代码：#includematlib.h同时用菜单命令Tools/Options/Directories设置matlib.h所在的路径，以便编译器能找到相应的头文件。第七章 频域处理 4）初始化MatrixC+库在 调 用 MatrixC+库 函 数 之 前，要 用initM(MATCOM_VERSION)函数来初始化类库调用，并用exitM()函数来结束类库调用，故还应在.cpp中加入下列代码：initM(MATCOM_VERSION);/用MatrixC+库函数的调用完成矩阵操作exitM();当然可以在一个类的构造函数中添加initM(MATCOM_VERSION);第七章 频域处理 以完成类库的初始化工作。可在该类的析构函数中添加exitM();以完成结束类库调用的工作。MATCOM_VERSION参数是一个在matlib.h头文件中定义的值为4501的常量（不同版本其值有所差异），它保证了动态链接库与matlib.h的匹配。第七章 频域处理 5）创建矩阵变量添加了MatrixC+库后，可用关键字Mm来定义中矩阵类型变量，例如：Mma,b,x;这条VC+语句在当前工程文件中创建了a、b、x三个矩阵变量，但这时a、b和x还只是空的矩阵，它们的矩阵单元还不包含任何的值。第七章 频域处理 6）访问矩阵单元矩阵单元的访问包括读和写操作，主要有：通过.r()来访问矩阵的实部，通过.i()来访问矩阵的虚部；通过BR()函数把一个数据列表赋给矩阵各个单元；通过.addr()或.addi()返回矩阵变量的实部或虚部的内存指针来完成对矩阵单元的访问。第七章 频域处理 3.安装与安装与MatrixC+库对应的动态链接库库对应的动态链接库当采用MatrixC+数学库完成矩阵的运算时，要求系统中有相关的ago4501.dll,v4501v.dll,opengl32.dll和glu32.dll四个动态链接库文件。如果系统中安装有Matcom4.5，这些动态链接库文件将会自动安装在windowsystem或windowsystem32下，如果系统没有安装Matcom4.5，只需将这四个动态链接库文件拷贝到相应目录下即可。第七章 频域处理 7.6.2 创建图像数据矩阵创建图像数据矩阵为利用矩阵运算完成图像变换，首先应将其图像数据赋给矩阵变量。设指针pBits指向读入内存中的图像数据首地址，可用如下代码将图像数据赋给矩阵变量（为简单起见，设图像为8位灰度图像）。1.在类定义头文件在类定义头文件.h中定义矩阵变量中定义矩阵变量Mmm_matBits;第七章 频域处理 2.在类实现文件在类实现文件.cpp中完成对矩阵变量的赋值操作中完成对矩阵变量的赋值操作/图像数据总字节数，其中nWidthBytes为每行图像数据占用的字节数DWORDSizeImage=nWidthBytes*nHeight;/通过创建一个1SizeImage的一维0矩阵，为矩阵变量分配内存m_matBits=zeros(1,SizeImage);/将矩阵由double型转换为unsignedchar型，以便将图像数据赋给矩阵m_matBits=muint8(m_matBits);/首先利用m_matBits.addr()获得矩阵实部的内存指针，然后利用memcpy()将图像数据一次性第七章 频域处理/赋给矩阵变量memcpy(m_matBits.addr(),pBits,SizeImage);/由于Mm型的矩阵在内存中是按先列后行的顺序存储的，因此，应先用MatrixC+的库/函数reshape()将一维矩阵m_matBits变为nHeightnWidthBytes的二维矩阵，再用rot90()函数/将其逆时针旋转90，将矩阵变为nWidthBytesnHeight的二维矩阵m_matBits=rot90(reshape(m_matBits,nWidthBytes,nHeight);/若图像宽度与其字节宽度不同，则将由系统补齐的每行字节数为4的整数倍的各列0删除，以减小矩阵大小加快处理速度第七章 频域处理 if(nWidthBytes!=nWidth)m_matBits=m_matBits(c_p,colon(1,1,nWidth);/将矩阵由unsignedchar型转换为double型，以便进行运算m_matBits=mdouble(m_matBits);第七章 频域处理 7.6.3 将矩阵数据赋给图像数据区将矩阵数据赋给图像数据区为使数据符合图像文件格式，还应将处理后的矩阵数据赋给图像数据区，其代码如下：/检查m_matBits矩阵是否为空if(isempty(m_matBits)=1)returnFALSE;/将矩阵数据范围限定于0255之间m_matBits=mabs(m_matBits);MmL=m_matBits255.0;MmNotL=!L;m_matBits=times(m_matBits,NotL)+L*255.0;第七章 频域处理/将矩阵由double型转换为unsignedchar型，以便将图像数据赋给矩阵m_matBits=muint8(m_matBits);/计算nWidth/4的余数intnTmp=(int)rem(nWidth,4);intnPadding;/为矩阵各列补0，以使其为符合4的整数倍if(nTmp0)nPadding=4-nTmp;m_matBits=cat(2,m_matBits,repmat(muint8(0),(BR(size(m_matBits,1),nPadding);/将矩阵顺时针旋转90m_matBits=rot90(m_matBits,-1);/将图像数据赋绘矩阵memcpy(pBits,m_matBits.addr(),(nWidthBytes*nHeight)*sizeof(unsignedchar);第七章 频域处理 7.6.4 利用矩阵运算进行图像变换利用矩阵运算进行图像变换图像数据矩阵建立后，便可利用MatrixC+的库函数通过各种矩阵运算完成图像变换的工作。1.离散傅立叶变换离散傅立叶变换MatrixC+的库提供了函数dft()、fft()和fft2()分别用于完成一维DFT、一维FFT和二维FFT运算，函数ifft()和ifft2()用于完成一维FFT和二维FFT的逆傅立叶变换。在图像处理中常用的函数是fft2()和ifft2()，利用这两个函数便可分别完成图像的正、反傅立叶变换。其核心代码如下：第七章 频域处理/调用MatrixC+库函数fft2()完成二维离散傅立叶变换Mmff=fft2(m_matBits);/调用MatrixC+库函数fftshift()将频域中心移到矩阵中心ff=fftshift(ff);/调用MatrixC+库函数mabs()计算频谱m_matBits=mabs(ff);第七章 频域处理 2.离散余弦变换离散余弦变换在MatrixC+库中并未提供离散余弦变换函数，因此需要编写相应的变换函数。由式（7-61）和式（7-63）可知，利用离散傅立叶正反变换便可完成离散余弦正反变换，实际上Matlab提供的dct（）、idct（）、dct2（）和idct2（）函数就是如此完成离散余弦正反变换的。仔细阅读相应的.m文件便可得出其实现变换的过程，这样就可以编写用以完成一维DCT函数dct()的代码。完成一维离散余弦变换函数dct()后，利用二维DCT的可分离性，便可方便的编制出用于二维离散余弦变换函数dct2()。其核心代码为第七章 频域处理 b=transpose(dct(transpose(dct(arg1);其中，arg1为待变换矩阵。同理，可编制出一维IDCT函数变换idct()和二维IDCT函数idct2()。与图像的离散傅立叶变换类似，实现图像的离散余弦变换的代码如下：Mmff=dct2(m_matBits);/调用MatrixC+库函数mabs()计算频谱m_matBits=mabs(ff);第七章 频域处理 3.离散沃尔什离散沃尔什-哈达玛变换哈达玛变换同样，在MatrixC+库中并未提供离散沃尔什-哈达玛变换的函数，但它提供了用于产生哈达玛矩阵的函数hadamard()，其语法为H=hadamard(N);其中的元素为1或1，注意该矩阵只有当n或n/12或n/20为2的整次幂时才存在。利用矩阵形式的WHT，可由如下代码完成离散沃尔什-哈达玛变换：第七章 频域处理/获得矩阵的行数和列数MmmSize=size(m_matBits);/生成哈达玛矩阵MmHColum=hadamard(mSize.r(1,1);MmHRow=hadamard(mSize.r(1,2);/进行离散哈达玛变换Mmff=HRow*m_matBits*HColum;/调整系数以便以图像方式显示结果ff=ff*1000/(mSize.r(1,1)*mSize.r(1,2);/调用MatrixC+库函数mabs()计算频谱m_matBits=mabs(ff);第七章 频域处理 7.7 小波变换简介小波变换简介 7.7.1 小波变换的理论基础小波变换的理论基础信号分析是为了获得时间和频率之间的相互关系。傅立叶变换提供了有关频率域的信息，但有关时间的局部化信息却基本丢失。与傅立叶变换不同，小波变换是通过缩放母小波（Motherwavelet）的宽度来获得信号的频率特征，通过平移母小波来获得信号的时间信息。对母小波的缩放和平移操作是为了计算小波系数，这些小波系数反映了小波和局部信号之间的相关程度。第七章 频域处理 1.连续小波变换（连续小波变换（CWT）像傅立叶分析一样，小波分析就是把一个信号分解为将母小波经过缩放和平移之后的一系列小波，因此小波是小波变换的基函数。小波变换可以理解为用经过缩放和平移的一系列小波函数代替傅立叶变换的正弦波和余弦波进行傅立叶变换的结果。图7-13表示了正弦波和小波的区别，由此可以看出，正弦波从负无穷一直延续到正无穷，正弦波是平滑而且是可预测的，而小波是一类在有限区间内快速衰减到0的函数，其平均值为0，小波趋于不规则、不对称。第七章 频域处理 图7-13正弦波和小波（a）正弦波曲线；(b)小波曲线第七章 频域处理 从小波和正弦波的形状可以看出，变化剧烈的信号，用不规则的小波进行分析比用平滑的正弦波更好，即用小波更能描述信号的局部特征。连续小波变换（ContinuousWaveletTransform，CWT）用下式表示：(7-79)式（7-79）表示小波变换是信号f(x)与被缩放和平移的小波函数()之积在信号存在的整个期间里求和的结果。CWT的变换结果是许多小波系数C，这些系数是缩放因子（scale）和平移（positon）的函数。第七章 频域处理 基本小波函数()的缩放和平移操作含义如下：(1)缩放。简单地讲，缩放就是压缩或伸展基本小波，缩放系数越小，则小波越窄，如图7-14所示。图7-14小波的缩放操作第七章 频域处理(2)平移。简单地讲，平移就是小波的延迟或超前。在数学上，函数f(t)延迟k的表达式为f(t-k)，如图7-15所示。图7-15小波的平移操作(a)小波函数(t)；(b)位移后的小波函数(t-k)第七章 频域处理 CWT计算主要有如下五个步骤：第一步：取一个小波，将其与原始信号的开始一节进行比较。第二步：计算数值C，C表示小波与所取一节信号的相似程度，计算结果取决于所选小波的形状，如图7-16所示。第三步：向右移动小波，重复第一步和第二步，直至覆盖整个信号，如图7-17所示。第四步：伸展小波，重复第一步至第三步，如图7-18所示。第七章 频域处理 图7-16计算系数值C第七章 频域处理 图7-17计算平移后系数值C 第七章 频域处理 图7-18计算尺度后系数值C第七章 频域处理 第五步：对于所有缩放，重复第一步至第四步。小波的缩放因子与信号频率之间的关系是：缩放因子scale越小，表示小波越窄，度量的是信号的细节变化，表示信号频率越高；缩放因子scale越大，表示小波越宽，度量的是信号的粗糙程度，表示信号频率越低。第七章 频域处理 2.离散小波变换（离散小波变换（DWT）在每个可能的缩放因子和平移参数下计算小波系数，其计算量相当大，将产生惊人的数据量，而且有许多数据是无用的。如果缩放因子和平移参数都选择为2j（j0且为整数）的倍数，即只选择部分缩放因子和平移参数来进行计算，就会使分析的数据量大大减少。使用这样的缩放因子和平移参数的小波变换称为双尺度小波变换（DyadicWaveletTransform），它是离散小波变换（DiscreteWaveletTransform，DWT）的一种形式。通常离散小波变换就是指双尺度小波变换。第七章 频域处理 执行离散小波变换的有效方法是使用滤波器，该方法是Mallat于1988年提出的，称为Mallat算法。这种方法实际上是一种信号分解的方法，在数字信号处理中常称为双通道子带编码。用滤波器执行离散小波变换的概念如图7-19所示。S表示原始的输入信号，通过两个互补的滤波器组，其中一个滤波器为低 通 滤 波 器，通 过 该 滤 波 器 可 得 到 信 号 的 近 似 值A（Approximations），另一个为高通滤波器，通过该滤波器可得到信号的细节值D（Detail）。第七章 频域处理 图7-19小波分解示意图第七章 频域处理 在小波分析中，近似值是大的缩放因子计算的系数，表示信号的低频分量，而细节值是小的缩放因子计算的系数，表示信号的高频分量。实际应用中，信号的低频分量往往是最重要的，而高频分量只起一个修饰的作用。如同一个人的声音一样，把高频分量去掉后，听起来声音会发生改变，但还能听出说的是什么内容，但如果把低频分量删除后，就会什么内容也听不出来了。第七章 频域处理 由图7-19可以看出离散小波变换可以表示成由低通滤波器和高通滤波器组成的一棵树。原始信号经过一对互补的滤波器组进行的分解称为一级分解，信号的分解过程也可以不断进行下去，也就是说可以进行多级分解。如果对信号的高频分量不再分解，而对低频分量进行连续分解，就可以得到信号不同分辨率下的低频分量，这也称为信号的多分辨率分析。如此进行下去，就会形成图7-20所示的一棵比较大的分解树，称其为信号的小波分解树（WaveletDecompositionTree）。实际中，分解的级数取决于要分析的信号数据特征及用户的具体需要。第七章 频域处理 图7-20多级信号分解示意图（a）信号分解；(b)小波分数；（c）小波分解树第七章 频域处理 对于一个信号，如采用图7-19所示的方法，理论上产生的数据量将是原始数据的两倍。于是，根据奈奎斯特（Nyquist）采样定理，可用下采样的方法来减少数据量，即在每个通道内（高通和低通通道）每两个样本数据取一个，便可得到离散小波变换的系数（Coefficient），分别用cA和cD表示，如图7-21所示。图中表示下采样。第七章 频域处理 图7-21小波分解下采样示意图第七章 频域处理 3.小波重构小波重构将信号的小波分解的分量进行处理后，一般还要根据需要把信号恢复出来，也就是利用信号的小波分解的系数还原出原始信号，这一过程称为小波重构（WaveletReconstruction）或叫做小波合成（WaveletSynthesis）。这一合成过程的数学运算叫做逆离散小波变换（Inverse Discrete Wavelet Transform，IDWT）。第七章 频域处理 图7-22小波重构算法示意图第七章 频域处理 1）重构近似信号与细节信号由图7-22可知，由小波分解的近似系数和细节系数可以重构出原始信号。同样，可由近似系数和细节系数分别重构出信号的近似值或细节值，这时只要近似系数或细节系数置为零即可。图7-23是对第一层近似信号或细节信号进行重构的示意图。第七章 频域处理 图7-23重构近似和细节信号示意（a）重构近似信号；(b)重构细节信号第七章 频域处理 2）多层重构在图7-23中，重构出信号的近似值A1与细节值D1之后，则原信号可用A1D1S重构出来。对应于信号的多层小波分解，小波的多层重构如图7-24所示。由图7-24可见重构过程为：A3D3A2；A2D2A1；A1+D1S。信号重构中，滤波器的选择非常重要，关系到能否重构出满意的原始信号。低通分解滤波器（L）和高通分解滤波器（H）及重构滤波器组（L和H）构成一个系统，这个系统称为正交镜像滤波器（QuadratureMirrorFilters，QMF）系统，如图7-25所示。第七章 频域处理 图7-24多层小波重构示意图第七章 频域处理 图7-25多层小波分解和重构示意图第七章 频域处理 4.小波包分析小波包分析小波分析是将信号分