第七章+频域处理+.ppt

1、第七章 频域处理,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-1）可用复数表示为,完成这种变换，一般采用的方法是线性正交变换。,(7-1),(7-2),

2、7.2 傅 立 叶 变 换,7.2.1 连续函数的傅立叶变换若把一个一维输入信号作一维傅立叶变换，该信号就被变换到频域上的一个信号，即得到了构成该输入信号的频谱，频谱反映了该输入信号由哪些频率构成。这是一种分析与处理一维信号的重要手段。当一个一维信号f(x)满足狄里赫莱条件，即f(x)（1） 具有有限个间断点； （2） 具有有限个极值点； （3） 绝对可积。,则其傅立叶变换对（傅立叶变换和逆变换）一定存在。在实际应用中，这些条件一般总是可以满足的。一维傅立叶变换对的定义为,式中： ，x称为时域变量，u称为频域变量。,以上一维傅立叶变换可以很容易地推广到二维，如果二维函数f(x, y)满足狄里赫

3、莱条件，则它的二维傅立叶变换对为,式中：x, y为时域变量；u, v为频域变量。,7.2.2 离散傅立叶变换要在数字图像处理中应用傅立叶变换， 还需要解决两个问题：一是在数学中进行傅立叶变换的f(x)为连续（模拟）信号， 而计算机处理的是数字信号（图像数据）；二是数学上采用无穷大概念，而计算机只能进行有限次计算。通常， 将受这种限制的傅立叶变换称为离散傅立叶变换（Discrete Fourier Transform，DFT)。设f(x)|f(0), f(1), f(2), , f(N-1)为一维信号f(x)的N个抽样， 其离散傅立叶变换对为,(7-7),(7-8),式中：x，u=0， 1， 2

4、， , 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决定了每个傅立叶变换结果的频率。,通常傅立叶变换为复数形式， 即,(7-11),式中，R(u)和I(u)分别是F(u)的实部和虚部。式(7-11)也可

5、表示成指数形式： F(u)=|F(u) |ej(u),(7-12),其中,(7-13),(7-14),通常称|F(u) |为f(x)的频谱或傅立叶幅度谱，(u)为f(x)的相位谱。频谱的平方称为能量谱或功率谱，它表示为,(7-15),考虑到两个变量，就很容易将一维离散傅立叶变换推广到二维。二维离散傅立叶变换对定义为,(7-16),(7-17),式中：u, x=0, 1, 2, , M-1；v, y=0, 1, 2, , N-1；x, y为时域变量，u, v为频域变量。像一维离散傅立叶变换一样，系数1/MN可以在正变换或逆变换中，也可以在正变换和逆变换前分别乘以系数 ，只要两式系数的乘积等于1M

6、N即可。二维离散函数的傅立叶频谱、 相位谱和能量谱分别为,(7-18),(7-19),(7-20),式中，R(u, v)和I(u, v)分别是F(u, v)的实部和虚部。,7.2.3 离散傅立叶变换的性质,表7-1 二维离散傅立叶变换的性质,1. 可分离性由可分离性可知，一个二维傅立叶变换可分解为两步进行， 其中每一步都是一个一维傅立叶变换。先对f(x, y)按行进行傅立叶变换得到F(x, v)，再对F(x, v)按列进行傅立叶变换，便可得到f(x, y)的傅立叶变换结果，如图7-4所示。显然对f(x, y)先按列进行离散傅立叶变换， 再按行进行离散傅立叶变换也是可行的。,图7-4 用两次一维

7、DFT计算二维DFT,2. 平移性质平移性质表明，只要将f(x, y)乘以因子(1)x+y，再进行离散傅立叶变换，即可将图像的频谱原点（0，0）移动到图像中心（M2, N2）处。图7-5是简单方块图像平移的结果。,图7-5 傅立叶频谱平移示意图 (a) 原图像；（b）无平移的傅立叶频谱；（c）平移后的傅立叶频谱,3. 旋转不变性由旋转不变性可知，如果时域中离散函数旋转0角度，则在变换域中该离散傅立叶变换函数也将旋转同样的角度。离散傅立叶变换的旋转不变性如图7-6所示。,图7-6 离散傅立叶变换的旋转不变性 （a） 原始图像； （b） 原始图像的傅立叶频谱； （c） 旋转45后的图像； （d）

8、图像旋转后的傅立叶频谱,7.2.4 快速离散傅立叶变换离散傅立叶变换计算量非常大，运算时间长。可以证明其运算次数正比于N2，特别是当N较大时，其运算时间将迅速增长， 以至于无法容忍。为此，研究离散傅立叶变换的快速算法（Fast Fourier Transform， FFT）是非常有必要的。下面介绍一种称为逐次加倍法的快速傅立叶变换算法（FFT），它是1965年Cooley和Tukey首先提出的。采用该FFT算法，其运算次数正比于NlbN，当N很大时计算量可以大大减少。例如，FFT的运算次数和DFT的运算次数之比，当N=1024时，比值为1102.4；当N=4096时，比值可达1341.3。,由

9、于二维离散傅立叶变换具有可分离性， 即它可由两次一维离散傅立叶变换计算得到，因此，仅研究一维离散傅立叶变换的快速算法即可。先将式（7-7）写成,(7-21),式中，W=e-j2N ，称为旋转因子。,这样，可将式(7-21)所示的一维离散傅立叶变换（DFT）用矩阵的形式表示为,式中，由Wux构成的矩阵称为W阵或系数矩阵。,(7-22),观察DFT的W阵，并结合W的定义表达式W=e-j2N，可以发现系数W是以N为周期的。这样，W阵中很多系数就是相同的， 不必进行多次重复计算，且由于W的对称性，即,因此可进一步减少计算工作量。例如，对于N=4， W阵为,(7-23),由W的周期性得：W4W0，W6W

10、2，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),

11、于是式(7-28)变为,(7-31),进一步考虑W的对称性和周期性可知 和 ， 于是,(7-32),由此，可将一个N点的离散傅立叶变换分解成两个N2短序列的离散傅立叶变换，即分解为偶数和奇数序列的离散傅立叶变换Fe(u)和Fo(u) 。,在此，以计算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)

12、和Fo(1)所构成的结构为蝶形运算单元, 其表示的运算为,(7-34),图7-7 蝶形运算单元,由于Fe(u)和Fo(u)都是4点的DFT，因此，如果对它们再按照奇偶进行分组， 则有,（7-35a）,（7-35b）,图7-8 4点DFT分解为2点DFT的蝶形流程图,图7-9 8点DFT的蝶形流程图,图7-10 8点DFT逐级分解框图,表7-2 自然顺序与码位倒序（N=8）,上述FFT是将f(x)序列按x的奇偶进行分组计算的，称之为时间抽选FFT。如果将频域序列的F(u)按u的奇偶进行分组计算， 也可实现快速傅立叶计算， 这称为频率抽选FFT。至此，读者应该对傅立叶变换的理论基础及其实现方式有所

13、了解。对于计算机专业的学生而言，每个人都应该尝试编写快速傅立叶变换的程序。当然，有关傅立叶变换的算法还有很多， 网上的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）

14、(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)的每

15、一列进行一维变换得到F(y, u)，再沿F(y, u)每一行取一维变换得到F(u, v)，其最终结果是一样的。该结论对反变换核也适用。,7.3.2 图像变换的矩阵表示数字图像都是实数矩阵， 设f(x, y)为MN的图像灰度矩阵， 通常为了分析、推导方便，可将可分离变换写成矩阵的形式： F=PfQ F=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变换之外，还采用许多其他的正交变换

16、。例如：离散余弦变换、沃尔什-哈达玛变换、K-L变换等，下面将对常用的变换作一简要介绍。,7.4 离散余弦变换（DCT）,离散余弦变换（Discrete Cosine Transform， DCT）的变换核为余弦函数。DCT除了具有一般的正交变换性质外， 它的变换阵的基向量能很好地描述人类语音信号和图像信号的相关特征。因此，在对语音信号、图像信号的变换中，DCT变换被认为是一种准最佳变换。近年颁布的一系列视频压缩编码的国际标准建议中，都把DCT作为其中的一个基本处理模块。除此之外， DCT还是一种可分离的变换。,7.4.1 一维离散余弦变换一维DCT的变换核定义为,式中，x, u=0, 1,

17、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,

18、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,

19、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-1 x=N， N+1, , 2N-1,(7-59),按照一维DCT的定义，fe(x)的DCT为,(7-60),式中，Re表示取复数的实部。,由于 为fe(x)的2N点DFT。因此，在作

20、DCT时，可把长度为N的f(x)的长度延拓为2N点的序列fe(x)，然后对fe(x)作DFT，最后取DFT的实部便可得到DCT的结果。同理对于离散余弦逆变换IDCT，可首先将F(u)延拓为,由式(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)点对应于高频成分（注： 此频谱图中未作频谱中心平移）。 由于DF

21、T和IDFT已有快速算法FFT和IFFT，因此可用它们实现快速DCT和IDCT算法FCT及IFCT。不过，由于FFT及IFFT中要涉及到复数运算， 因此这种FCT及IFCT算法并不是最佳的。,图7-11 DFT和DCT的频谱分布 （a）DFT频谱分布； （b） DCT频谱分布,7.5 离散沃尔什-哈达玛变换（WHT）,7.5.1 一维离散沃尔什-哈达玛变换1. 沃尔什函数沃尔什函数是1923年由美国数学家沃尔什提出的。沃尔什函数系是一个完备正交函数系，其值只能取1和1。从排列次序上可将沃尔什函数分为三种定义方法：一种是按照沃尔什排列来定义（按列率排序）；另一种是按照佩利排列来定义（按自然排序）

22、；第三种是按照哈达玛排列来定义。由于哈达玛排序的沃尔什函数是由2n（n=0，1，2，）阶哈达玛矩阵（Hadamard Matrix）得到的，而哈达玛矩阵的最大优点在于它具有简单的递推关系， 即高阶矩阵可用两个低阶矩阵的克罗内克积求得，因此在此只介绍哈达玛排列定义的沃尔什变换。,N=2n(n=0, 1, 2, )阶哈达玛矩阵每一行的符号变化规律对应于某一个沃尔什函数的符号变化规律，即N=2n(n=0, 1, 2, )阶哈达玛矩阵的每一行对应于一个离散沃尔什函数，哈达玛矩阵与沃尔什函数系不同之处仅仅是行的次序不同。2n阶哈达玛矩阵有如下形式：,(7-64),(7-65),(7-66),哈达玛矩阵的

23、阶数是按N2n(n0, 1, 2, )规律排列的，阶数较高的哈达玛矩阵，可以利用矩阵的克罗内克积运算，由低阶哈达玛矩阵递推得到，即,（7-67）,矩阵的克罗内克积(Kronecker Product)运算用符号记作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阶哈达玛矩阵。

24、,由哈达玛矩阵的特点可知，沃尔什-哈达玛变换的本质上是将离散序列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是一幅数字图像及对其进行二维W

25、HT变换的结果。,图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进行变换而得到最终结果，由此便可实现二维的FW

26、HT。,综上所述，WHT是将一个函数变换成取值为1或1的基本函数构成的级数，用它来逼近数字脉冲信号时要比FFT有利。同时， WHT只需要进行实数运算，存储量比FFT要少得多， 运算速度也快得多。因此，WHT在图像传输、 通信技术和数据压缩中被广泛使用。,7.6 用MatrixC+库实现图像 变换的VC+编程,7.6.1 Matrix简介及其与VC+工程的集成1. Matrix简介MatrixC+数学库是MathTools公司利用Matcom技术开发的一个面向专业从事工程技术和科学计算人员的矩阵运算动态链接库。这个C+库提供了绝大多数的关于矩阵类、矩阵操作、 数值计算等函数的定义，它提供了一个双

27、精度Matrix类型Mm， 它可以定义的矩阵是复数矩阵、实数矩阵、稀疏矩阵甚至n维矩阵。同时还提供了近千个与矩阵运算相关的函数，并对如、 *、 /、 ()等操作符进行了重载。库中所提供的函数涉及线性代数、多项式数学、信号处理、文件输入输出、绘图等各个方面。充分利用这些库函数， 便可完成数字图像变换的各种操作。,2. 在C+工程中集成MatrixC+数学库为了将MatrixC+库集成到C+的工程文件中，需要做如下的操作： （1） 将库文件v4501v.lib库文件添加到C+的工程中； （2） 将库声明头文件matlib.h包含到工程中； （3） 初始化库； （4） 创建矩阵变量； （5） 访问矩

28、阵单元； （6） 调用MatrixC+库函数完成矩阵操作。,为了使读者掌握如何在VC+中使用MatrixC+数学库的方法，下面详细介绍从创建工程到最后调用MatrixC+库函数的全过程。1） 创建一个新的VC+工程HTSS在VC+中，采用默认选项创建一个MDI MFC工程。当然MatrixC+库可以用在任何其他类型的工程中，如Console 、 OWL、ActiveX等。,2） 将库文件添加到工程中为使C+编译器能够将MatrixC+库链接到最后的执行文件中去，必须将MatrixC+库加入到C+的工程文件中。用菜单命令Project/Add to Project/Files将C:matcom4

29、5libv4501v.lib（如果v4501v.lib存在不同的文件夹中，请选择正确的路径）加入到C+工程中。,3） 包含matlib.h头文件在要使用MatrixC+库的.h或.cpp文件的头部， 将matlib.h头文件包含到工程中，即在使用MatrixC+库的.h或.cpp文件头部添加如下代码： #include “matlib.h“同时用菜单命令Tools/Options/Directories设置matlib.h所在的路径，以便编译器能找到相应的头文件。,4）初始化MatrixC+库在调用MatrixC+库函数之前，要用“initM(MATCOM_VERSION)“函数来初始化类库调

30、用，并用“exitM()“函数来结束类库调用， 故还应在.cpp 中加入下列代码： initM(MATCOM_VERSION); /用MatrixC+库函数的调用完成矩阵操作exitM(); 当然可以在一个类的构造函数中添加initM(MATCOM_VERSION);,以完成类库的初始化工作。可在该类的析构函数中添加 exitM(); 以完成结束类库调用的工作。MATCOM_VERSION参数是一个在matlib.h头文件中定义的值为4501的常量（不同版本其值有所差异）， 它保证了动态链接库与matlib.h的匹配。,5） 创建矩阵变量添加了MatrixC+库后，可用关键字Mm来定义中矩阵类

31、型变量， 例如： Mm a , b, x; 这条VC+语句在当前工程文件中创建了a、b、x三个矩阵变量，但这时a、b和x还只是空的矩阵， 它们的矩阵单元还不包含任何的值。,6） 访问矩阵单元矩阵单元的访问包括读和写操作，主要有：通过.r()来访问矩阵的实部，通过.i()来访问矩阵的虚部； 通过BR()函数把一个数据列表赋给矩阵各个单元；通过.addr()或.addi()返回矩阵变量的实部或虚部的内存指针来完成对矩阵单元的访问。,3. 安装与MatrixC+库对应的动态链接库当采用MatrixC+数学库完成矩阵的运算时， 要求系统中有相关的ago4501.dll, v4501v.dll, ope

32、ngl32.dll和glu32.dll四个动态链接库文件。如果系统中安装有Matcom4.5，这些动态链接库文件将会自动安装在windowsystem或windowsystem32下， 如果系统没有安装Matcom4.5， 只需将这四个动态链接库文件拷贝到相应目录下即可。,7.6.2 创建图像数据矩阵为利用矩阵运算完成图像变换， 首先应将其图像数据赋给矩阵变量。设指针pBits指向读入内存中的图像数据首地址， 可用如下代码将图像数据赋给矩阵变量（为简单起见， 设图像为8位灰度图像）。1. 在类定义头文件.h中定义矩阵变量Mm m_matBits;,2. 在类实现文件.cpp中完成对矩阵变量的赋

33、值操作/图像数据总字节数， 其中nWidthBytes为每行图像数据占用的字节数DWORD SizeImage = nWidthBytes * nHeight; /通过创建一个1SizeImage的一维0矩阵， 为矩阵变量分配内存m_matBits = zeros(1, SizeImage); /将矩阵由double型转换为unsigned char型， 以便将图像数据赋给矩阵m_matBits = muint8(m_matBits); /首先利用m_matBits.addr()获得矩阵实部的内存指针， 然后利用memcpy()将图像数据一次性,/赋给矩阵变量memcpy(m_matBits.

34、addr(), pBits, SizeImage); /由于Mm型的矩阵在内存中是按先列后行的顺序存储的， 因此， 应先用MatrixC+的库/函数reshape()将一维矩阵m_matBits变为nHeightnWidthBytes的二维矩阵， 再用rot90()函数/将其逆时针旋转90， 将矩阵变为nWidthBytesnHeight的二维矩阵m_matBits = rot90(reshape(m_matBits, nWidthBytes, nHeight); /若图像宽度与其字节宽度不同， 则将由系统补齐的每行字节数为4的整数倍的各列0删除，以减 小矩阵大小加快处理速度,if(nWidt

35、hBytes != nWidth) m_matBits = m_matBits(c_p, colon(1, 1, nWidth); /将矩阵由unsigned char型转换为double型， 以便进行运算 m_matBits = mdouble(m_matBits);,7.6.3 将矩阵数据赋给图像数据区为使数据符合图像文件格式，还应将处理后的矩阵数据赋给图像数据区，其代码如下：,/检查m_matBits矩阵是否为空 if(isempty(m_matBits) = 1) return FALSE; /将矩阵数据范围限定于0255之间 m_matBits = mabs(m_matBits);

36、Mm L = m_matBits 255.0; Mm NotL = !L; m_matBits = times(m_matBits, NotL) + L * 255.0;,/将矩阵由double型转换为unsigned char型， 以便将图像数据赋给矩阵 m_matBits = muint8(m_matBits); /计算nWidth/4的余数 int nTmp = (int)rem(nWidth, 4); int nPadding; /为矩阵各列补0， 以使其为符合4的整数倍 if(nTmp 0) nPadding = 4 - nTmp; m_matBits = cat(2, m_matB

37、its, repmat(muint8(0), (BR(size(m_matBits, 1), nPadding); /将矩阵顺时针旋转90 m_matBits = rot90(m_matBits, -1); /将图像数据赋绘矩阵 memcpy(pBits, m_matBits.addr(), (nWidthBytes * nHeight)*sizeof(unsigned char);,7.6.4 利用矩阵运算进行图像变换图像数据矩阵建立后，便可利用MatrixC+的库函数通过各种矩阵运算完成图像变换的工作。 1. 离散傅立叶变换MatrixC+的库提供了函数dft()、fft()和fft2()

38、分别用于完成一维DFT、 一维FFT和二维FFT运算，函数ifft()和ifft2()用于完成一维FFT和二维FFT的逆傅立叶变换。在图像处理中常用的函数是fft2()和ifft2()，利用这两个函数便可分别完成图像的正、反傅立叶变换。其核心代码如下：,/调用MatrixC+库函数fft2()完成二维离散傅立叶变换 Mm ff = fft2(m_matBits); /调用MatrixC+库函数fftshift()将频域中心移到矩阵中心 ff = fftshift(ff); /调用MatrixC+库函数mabs()计算频谱 m_matBits = mabs(ff);,2. 离散余弦变换在Matr

39、ixC+库中并未提供离散余弦变换函数，因此需要编写相应的变换函数。由式（7-61）和式（7-63）可知， 利用离散傅立叶正反变换便可完成离散余弦正反变换，实际上Matlab提供的dct（）、 idct（）、 dct2（）和idct2（）函数就是如此完成离散余弦正反变换的。仔细阅读相应的.m文件便可得出其实现变换的过程， 这样就可以编写用以完成一维DCT函数dct()的代码。完成一维离散余弦变换函数dct()后，利用二维DCT的可分离性，便可方便的编制出用于二维离散余弦变换函数dct2()。其核心代码为,b = transpose(dct(transpose(dct(arg1); 其中， arg

40、1为待变换矩阵。同理， 可编制出一维IDCT函数变换idct()和二维IDCT函数idct2()。与图像的离散傅立叶变换类似， 实现图像的离散余弦变换的代码如下： Mm ff = dct2(m_matBits); /调用MatrixC+库函数mabs()计算频谱m_matBits = mabs(ff); ,3. 离散沃尔什-哈达玛变换同样，在MatrixC+库中并未提供离散沃尔什-哈达玛变换的函数，但它提供了用于产生哈达玛矩阵的函数hadamard()， 其语法为H=hadamard(N); 其中的元素为1或1， 注意该矩阵只有当n或n/12或n/20为2的整次幂时才存在。利用矩阵形式的WHT

41、， 可由如下代码完成离散沃尔什-哈达玛变换：,/获得矩阵的行数和列数 Mm mSize = size(m_matBits); /生成哈达玛矩阵 Mm HColum = hadamard(mSize.r(1, 1); Mm HRow = hadamard(mSize.r(1, 2); /进行离散哈达玛变换 Mm ff = HRow * m_matBits * HColum; /调整系数以便以图像方式显示结果 ff = ff * 1000/(mSize.r(1, 1)*mSize.r(1, 2); /调用MatrixC+库函数mabs()计算频谱 m_matBits = mabs(ff);,7.7

42、 小波变换简介,7.7.1 小波变换的理论基础信号分析是为了获得时间和频率之间的相互关系。傅立叶变换提供了有关频率域的信息，但有关时间的局部化信息却基本丢失。与傅立叶变换不同，小波变换是通过缩放母小波（Mother wavelet）的宽度来获得信号的频率特征， 通过平移母小波来获得信号的时间信息。对母小波的缩放和平移操作是为了计算小波系数，这些小波系数反映了小波和局部信号之间的相关程度。,1. 连续小波变换（CWT）像傅立叶分析一样，小波分析就是把一个信号分解为将母小波经过缩放和平移之后的一系列小波，因此小波是小波变换的基函数。小波变换可以理解为用经过缩放和平移的一系列小波函数代替傅立叶变换的

43、正弦波和余弦波进行傅立叶变换的结果。图7-13表示了正弦波和小波的区别，由此可以看出，正弦波从负无穷一直延续到正无穷，正弦波是平滑而且是可预测的， 而小波是一类在有限区间内快速衰减到0的函数，其平均值为0， 小波趋于不规则、不对称。,图7-13 正弦波和小波 （a） 正弦波曲线； (b) 小波曲线,从小波和正弦波的形状可以看出，变化剧烈的信号，用不规则的小波进行分析比用平滑的正弦波更好，即用小波更能描述信号的局部特征。连续小波变换（Continuous Wavelet Transform， CWT）用下式表示：,(7-79),式（7-79）表示小波变换是信号f(x)与被缩放和平移的小波函数()

44、之积在信号存在的整个期间里求和的结果。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计算主要有如下五个步骤： 第一步： 取一个小波， 将其与原始信号的

45、开始一节进行比较。 第二步： 计算数值C， C表示小波与所取一节信号的相似程度，计算结果取决于所选小波的形状， 如图7-16所示。第三步：向右移动小波，重复第一步和第二步，直至覆盖整个信号，如图7-17所示。第四步： 伸展小波， 重复第一步至第三步， 如图7-18所示。,图7-16 计算系数值C,图7-17 计算平移后系数值C,图7-18 计算尺度后系数值C,第五步：对于所有缩放，重复第一步至第四步。小波的缩放因子与信号频率之间的关系是：缩放因子scale越小，表示小波越窄，度量的是信号的细节变化，表示信号频率越高；缩放因子scale越大， 表示小波越宽，度量的是信号的粗糙程度，表示信号频率越

46、低。,2. 离散小波变换（DWT）在每个可能的缩放因子和平移参数下计算小波系数，其计算量相当大， 将产生惊人的数据量，而且有许多数据是无用的。如果缩放因子和平移参数都选择为2j（j0且为整数）的倍数， 即只选择部分缩放因子和平移参数来进行计算， 就会使分析的数据量大大减少。使用这样的缩放因子和平移参数的小波变换称为双尺度小波变换（Dyadic Wavelet Transform），它是离散小波变换（Discrete Wavelet Transform， DWT）的一种形式。通常离散小波变换就是指双尺度小波变换。,执行离散小波变换的有效方法是使用滤波器， 该方法是Mallat于1988年提出的，

47、称为Mallat算法。这种方法实际上是一种信号分解的方法， 在数字信号处理中常称为双通道子带编码。用滤波器执行离散小波变换的概念如图7-19所示。S表示原始的输入信号， 通过两个互补的滤波器组， 其中一个滤波器为低通滤波器，通过该滤波器可得到信号的近似值A（Approximations），另一个为高通滤波器， 通过该滤波器可得到信号的细节值D（Detail）。,图7-19 小波分解示意图,在小波分析中，近似值是大的缩放因子计算的系数，表示信号的低频分量，而细节值是小的缩放因子计算的系数，表示信号的高频分量。实际应用中，信号的低频分量往往是最重要的，而高频分量只起一个修饰的作用。如同一个人的声音

48、一样， 把高频分量去掉后，听起来声音会发生改变，但还能听出说的是什么内容，但如果把低频分量删除后，就会什么内容也听不出来了。,由图7-19可以看出离散小波变换可以表示成由低通滤波器和高通滤波器组成的一棵树。原始信号经过一对互补的滤波器组进行的分解称为一级分解，信号的分解过程也可以不断进行下去，也就是说可以进行多级分解。如果对信号的高频分量不再分解，而对低频分量进行连续分解，就可以得到信号不同分辨率下的低频分量， 这也称为信号的多分辨率分析。如此进行下去， 就会形成图7-20所示的一棵比较大的分解树， 称其为信号的小波分解树（Wavelet Decomposition Tree）。实际中， 分解的级数取决于要分析的信号数据特征及用户的具体需要。,图7-20 多级信号分解示意图 （a） 信号分解； (b) 小波分数； （c）小波分解树,图7-21 小波分解下采样示意图,3. 小波重构将信号的小波分解的分量进行处理后，一般还要根据需要把信号恢复出来，也就是利用信号的小波分解的系数还原出原始信号，这一过程称为小波重构（Wavelet Reconstruction）或叫做小波合成（Wavelet Synthesis）。这一合成过程的数学运算叫做逆离散小波变换（Inverse Discrete Wavelet Transform， IDWT）。,