精华内容
下载资源
问答
  • 多分辨率分析(MRA,Multi-resolution Analysis)又称为尺度分析,是建立在函数空间概念的理论,创建者S.Mallat是在研究图像处理问题时建立这套理论,并提出了著名的Mallat算法。MRA不仅为正交小波基的构造提供了...
  • 另一方面,离散小波变换(DWT)为原始信号的分析和合成提供了足够的信息,同时大大减少了计算时间。 与CWT相比,DWT易于实施。 DWT的基本概念以及它的属性和用于计算它的算法将在本节中介绍。 与前几章一样,提供了...

    1. 为什么需要离散小波变换

    尽管离散连续小波变换可以通过计算机计算连续小波变换,但这并不是真正的离散变换。 实际上,小波序列只是CWT(连续小波变换)的一个采样版本,就信号的重构而言,它提供的信息是高度冗余的。 另一方面,这种冗余需要大量的计算时间和资源。 另一方面,离散小波变换(DWT)为原始信号的分析和合成提供了足够的信息,同时大大减少了计算时间。

    与CWT相比,DWT易于实施。 DWT的基本概念以及它的属性和用于计算它的算法将在本节中介绍。 与前几章一样,提供了一些示例来帮助解释DWT。

    2. 离散小波变换(DWT)

    DWT的基础可以追溯到1976年,当时Croiser,Esteban和Galand设计了一种分解离散时间信号的技术。 同年,Crochiere,Weber和Flanagan在语音信号编码方面做了类似的工作。 他们将分析方案命名为子带编码。 1983年,伯特(Burt)定义了一种非常类似于子带编码的技术,并将其命名为金字塔编码,也称为多分辨率分析。 1989年晚些时候,Vetterli和Le Gall对子带编码方案进行了一些改进,消除了金字塔编码方案中现有的冗余。 下面说明子带编码。 有关离散小波变换和多分辨率分析理论的详细介绍,可以在该主题的许多文章和书籍中找到,这超出了本教程的范围。

    3. 子带编码和多分辨率分析

    主要思想与CWT中的思想相同。 使用数字滤波技术可获得数字信号的time-scale表示。 回想一下,CWT是小波在不同尺度上与信号之间的相关性,该scale(或频率)被用作相似性的度量。 通过更改分析窗口的scale,及时移动窗口,乘以信号以及对所有时间进行积分,可以计算出连续小波变换。 在离散情况下,使用不同截止频率的滤波器来分析不同scale的信号。 信号通过一系列高通滤波器以分析高频,并通过一系列低通滤波器以分析低频。

    信号的分辨率(衡量信号中详细信息的数量)通过滤波操作进行更改,而scale则通过上采样和下采样(子采样)操作进行更改。 对信号进行二次采样对应于降低采样率或删除信号的某些采样。 例如,二次抽样是指将信号的每隔一个样本丢弃一次。 以因子n进行二次采样可将信号中的样本数量减少n次。

    信号的上采样对应于通过向信号添加新的采样来增加信号的采样率。 例如,以2进行上采样是指在信号的每两个采样之间添加一个新的采样,通常为零或内插值。 对信号进行n倍的上采样会使信号中的样本数增加n倍。

    尽管不是唯一的选择,但是DWT系数通常是在二元网格上从CWT采样的,即 s 0 = 2 s_0 = 2 s0=2 t 0 = 1 t_0 = 1 t0=1,产生 s = 2   j s = 2^{\, j} s=2j t = k ⋅ 2   j t = k \cdot 2^{\, j} t=k2j,如第3部分所述。由于信号是离散时间函数,因此在以下讨论中,术语函数和序列将互换使用。 该序列将由 x [ n ] x [n] x[n]表示,其中n是整数。

    该过程开始于使该信号(序列)通过具有脉冲响应 h [ n ] h [n] h[n]的半带数字低通滤波器。 对信号进行滤波对应于信号的卷积与滤波器的脉冲响应的数学运算。 离散时间的卷积运算定义如下:

    x [ n ] ∗ h [ n ] = ∑ k = − ∞ ∞ x [ k ] ⋅ h [ n − k ] x[n] * h[n] = \sum\limits_{k = -\infty}^\infty x[k] \cdot h[n - k] x[n]h[n]=k=x[k]h[nk]···················································(4.1)

    半带低通滤波器可去除所有高于信号最高频率一半的频率。例如,如果信号最大具有1000 Hz分量,则半带低通滤波将去除500 Hz以上的所有频率。

    此时,频率单位特别重要。在离散信号中,频率以弧度表示。因此,就径向频率而言,信号的采样频率等于2p弧度。因此,如果以奈奎斯特速率(是信号中最大频率的两倍)对信号进行采样,则信号中存在的最高频率分量将为p弧度。即,奈奎斯特速率对应于离散频域中的p rad / s。因此,使用Hz不适用于离散信号。但是,每当需要澄清讨论时都使用Hz,因为用Hz来考虑频率是很常见的。应始终记住,离散时间信号的频率单位为弧度。

    使信号通过半带低通滤波器后,根据奈奎斯特规则,可以消除一半采样,因为信号现在的最高频率为 p 2 \frac {p} {2} 2p弧度,而不是p弧度。简单地丢弃所有其他采样将对信号进行二次采样,信号的点数将减半。现在,信号的大小增加了一倍。请注意,低通滤波会删除高频信息,但不会改变音阶。仅子采样过程会更改比例。另一方面,分辨率与信号中的信息量有关,因此,它受滤波操作的影响。半带低通滤波消除了一半的频率,这可以解释为丢失了一半的信息。因此,在滤波操作之后,分辨率降低了一半。但是请注意,滤波后的二次采样操作不会影响分辨率,因为从信号中删除一半的频谱分量仍然会使一半的样本数量成为多余。可以丢弃一半样本,而不会丢失任何信息。总之,低通滤波将分辨率减半,但比例不变。由于样本数量的一半是冗余的,因此该信号将被2二次采样。这使规模翻倍。

    该过程可以用数学表示为:

    y [ n ] = ∑ k = − ∞ ∞ h [ k ] ⋅ x [ 2 n − k ] y[n] = \sum\limits_{k = -\infty}^\infty h[k] \cdot x[2n - k] y[n]=k=h[k]x[2nk]·························································(4.2)

    话虽如此,我们现在看一下DWT的实际计算方法:DWT通过将信号分解为粗略的近似值和详细信息来分析具有不同分辨率的不同频带上的信号。 DWT使用两组函数,分别称为缩放函数和小波函数,分别与低通和高通滤波器关联。 通过对时域信号进行连续的高通和低通滤波,可以简单地将信号分解为不同的频带。 原始信号 x [ n ] x [n] x[n]首先通过半带高通滤波器 g [ n g [n g[n]和低通滤波器 h [ n ] h [n] h[n]。 滤波之后,根据奈奎斯特规则可以消除一半样本,因为信号现在的最高频率为 p 2 \frac {p} {2} 2p弧度,而不是p。 因此,只需丢弃其他所有采样,即可对信号进行2采样。 这构成了一个分解级别,可以用数学方式表示为:

    y h i g h [ k ] = ∑ n x [ n ] ⋅ g [ 2 k − n ] y_{high}[k] = \sum\limits_{n} x[n] \cdot g[2k - n] yhigh[k]=nx[n]g[2kn]

    y l o w [ k ] = ∑ n x [ n ] ⋅ h [ 2 k − n ] y_{low}[k] = \sum\limits_{n} x[n] \cdot h[2k - n] ylow[k]=nx[n]h[2kn]························································(4.3)

    其中 y h i g h [ k ] y_ {high} [k] yhigh[k] y l o w [ k ] y_ {low} [k] ylow[k]分别是二次采样后的高通和低通滤波器的输出。

    这种分解将时间分辨率减半,因为现在只有一半数量的样本可以表征整个信号。 但是,此操作使频率分辨率提高了一倍,因为信号的频带现在仅跨越先前频带的一半,从而有效地将频率的不确定性降低了一半。 可以重复上述过程,也称为子带编码,以进行进一步分解。 在每个级别上,滤波和二次采样将导致样本数量减少一半(因此时间分辨率减少一半),跨度频段减少一半(频率分辨率增加一倍)。 图4.1说明了此过程,其中 x [ n ] x [n] x[n]是要分解的原始信号, h [ n ] h [n] h[n] g [ n ] g [n] g[n]分别是低通和高通滤波器。 每个级别的信号带宽在图中标记为“ f f f”。

    在这里插入图片描述例如,假设原始信号 x [ n ] x [n] x[n]具有512个采样点,跨越零到p rad / s的频带。在第一个分解级别上,信号通过高通和低通滤波器,然后通过2进行二次采样。高通滤波器的输出具有256个点(因此是时间分辨率的一半),但仅跨越 p 2 \frac {p } {2} 2p至rad / s(因此频率分辨率翻倍)。这256个样本构成DWT系数的第一级。低通滤波器的输出也有256个采样,但是它跨越了频带的另一半,频率从0到 p 2 r a d / s \frac {p} {2} rad / s 2prad/s。然后,该信号通过相同的低通和高通滤波器进行进一步分解。第二个低通滤波器的输出随后进行二次采样,采样范围为0到 p 4 r a d / s \frac {p} {4} rad / s 4prad/s,范围为128个采样,第二个高通滤波器的输出之后进行二次采样,采样范围为128个采样范围 p 4 \frac {p} {4} 4p p 2 r a d / s \frac {p} {2} rad / s 2prad/s的范围。第二高通滤波信号构成第二级DWT系数。该信号的时间分辨率为第一级信号的一半,但频率分辨率为第一级信号的两倍。换句话说,与原始信号相比,时间分辨率降低了4倍,而频率分辨率提高了4倍。然后再次对低通滤波器的输出进行滤波以进行进一步分解。这个过程一直持续到剩下两个样本为止。对于此特定示例,将有8个分解级别,每个分解级别具有前一个级别的样本数量的一半。然后,通过合并从最后一个分解级别开始的所有系数(在这种情况下为两个样本),获得原始信号的DWT。然后,DWT将具有与原始信号相同数量的系数。

    在原始信号中最突出的频率将在包含那些特定频率的DWT信号区域中显示为高振幅。 该变换与傅立叶变换的区别在于,这些频率的时间定位不会丢失。 但是,时间本地化将具有取决于它们出现在哪个级别的分辨率。 如果信号的主要信息位于高频(这是最常见的情况),则这些频率的时间定位将更加精确,因为它们具有更多的样本数量。 如果主要信息仅位于非常低的频率,则时间定位将不是很精确,因为很少有样本用于表示这些频率的信号。 该过程实际上在高频下提供了良好的时间分辨率,而在低频下提供了良好的频率分辨率。 遇到的大多数实际信号都是这种类型的。

    在原始信号中不是很突出的频带将具有非常低的幅度,并且可以丢弃DWT信号的那部分而不会造成任何重大信息损失,从而可以减少数据。图4.2说明了DWT信号的外观以及如何提供数据缩减的示例。图4.2_a显示了典型的512采样信号,已将其标准化为单位幅度。水平轴是样本数,而垂直轴是归一化幅度。图4.2_b显示了图4.2_a中信号的8级DWT。此信号中的最后256个采样对应于信号中的最高频段,前128个采样对应于第二个最高频段,依此类推。应当注意,仅对应于较低分析频率的前64个样本携带相关信息,而该信号的其余部分实际上没有信息。因此,除了前64个样本外,所有其他样本都可以丢弃,而不会丢失任何信息。这就是DWT如何提供一种非常有效的数据缩减方案。

    在这里插入图片描述
    我们将重温该示例,因为它为如何解释DWT提供了重要的见解。 但是,在此之前,我们需要结束对DWT的数学分析。

    离散小波变换的一个重要特性是高通和低通滤波器的脉冲响应之间的关系。 高通和低通滤波器不是彼此独立的,它们之间的关系是:

    g [ L − 1 − n ] = ( − 1 ) n ⋅ h [ n ] g[L - 1 - n] = (-1)^n \cdot h[n] g[L1n]=(1)nh[n]··································································(4.4)

    其中g [n]是高通滤波器,h [n]是低通滤波器,L是滤波器长度(以点数为单位)。 请注意,这两个过滤器是彼此的奇数索引交替反向版本。 从低通到高通的转换由(-1)n项提供。 满足此条件的滤波器通常用于信号处理,它们被称为正交镜像滤波器(QMF)。 这两个滤波和二次采样操作可以表示为:

    y h i g h [ k ] = ∑ n x [ n ] ⋅ g [ − n + 2 k ] y_{high}[k] = \sum\limits_{n} x[n] \cdot g[-n + 2k] yhigh[k]=nx[n]g[n+2k]

    y l o w [ k ] = ∑ n x [ n ] ⋅ h [ − n + 2 k ] y_{low}[k] = \sum\limits_{n} x[n] \cdot h[-n + 2k] ylow[k]=nx[n]h[n+2k]···························································(4.5)

    由于半带滤波器形成正交基,因此这种情况下的重建非常容易。 以相反的顺序执行上述过程以进行重建。 每个级别的信号均被上采样两次,并通过合成滤波器g [n]和h [n](分别为高通和低通),然后相加。 这里有趣的一点是,除了时间反转外,分析和合成过滤器彼此相同。 因此,重建公式变为(对于每一层)

    x [ n ] = ∑ k = − ∞ ∞ (   y h i g h [ k ] ⋅ g [ − n + 2 k ]   ) + (   y l o w [ k ] ⋅ h [ − n + 2 k ]   ) x[n] = \sum\limits_{k = -\infty}^\infty \left( \, y_{high}[k] \cdot g[-n + 2k] \, \right) + \left( \, y_{low}[k] \cdot h[-n + 2k] \, \right) x[n]=k=(yhigh[k]g[n+2k])+(ylow[k]h[n+2k])··········(4.6)

    但是,如果滤波器不是理想的半带,则无法实现完美的重构。 尽管不可能实现理想的滤波器,但在某些条件下仍可以找到提供完美重构的滤波器。 最著名的是Ingrid Daubechies开发的那些,被称为Daubechies的小波。

    注意,由于连续的2次二次采样,信号长度必须是2的幂,或者至少是2的幂的倍数,以便该方案有效。 信号的长度决定了信号可以分解为的电平数。 例如,如果信号长度为1024,则可能有十个分解级别。

    解释DWT系数有时可能会很困难,因为DWT系数的表示方式非常特殊。 为了使真实情况更短,将每个级别的DWT系数从最后一个级别开始进行级联。 为了阐明这个概念,可以举一个例子:

    假设我们在10 MHZ处采样了256个样本的长信号,我们希望获得其DWT系数。由于信号以10 MHz采样,因此信号中存在的最高频率分量为5 MHz。在第一级,信号通过低通滤波器h [n]和高通滤波器g [n],其输出被二次采样。高通滤波器的输出为第一级DWT系数。其中有128个,它们代表[2.5 5] MHz范围内的信号。这128个样本是最后绘制的128个样本。低通滤波器输出(也具有128个采样点,但跨越[0 2.5] MHz频带)将通过使它们通过相同的h [n]和g [n]进一步分解。第二个高通滤波器的输出是2级DWT系数,并且这64个样本位于图中的128个1级系数之前。通过再次使第二低通滤波器的输出经过滤波器h [n]和g [n],进一步分解它。第三个高通滤波器的输出为3级DWT系数。这32个样本位于图中的2级DWT系数之前。

    该过程继续进行,直到在级别9只能计算1个DWT系数为止。这个系数是DWT图中第一个绘制的系数。随后是2个8级系数,4个7级系数,8个6级系数,16个5级系数,32个4级系数,64个3级系数,128个2级系数,最后是256个1级系数。注意,在较低频率下使用的样本数量越来越少,因此,时间分辨率随着频率的降低而降低,但是由于频率间隔在低频时也会降低,因此频率分辨率会提高。显然,仅仅由于时间分辨率大大降低,前几个系数将不会携带大量信息。为了说明这种丰富而奇怪的DWT表示,让我们看一下现实世界中的信号。我们的原始信号是一个256样本长的超声信号,它以25 MHz的频率采样。该信号最初是使用2.25 MHz传感器产生的,因此信号的主要频谱成分为2.25 MHz。最后的128个样本对应[6.25 12.5] MHz范围。从图中可以看出,这里没有可用的信息,因此可以丢弃这些样本而不会丢失任何信息。前面的64个样本表示[3.12 6.25] MHz范围内的信号,该信号也不携带任何重要信息。小故障可能对应于信号中的高频噪声。前32个样本代表[1.5 3.1] MHz范围内的信号。如您所见,正如我们预期的那样,信号的大部分能量集中在这32个样本上。前16个采样对应[0.75 1.5] MHz,在此级别看到的峰值可能表示信号的较低频率包络。先前的样本可能未包含任何其他重要信息。可以肯定地说,我们可以使用第三级和第四级系数,也就是说,我们可以用16 + 32 = 48个样本来表示这个256个样本的长信号,这大大减少了数据量,这会使您的计算机非常满意。

    从小波变换的这一特殊属性中受益最大的领域是图像处理。如您所知,图像,尤其是高分辨率图像占用大量磁盘空间。实际上,如果本教程要花很长时间下载,那主要是由于图像。 DWT可用于减小图像尺寸,而不会损失很多分辨率。方法如下:

    对于给定的图像,您可以计算每行的DWT,并丢弃DWT中所有小于某个阈值的值。然后,我们仅保存那些高于每行阈值的DWT系数,并且当我们需要重建原始图像时,我们只需为每行填充与所丢弃系数数量一样多的零,并使用反DWT来重建每行原始图像的行。我们还可以分析不同频带的图像,并仅使用特定频带的系数来重建原始图像。我将尝试尽快放入示例图像,以说明这一点。

    越来越受到关注的另一个问题是不仅在低通侧而且在两侧都进行分解(子带编码)。换句话说,分别放大信号的低频带和高频带。可以将其可视化为具有图4.1的树结构的两侧。结果就是所谓的小波包。我们不会在此讨论小波包,因为它不在本教程的讨论范围之内。任何对小波包感兴趣的人,或者对DWT有更多了解的人,都可以在市场上任何可用的文本中找到该信息。

    到此,我们的迷你系列小波教程结束了。如果我可以帮助那些努力理解小波的人,那么我将考虑花费在本教程上的时间和精力。我想提醒一下,本教程既不是小波变换的完整也不是完整的介绍。它只是小波概念的概述,旨在为那些发现小波可用文本相当复杂的人提供第一个参考。可能存在许多结构和/或技术错误,如果您能向我指出这些错误,我将不胜感激。您的反馈对于本教程的成功至关重要。

    非常感谢您对Wavelet教程的关注。

    展开全文
  • 所谓的小波的小是针对傅里叶波而言,傅里叶波指的是在时域空间无穷震荡的...傅里叶变换和小波变换的区别: ■傅里叶变换:基础函数是正弦(或余弦)函数。反映的是图像的整体特征, 其频域分析具有很好的局部性,但空间


    所谓的小波的小是针对傅里叶波而言,傅里叶波指的是在时域空间无穷震荡的正弦(或余弦波)。
    “小波”(wavelet)是一种“尺度”很小的波动,并具有时间和频率特性。
    小波函数必须满足以下两个条件:
    (1)小波必须是振荡的;
    (2)小波是能量在时域非常集中的波,它的能量有限,都集中在某一点附近,即振幅只能在一个很短的一段区间上非0,即是局部化的,且积分的值为零。如下图所示
    在这里插入图片描述
    可参考: 通俗易懂讲解小波变换

    傅里叶变换和小波变换的区别:
    ■傅里叶变换:基础函数是正弦(或余弦)函数。反映的是图像的整体特征, 其频域分析具有很好的局部性,但空间(时间)域上没有局部化功能。傅立叶变换是将信号完全的放在频率域中分析,但无法给出信号在每一个时间点的变化情况,并且对时间轴上任何点的突变都会影响整个频率的信号。
    ■小波变换:基函数是小波,具有变化的频率和有限的持续时间。是空间(时间)和频率的局部变换,它通过伸缩平移运算对信号逐步进行多尺度细化,最终达到高频处时间细分,低频处频率细分,能自动适应时频信号分析的要求,从而可聚焦到信号的任意细节。是多分辨率理论的分析基础。
    多分辨率理论:将多种学科的技术有效地统一在一起,其优势很明显,某种分辨率下所无法发现的特性在另一种分辨率下将很容易被发现。
    在这里插入图片描述
    小波的性质

    • 可分离性、尺度可变性、平移性
    • 多分辨率的一致性
    • 正交性

    一、基础

    图像通常是由相似纹理和灰度级连成的区域,它们相结合就形成了物体。

    • 若物体的尺寸较小或对比度较低,通常以较高的分辨率进行研究
    • 若物体的尺寸较大或对比度较高,粗略的观察就已足够
    • 若较大物体和较小物体同时存在时,以不同的分辨率来研究它们将更具优势,即多分辨率分析。

    以下将介绍与多分辨率分析相关的三种图像处理技术。

    1.1 图像金字塔

    是以多分辨率来表示图像的一种结构,结构非常有效,且概念简单。是一系列以金字塔形状排列的、分辨率逐步降低的图像集合。
    金字塔底部是待处理图像的高分辨率表示,顶部则包含低分辨率的近似。向金字塔上层移动时,尺寸和分辨率逐步降低。
    在这里插入图片描述
    基础级:大小为2J x 2J 或 N x N
    顶点级:大小为1x1(单个像素)
    第j级大小为2j x 2j ,0 <= j <= J。
    但通常金字塔会截短到P+1级,即将级别限制到P来降低原图像的分辨率近似,也就是说不会到金字塔靠近顶端的位置。

    第j-1级近似输出用来建立近似值金字塔;作为金字塔基级的原始图像和它的P级减少的分辨率近似都能直接获取并调整;
    第j级的预测残差输出用于建立预测残差金字塔;近似值和预测残差金字塔都通过迭代计算获得。迭代算法如下:
    1、初始化,原始图象大小,j=J
    2、j-1级,以2为步长进行子抽样,计算输入图像减少的分辨率近似值—j-1级近似值,生成子抽样金字塔。
    3、对j-1 级近似值进行步长为2的内插,并进行过滤,生成与输入图像等分辨率的预测图像。
    4、输入图像和预测图像之间的差异,产生预测残差金字塔。
    5、重复2、3、4步骤。
    在这里插入图片描述
    金字塔低分辨率级别用于分析较大的结构或图像的整体内容;高分辨率图像适合于分析单个物体的特性。

    1.2 子带编码

    一幅图像被分解为一组频带受限的分量,称为子带。
    ■ 子带可以重组在一起无失真地重建原始图象
    ■ 每个子带通过对输入进行带通滤波而得到
    ■ 子带带宽小于原始图像带宽,子带可以进行无信息损失的抽样
    ■ 原始图象的重建可以通过内插、滤波、和叠加单个子带来完成
    在这里插入图片描述
    该系统由两组滤波器组构成,h0(n)和h1(n)是半波段滤波器,它们的理想传递特性为H0和H1,用于将输入序列分成两个半长序列flp(n)和fhp(n),表示输入的子带。

    分析滤波器组即为分解,综合滤波器组即为重构

    对图像处理来说:分解就是通过下采样提取图像的低频(近似)和高频(细节)信息。每一层(尺度)的分解都是对上一层分解中的低频信息进行再分解。重构是通过上采样将分解的低频和高频信息再合并为分解前的图像的近似。

    h0(n)为低通滤波器,其输出flp(n)称为f(n)的近似;h1(n)为高通滤波器,其输出flp(n)称为高频部分或细节部分
    综合滤波器g0(n)和g1(n)将flp(n)和fhp(n)合并。
    子带编码的目的:选择h0(n)、h1(n)、g0(n)、g1(n),以便使子带编码和解码系统的输入和输出是相同的,完成这一任务时,可以说最终系统采用了完美重建滤波器

    一维滤波器用于图像处理的二维可分离滤波器。首先用于一个维度( 如垂直方向),再应用于另一个维度(如水平方向),两个阶段都执行下采样(其中一次是在第二个滤波操作之前执行的)。
    在这里插入图片描述
    2↓表示以2为基进行下采样。输出的a(m,n),dV(m,n),dH(m,n),dD(m,n)分别称为输入图像的近似子带、垂直细节子带、水平字节子带和对角线细节子带。这些子带可分为4个更小的子带,更小的子带还可再分。[a为approximation的首字母,d为detail的首字母]

    使用上图所示的编码系统对花瓶进行4子带分离。
    在这里插入图片描述

    1.3 哈尔变换(Haar)

    哈尔变换的基函数是已知的最古老、也是最简单的正交小波。
    哈尔变换可用如下矩阵形式表示:
    在这里插入图片描述
    F是一个N x N图像矩阵,H是一个N x N哈尔变换矩阵,T是一个N x N变换结果。
    H包含哈尔基函数hk(z),z∈[0,1],k = 0,1,2,…,N-1,其中N = 2n。要生成矩阵H,要定义整数k,k = 2p + q - 1,其中0<= p <= n-1,当p = 0时,q = 0或1;当p≠0时,1<= q <= 2n
    在这里插入图片描述
    N x N哈尔变换矩阵的第i行包含了元素hi(z),其中z = 0/N,1/N,2/N,…,(N-1)/N。
    如2 x 2的哈尔变换矩阵为:
    在这里插入图片描述
    当N = 4时,且假设k,p,q的值如下:
    在这里插入图片描述
    时,4 x 4的变换矩阵H4为:
    在这里插入图片描述
    H2的行可用于定义一个2抽头完美重建滤波器组(见子带编码)的分析滤波器h0(n)和h1(n),以及最简单且最古老的小波变换的缩放比例和小波向量。

    二、多分辨率展开

    在多分辨率分析(MRA)中,尺度函数被用于建立一个函数或一幅图像的一系列近似,相邻两近似之间的分辨率相差2倍。使用称为小波的附加函数来对相邻近似之间的差进行编码。

    2.1 级数展开

    一个信号或函数f(x)可展开为函数的线性组合:
    在这里插入图片描述
    k为有限和或无限和的整数下标,αk为展开系数,φk(x)为展开函数。
    展开唯一,即任意给定的f(x)只有一组αk与之对应。则φk(x)称为基函数,集合{φk(x)}称为可表示这样一类函数的
    可展开的函数形成了一个函数空间,称为展开集合的闭合跨度,表示如下:
    在这里插入图片描述
    对任意函数空间V及相应的展开集合{φk(x)},都有一个表示为{φk ̃(x)}的对偶函数集合。展开系数αk即是对偶函数φk ̃(x)和函数f(x)的内积。
    在这里插入图片描述

    2.2 尺度函数

    考虑由实、平方可积函数φ(x)的整数平移和二值尺度组成的展开函数集合{φj,k(x)}
    在这里插入图片描述
    其中φj,k(x)对所有的j,k∈Z(整数集)和φ(x)属于L2(R)都成立(L2(R)表示度量的、平方可积的一维函数集合,R为实数集)。
    整数平移k决定了φj,k(x)沿x轴的位置;尺度j决定了φj,k(x)的宽度;2j/2控制函数的幅度。
    因为φj,k(x)的形状随j发生变化,故φ(x)称为尺度函数
    选择适当的φ(x),可使{φj,k(x)}张成L2(R)。
    对任意的j,将k上张成的子空间表示为:
    在这里插入图片描述
    增加j就会增加Vj的大小,进而允许子空间中包含具有更小变量或更细细节的函数。因为随着j的增大,用于表示子空间函数的φj,k(x)会变窄,且x有较小变化就可以分开。
    在这里插入图片描述
    hφ(n)为小波函数系数,hφ为小波向量。

    2.3 小波函数

    满足MRA要求的尺度函数被定义为小波函数ψ(x),它与其整数平移及二值尺度一起,跨越了任意连个相邻尺度子空间Vj和Vj+1之间的差。
    在这里插入图片描述
    对跨越图中Wj空间的所有k∈Z,定义小波集合{ψj,k(x)}
    在这里插入图片描述
    使用尺度函数,可以写出:
    在这里插入图片描述
    如果f(x)∈Wj,则
    在这里插入图片描述
    尺度函数和小波函数子空间由下式联系起来
    在这里插入图片描述
    ⨁表示空间并集,Vj+1中Vj的正交补集是Wj,且Vj中的所有成员与Wj中的所有成员都正交。故
    在这里插入图片描述
    所有可度量的、平方可积的函数空间可以表示为
    在这里插入图片描述
    若f(x)是V1而非V0的元素,则(1)式的展开包含使用V0尺度函数的f(x)的近似;来自W0的小波将对这种近似于实际函数之间的差进行编码。故可推广:
    在这里插入图片描述
    j0是任意开始尺度。
    因为小波空间存在于由相邻较高分辨率尺度函数跨越的空间中,故小波函数可表示成平移后的双倍分辨率尺度函数的加权和,表示如下:
    在这里插入图片描述
    hψ(n)为小波函数系数,hψ为小波向量。
    利用小波跨越正交补集空间和整数小波平移是正交的条件,可得hψ(n)和hφ(n)按下述方式相关:
    在这里插入图片描述

    三、小波变换

    3.1 一维小波变换

    3.1.1 小波级数展开

    与小波ψ(x)和尺度函数φ(x)相关的函数f(x)∈L2(R)的小波级数展开如下:
    在这里插入图片描述
    j0是任意的开始尺度,cj0(k)和dj(k)分别是尺度函数和小波函数下的展开系数αk的改写形式。cj0(k)称为近似和/或尺度系数,dj(k)称为细节和/小波系数。
    哈哈啊
    在这里插入图片描述
    即为被展开的函数和展开函数的内积。

    3.1.2 离散小波变换

    上一小节的小波系数展开将一个连续变量函数映射为一系列系数。若待展开的函数是离散的(即数字序列),则得到的系数就称为离散小波变换(DWT)。则序列f(n)的正向DWT系数(类似于上一节的cj0(k)和dj(k))如下:
    在这里插入图片描述
    其中φj0,k(n)ψj,k(n)是基函数φj0,k(x)ψj,k(x)的取样形式,n = 0,1,2,…,M-1。所以f(n)的展开如下:
    在这里插入图片描述

    3.1.3 连续小波变换

    离散小波变换的自然延伸是连续小波变换(CWT),连续小波变换将一个连续函数变换为两个连续变量(平移和尺度)的高冗余度函数。得到的变换易于解释并且对于时间-频率分析是有价值的。我们感兴趣的是离散图像,此处只是为了完整性。
    连续平方可积函数f(x)的连续小波变换与实数值小波ψ(x)的关系定义如下:
    在这里插入图片描述
    其中
    在这里插入图片描述
    s和τ分别为尺度参数和平移参数
    f(x)可由连续小波反变换求得:
    在这里插入图片描述
    其中
    在这里插入图片描述
    式中Ψ(μ)是ψ(x)的傅里叶变换。

    3.2 快速小波变换

    快速小波变换(FWT)是一种实现离散小波变换(DWT)的高效计算。类似于2子带的子带编码方案。

    正变换,即分解
    先定义尺度函数和小波函数,再计算尺度j的近似和细节系数。

    3.2.1 一维情况

    由前面描述的尺度函数和小波函数
    在这里插入图片描述
    其中hφ(n)为尺度函数系数,hφ为尺度向量;hψ(n)为小波函数系数,hψ为小波向量。
    用2j对x进行尺度化,用k对它进行平移,并令m = 2k + n,则
    尺度函数:
    在这里插入图片描述
    小波函数
    在这里插入图片描述
    近似系数和细节系数如下:
    在这里插入图片描述
    hφ(-n) = h0(n)(2带子带编码和解码系统中的分析部分),hψ(-n) = h1(n)
    在这里插入图片描述
    在这里插入图片描述
    反变换,即重构
    快速小波反变换
    反变换即综合滤波器,综合滤波器和分析滤波器之间的顺序是相反的。
    g0(n) =h0(-n) = hφ(n) =,g1(n) = h1(-n) = hψ(n)
    在这里插入图片描述
    在这里插入图片描述

    3.3 二维小波变换

    在二维情况下,需要一个二维尺度函数φ(x,y)和三个二维小波ψH(x,y),ψV(x,y),ψD(x,y)。每个二维小波都是两个一维函数的乘积。排除产生一维结果的乘积,如φ(x)ψ(x),4个剩下的乘积可产生尺度函数
    在这里插入图片描述
    和可分的“方向敏感”小波
    在这里插入图片描述
    这些小波度量函数的变化——图像的灰度变化——沿不同方向的变化:ψH度量沿列方向的变化(如水平边缘),ψV响应沿行方向的变化(如垂直边缘),ψD度量对应对角线方向的变化。
    定义尺度和平移基函数:
    在这里插入图片描述
    故大小为M*N的图像f(x,y)的离散小波变换是:
    在这里插入图片描述
    通常令j0 = 0,并选择N = M = 2J,故j = 0,1,2,…,J-1和m = n = 0,1,2,…,2j-1。
    f(x,y)可通过离散小波反变换得到:
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述

    3.4 快速小波变换的实现

    3.4.1 使用小波工具箱的快速小波变换

    wfilters 小波滤波器

    [Lo_D,Hi_D,Lo_R,Hi_R] = wfilters(wname)
    

    输出Lo_D,Hi_D,Lo_R和Hi_R是行向量,分别返回与wname相关低通分解、高通分解、低通重构和高通重构滤波器。wname是正交或双正交小波的名字,决定返回的滤波器系数,下表中一一列出
    在这里插入图片描述

    [F1,F2] = wfilters(wname,type)
    

    返回一对与wname相关的type类型滤波器。type可取’d’、‘r’、‘l’、‘h’四种值。
    分解或低通滤波器在F1中返回,重构或高通滤波器放在F2中。
    如:wfilters(‘db6’,‘h’) 返回一对与db6相关的高通滤波器(包括高通分解和高通重建)
    在这里插入图片描述
    waveinfo 查看小波家族(即wname中的内容)

    waveinfo(wfamily)
    

    wavefun 小波和尺度函数

    [phi,psi,xval] = wavefun(wname,iter)
    % 对正交变换。返回与wname相关的小波和尺度函数的近似向量phi和psi,近似向量是在网格点xval处估计的。正整数iter控制计算的迭代次数,可决定近似值的精度。
    [phi1,psi1,phi2,psi2,xval] = wavefun(wname,iter)
    % 对双正交变换。phi1和psi1对应分解函数,phi2和psi2对应重构函数。
    [psi,xval] = wavefun(wname,iter)
    % 对没有尺度函数的小波,如Morlet小波、Mexcian小波、Gaussian derivatives小波和复小波。
    
    % 二维小波和尺度函数(只用于正交小波)
    [S,W1,W2,W3,XYVAL] = wavefun2('wname',ITER)   % S为尺度函数,W1,W2,W3为三个“方向敏感”的小波函数
    [S,W1,W2,W3,XYVAL] = wavefun2('wname',ITER,'plot')  % 计算并画出函数图像
    

    dwt2 二维离散小波变换(单尺度)

    [cA,cH,cV,cD]=dwt2(X,'wname')   %使用指定的小波基函数对矩阵X进行二维离散小波变换
    [cA,cH,cV,cD]=dwt2(X,Lo_D,Hi_D)  %使用指定的低通滤波器Lo_D和高通滤波器Hi_D分解信号
                                     % cA--近似分量(低频分量) cH--水平方向细节分量 cV--垂直方向细节分量;cD--对角方向细节分量
    

    示例

    A=imread('C:\Users\win\Desktop\girl.jpg');
    [cA,cH,cV,cD]=dwt2(A,'haar');%使用haar小波
    figure,imshow(A);title('原图');
    figure,subplot(2,2,1);imshow(cA/255);title('低频分量');
    subplot(2,2,2),imshow(cH),title('水平细节分量');
    subplot(2,2,3),imshow(cV),title('垂直细节分量');
    subplot(2,2,4),imshow(cD),title('对角线细节分量');
    

    在这里插入图片描述
    在这里插入图片描述
    wavedec2 二维多尺度分解

    wavedec2  可直接处理彩色图像
    [C,S] = wavedec2(X,N,wname)  % 对图像X用wname小波基函数实现N尺度(级数)分解。
    [C,S] = wavedec2(X,N,Lod,Hid) % 用特定的低通和高通滤波器(Lod和Hid)进行小波分解。
    

    对第一种情况进行解析
    经过小波分解之后得到的所有图像都被称为小波系数。C为小波分解向量,即各层分解系数,S为各层分解系数长度,也就是大小.。
    C的结构:
    C=[A(N)|H(N)|V(N)|D(N)|H(N-1)|V(N-1)|D(N-1)|H(N-2)|V(N-2)|D(N-2)|…|H(1)|V(1)|D(1)]
    C是一个行向量,即:1*(size(X)),(e.g,X=256256,then c大小为:1(256256)=165536)
    A(N)代表第N层低频系数,即图像第N层的近似,尺度最小,在金字塔中就是每层的下采样的图像,在金字塔中,在第N-1层下采样到N层,N层的图像维度(尺度)变小了,也就意味着在下采样过程中丢失了信息,丢失的信息实质是高频信息,这些信息在小波分解中可以通过HVD这些高频分量来保存。
    H(N)|V(N)|D(N)代表第N层高频系数,分别是水平,垂直,对角高频,以此类推,到H(1)|V(1)|D(1).故C包含3N+1个部分。
    每一次的小波分解都是在近似图像上进行分解。a为近似图像
    在这里插入图片描述
    (l)LL子带是图像的近似表示。
    (2)HL子带是水平子带
    (3)LH子带是垂直子带
    (4)HH子带是对角子带
    S的结构:
    S是一个大小为(N+2)*2的记录数组,储存各层分解系数大小。其形式为:
    S = [sa(N,:) sd(N,:) sd(N-1,:) … sd(1,:) sX]
    sa(N,:) sd(i,:) sX分别包含第N级的近似矩阵A(N)的水平和垂直维数,第i尺度的细节(水平、垂直和对角线)和原图像X的大小。
    即第一行是A(N)的大小,第二行是H(N)|V(N)|D(N)|的大小,第三行是 H(N-1)|V(N-1)|D(N-1)的大小,倒数第二行是H(1)|V(1)|D(1)大小,最后一行是X的大小。S中的元素是以列向量的形式组织的。
    示例

    X=imread('C:\Users\win\Desktop\girl.jpg');
    [c,s]=wavedec2(X,2,'db1');%进行2尺度二维离散小波分解。分解小波函数:db1
    [cH1,cV1,cD1]=detcoef2('all',c,s,1);%尺度1的所有方向的高频系数
    [cH2,cV2,cD2]=detcoef2('all',c,s,2);%尺度2的所有方向的高频系数
    cA1=appcoef2(c,s,'db1',1);%尺度1的低频系数
    cA2=appcoef2(c,s,'db1',2);%尺度2的低频系数
    
    figure;
    subplot(2,2,[1,2]),imshow(X);title('原图');
    subplot(2,2,3),imshow(uint8(cA1));title('尺度1的低频系数图像');
    subplot(2,2,4),imshow(uint8(cA2));title('尺度2的低频系数图像');
    figure;
    subplot(2,3,1),imshow(uint8(cH1));title('尺度1水平方向高频系数');
    subplot(2,3,2),imshow(uint8(cV1));title('尺度1垂直方向高频系数');
    subplot(2,3,3),imshow(uint8(cD1));title('尺度1斜线方向高频系数');
    subplot(2,3,4),imshow(uint8(cH2));title('尺度2水平方向高频系数');
    subplot(2,3,5),imshow(uint8(cV2));title('尺度2垂直方向高频系数');
    subplot(2,3,6),imshow(uint8(cD2));title('尺度2斜线方向高频系数');
    

    在这里插入图片描述
    在这里插入图片描述
    快速小波 变换可能会导致边界失真,为了将失真减到最小,需要对边界进行填充扩展。可使用dwtmode(离散小波变换扩展模式)来扩展及填充待处理的图像。

    dwtmode(mode)  % 设置离散小波变换和小波包变换的图像扩展模式
    st = dwtmode('status') % 显示当前模式并返回到st
    st = dwtmode('status','nodisp')  % 返回当前模式st,并且在MATLAB命令窗口中不显示任何状态或警告文本。
    dwtmode('save',mode)  % 保存mode作为新的默认模式
    

    在这里插入图片描述
    例: [cA,cH,cV,cD] = dwt2(x,‘db4’,‘mode’,‘per’);

    3.4.2 不使用小波工具箱的快速小波变换

    将使用wavefilter和wavefast函数来代替小波工具箱函数wfilters和wavedec2函数。
    wavefilter.m

    function [varargout] = wavefilter(wname,type)
    % wavefilter 生成小波分解和重构滤波器
    % 返回分解和/或重构滤波器
    % [ld,hd,lr,hr] = wavefilter('haar') 得到基于haar小波的低和高分解滤波器(ld,hd)和重构滤波器(lr,hr)
    % [ld,hd] = wavefilter('haar','d') 得到分解滤波器
    % [lr,hr] = wavefilter('haar','r') 得到重构滤波器
    % type:d(分解)和r(重构)
    
    % 检查输入和输出参数
    narginchk(1,2); % narginchk 来检查输入参数的个数在期望的范围内.输入参数nargin个数小于1或大于2则产生错误
    if (nargin == 1 && nargout ~= 4) || (nargin ==2 && nargout ~= 2)
        error('输出参数数量无效')
    end
    if nargin == 1 && ~ischar(wname)
        error('wname必须是字符串');
    end
    if nargin == 2 && ~ischar(type)
        error('type必须是字符串');
    end
    
    % 生成滤波器
    switch lower(wname)
        case {'haar','db1'}   % 标准正交滤波器
            ld = [1 1] / sqrt(2);
            hd = [-1 1] /sqrt(2);
            lr = ld;
            hr = -hd;
        case'db4' % 标准正交滤波器
            ld=[-1.059740178499728e-002 3.288301166698295e-00 ...
                3.084138183598697e-002 -1.870348117188811e-001...
                -2.79837641698385e-002 6.308807679295904-001 ...
                7.148465705525415e-001 2.303778133088552e-001];
            t=(0:7);
            hd=ld; hd(end:-1:1)=cos(pi*t).*ld;
            lr=ld; lr(end:-1:1)=ld;
            hr=cos(pi*t).*ld;
            
        case'sym4'  % 标准正交滤波器
            ld=[-7.576571478927333e-002 -2.963552764599851e-002...
                4.976186676320155e-001 8.037387518059161e-001...
                2.978577956052774e-001 -9.921954357684722e-002...
                -1.260396726203783e-002 3.222310060404270e-002];
            t=(0:7);
            hd=ld; hd(end:-1:1)=cos(pi*t).*ld;
            lr=ld; lr(end:-1:1)=ld;
            hr=cos(pi*t).*ld;
            
        case'bior6.8'  % 双正交滤波器
            ld=[0 1.908831736481291e-003 -1.914286129088767e-003...
                -1.699063986760234e-002 1.193456527972926e-002...
                4.973290349094079e-002 -7.726317316720414e-002...
                -9.405920439573646e-002 4.207962846098268e-001...
                8.259229974584023e-001 4.207962846098268e-001...
                -9.405920920349573646e-002 -7.726317316720414e-002...
                4.973290349094079e-002 1.19345652792926e-002...
                -1.699063986760234e-002 -1.914286129088767e-003...
                1.908831736481291e-003];
            hd=[0 0 0 1.442628250562444e-002 -1.446750489679015e-002...
                -7.872200106262882e-002 4.036797903033992e-002...
                4.178491091502746e-001 -7.589077294536542e-001...
                4.178491091502746e-001 4.036797903033992e-002...
                -7.872200106262882e-002 -1.446750489679015e-002...
                1.442628250562444e-002 0 0 0 0];
            t=(0:17);
            lr=cos(pi*(t+1)).*hd;
            hr=cos(pi*t).*ld;
            
        case'jpeg.7'  % 双正交滤波器
            ld=[0 0.2674875741080976 -0.01686411844287495...
                -0.07822326652898785 0.2668641184428723...
                0.6029490182363579 0.2668641184428723...
                -0.07822326652898785 -0.01686411844287495...
                0.02674875741070976];
            hd=[0 -0.9127176311424948 0.05754352622849957...
                0.5912717631142470 -1.115087052456994...
                0.5912717631142470 0.05754352622849957...
                -0.09127176311424948 0 0];
            t=(0:9);
            lr=cos(pi*(t+1).*hd);
            hr=cos(pi*t).*ld;
            
        otherwise
            error('未识别小波名称');
    end
    
    %输出需要的滤波器
    if (nargin==1)
        varargout(1:4)={ld,hd,lr,hr};
    else 
        switch lower(type(1))
            case'd'
                varargout={ld,hd};
            case'r'
                varargout={lr,hr};
            otherwise
                error('未识别滤波器类型.');
        end
    end  
    

    wavefast.m
    注: wavefast用于处理灰度图像

    function [c,s] = wavefast(x,n,varargin)
    % wavefast 多尺度2维快速小波变换(即小波分解)
    % [C,L] = wavefast(x,n,lp,hp) 基于分解滤波器lp和hp,图像x的二维n尺度下的快速小波变换
    % [C,L] = wavefast(x,n,wname) 小波wname使用wavefilter得到的分解滤波器(lp和hp)
    % n必须小于或等于图像的尺寸最大值的log2.滤波器的lp和hp必须是偶数
    % 为减少边界失真,x是对称扩展的。若x = [c1 c2 c3 … cn](1维情况下),则对称扩展为[… c3 c2 c1 c1 c2 c3 … cn cn cn-1 cn-2 …]
    % 输出:
    % 矩阵C是一个系数分解向量   C = [a(n) h(n) v(n) d(n) h(n-1) … v(1) d(1)]
    % a,h,v,d分别是近似、水平、竖直和对角线系数。C有3n+1部分,其中n是小波分解数
    % 矩阵S是一个(n + 2) * 2 的记录矩阵
    % S = [sa(n,:);sd(n,:);sd(n-1,:);… sd(1,:);sx] 其中sa和sd是近似和细节尺寸
    
    % 检查输入参数的合理性
    narginchk(3,4);
    if nargin == 3
        if ischar(varargin{1}) % 判断wname是否是字符串
            [lp,hp] = wavefilter(varargin{1},'d');
        else
            error('错误的小波名称');
        end
    else
        lp = varargin{1};hp = varargin{2};
    end
    fl = length(lp);sx = size(x);
    if (ndims(x) ~= 2)||(min(sx) < 2)|| ~isreal(x)|| ~isnumeric(x)  % ndims:数组维度数目
        error('x必须是实数矩阵');
    end
    if (ndims(lp) ~= 2)|| ~isreal(lp)|| ~isnumeric(lp)...
            || (ndims(hp) ~= 2)|| ~isreal(hp)|| ~isnumeric(hp)...
            || (f1 ~= length(hp))|| rem(fl,2) ~= 0     % rem(x,y):求整除x/y的余数
        error('lp和hp必须是偶数并且长度相等, ... 数值型滤波器向量.');
    end
    if ~isreal(n)||~isnumeric(n)||(n<1)||(n > log2(max(sx)))
        error('n必须是个实数标量介于1到 ... log2(max(size(x)))');
    end
    
    % 设置初始输出数据结构和初始近似值
    c = [];s = sx;app = double(x);
    % 对于每个分解
    for i = 1:n
        % 对称地扩展近似值
        [app,keep] = symextend(app,fl);
        % 用hp和下采样计算卷积行。然后用hp和lp得到的卷积列去获取对角线和垂直的系数。
        rows = symconv(app,hp,'row',fl,keep);
        coefs = symconv(rows,hp,'col',fl,keep);
        c = [coefs(:)' c]; s = [size(coefs);s];
        coefs = symconv(rows,lp,'col',fl,keep);
        c = [coefs(:)' c];
        % 用hp和采样率计算卷积行。然后用hp和lp得到的卷积列去获取垂直的和下一个近似的系数。
        rows=symconv(app,lp,'row',fl,keep);
        coefs=symconv(rows,hp,'col',fl,keep);
        c=[coefs(:)' c];
        app=symconv(rows,lp,'col',fl,keep);
    end
    
    %追加最后的近似值结构
    c=[app(:)' c]; s=[size(app);s];
    
    %---------------------------------------------------------------%
    function [y,keep] = symextend(x,fl)
    % 计算卷积和下采样后的系数放在keep中,然后在两个维度上扩展x
    keep = floor((fl + size(x) - 1) / 2);
    y = padarray(x,[(fl - 1) (fl - 1)],'symmetric','both'); 
        % fp = padarray(f,[r,c],method,direction)对图像进行填充 
        % f为输入图像,fp为填充后的图像,[r,c]给出填充f的行数和列数,method表示填充方法,direction表示填充的方向
        % 若参数中不包括direction,则默认为'both'(在每一维的第一个元素前和最后一个元素后填充)
    %----------------------------------------------------------------%
    function y = symconv(x,h,type,fl,keep)
    % 用h和下采样,卷积x的行或列,并提取对称扩展的中心部分 
    if strcmp(type,'row')
        y = conv2(x,h); % conv2 :二维卷积
        y = y(:,1:2:end); % 抛弃偶数索引行或列(即以2下取样)
        y = y(:,fl / 2 + 1:fl / 2 + keep(2)); % 提取每行或列的中心keep元素
    else
        y = conv2(x,h');
        y = y(1:2:end,:);
        y = y(fl / 2 + 1:fl / 2 + keep(1),:);
    end
    

    当产生实际上相同的结果时,自定义函数wavefast的速度要比对应的wavedec2函数几乎快两倍。

    3.5 小波分解结构的处理

    3.5.1 用于执行变换分解向量C的小波工具箱函数

    >> X = magic(8);
    >> [C,S] = wavedec2(X,3,'haar');
    >> size(C)
    ans =
         1    64
    >> S
    S =
         1     1
         1     1
         2     2
         4     4
         8     8
    

    使用函数wavedec2进行3尺度haar小波分解。C的大小为1 x 64,S的大小为5 x 2。因为是N = 3级(尺度)分解,故C中会有3N + 1 = 10个近似和细节子矩阵的元素。基于S,这些子矩阵包括:(a)分解级别3:一个大小为1 x 1的近似矩阵和3个大小为1 x 1的细节矩阵;(b)分解级别2:3个大小为2 x 2的细节矩阵;©分解级别1:3个大小为4 x 4的细节矩阵。最后一行包含了原图像X的尺寸大小。

    appcoef2 提取二维近似(低频)系数

    A = appcoef2(C,S,wname) % 使用二维小波分解结构和wname指定的小波,提取近似系数
    A = appcoef2(C,S,LoR,HiR) % 使用低通重构和高通重构滤波器提取近似系数
    A = appcoef2(___,N) % 返回第N尺度的近似系数
    

    detcoef2 提取二维细节(高频)系数

    D = detcoef2(O,C,S,N)  % 从小波分解结构中提取水平、竖直或对角线细节系数,O = ‘h’(或'v'或'd'),N为尺度,必须是整数,且1≤N≤size(S,1)-2。
    [H,V,D] = detcoef2('all',C,S,N)  % 返回N尺度上的水平、竖直和对角线细节系数。等同于detcoef2('a',C,S,N)
    D = detcoef2('compact',C,S,N)  % 返回N尺度上的细节系数,并以行的形式存储。
    

    wthcoef2 二维小波系数阈值化

    NC = wthcoef2(type,C,S,N,T,SORH)  % type为近似系数‘a’或细节系数‘h’、‘v’、‘d’。N是分解等级的一个向量,该向量将基于向量T中的相应阈值进行阈值处理。SORH针对软硬阈值处理分别设置为's'或'h'。若省略T,则所有满足type和N规范的系数都将赋值为零。
    

    3.5.2 不使用小波工具箱编辑(提取、置零或修改)小波分解系数

    wavework.m

    function [varargout] = wavework(opcode,type,c,s,n,x)
    % wavework 用于编辑小波分解系数
    % 通过type和尺度n获取分解系数,基于opcode进行访问或修改
    
    % opcode      
    % 'copy'     [varargout] = Y  获得需要的类型type第n尺度的系数矩阵
    % 'cut'      [varargout] = [NC,Y]  将type类型下c中的分解系数Y置零,保存在新的分解向量NC中
    % 'paste'    [varargout] = [NC]   % 用x的元素代替type类型下c中的分解系数,并将新的分解结构放到NC中
    
    % type   近似系数'a'或细节系数'h','v','d'
    % [c,s]  是小波工具箱分解结构
    % n      是分解尺度(分解级数)(若type = 'a'则忽略)
    % x      是用于复制的二维系数矩阵
    
    % 检查输入参数的合理性
    narginchk(4,6);
    
    if (~ismatrix(c)) || (size(c,1)~=1)
        error('C必须是一个行向量.')
    end
    
    if (~ismatrix(s) || ~isreal(s) || ~isnumeric(s) || (size(s,2)) ~=2)
        error('S必须是一个实数数值型的两列数组');
    end
    
    elements=prod(s,2); %c的每个系数子矩阵中的元素数。 prod:数组元素的乘积,结果为包含每一行乘积的列向量
    if (length(c)<elements(end)) || ...
            ~(elements(1)+3*sum(elements(2:end-1))>=elements(end))
        error(['[c,s] 必须是一个标准的小波分解' '结构.']);
    end
    
    if nargin<5
        n=1;%默认尺度(分级)为1
    end
    nmax=size(s,1)-2; % 在[c,s]下最大尺度(分级)
    aflag=(lower(type(1))=='a');
    if ~aflag && (n>nmax)
        error('N超过了[c,s]下的分解.'); % 尺度大过了分解的最大尺度
    end
    
    switch lower(type(1))  %指向type和n的相关系数的一对指针
        case 'a' 
            nindex=1;
            start=1; % 开始索引
            stop=elements(1); %结束索引 近似矩阵中的元素的数量
            ntst=nmax;
        case{'h','v','d'}
            switch type
                case 'h',offset=0; %细节的补偿值
                case 'v',offset=1;
                case 'd',offset=2;
            end
            nindex=size(s,1)-n; %索引到详细信息
            start=elements(1)+3*sum(elements(2:nmax-n+1))+...
                   offset*elements(nindex)+1;
            stop=start+elements(nindex)-1;
            ntst=n;
        otherwise
            error('TYPE 必须是 ''a'',''h'',''v'',或''d''.');
    end
    
    switch lower(opcode)    %执行opcode要求的操作
        case{'copy','cut'}  % copy:提取type类型下的c的分解系数并复制到y中; cut:提取type类型下c的分解系数,将其置零并复制到y中。而y是一个已预分配好的二维矩阵,其大小由s决定
            y=repmat(0,s(nindex,:));
            y(:)=c(start:stop);nc=c;
            if strcmpi(opcode(1:3),'cut')
                nc(start:stop)=0;varargout={nc,y};
            else
                varargout={y};
            end
        case'paste' 
            if numel(x) ~= elements(end-ntst)
                error('X is not sized for the requested paste.');
            else
                nc=c;nc(start:stop)=x(:);varargout={nc};
            end
        otherwise
            error('不能识别OPCODE.')
    end
    

    wavecopy、wavecut和wavepaste分别使用wavework处理(提取,置零或修改)小波分解结构中的C

    wavecopy.m

    function y = wavecopy(type,c,s,n)
    % wavecopy 获取第n尺度(级别)小波分解结构的系数矩阵
    % 返回基于type和n的系数矩阵
    % type  :近似系数'a'或细节系数'h','v,'d'
    % [c,s]  是小波分解结构
    % n是分解尺度(级数)(若type = 'a'则忽略)
    
    % 检查输入参数的合理性
    narginchk(3,4);
    if nargin == 4
        y = wavework('copy',type,c,s,n);
    else
        y = wavework('copy',type,c,s);
    end
    

    wavecut.m

    function [nc,y] = wavecut(type,c,s,n)
    % wavecut 小波分解结构中的零系数
    % 返回细节或近似系数矩阵y置零的新的分解向量nc(基于type和n)。
    % type  :近似系数'a'或细节系数'h','v,'d'
    % [c,s]  是小波分解结构
    % n是分解尺度(级数)(若type = 'a'则忽略)
    
    % 检查输入参数的合理性
    narginchk(3,4);
    if nargin == 4
        [nc,y] = wavework('cut',type,c,s,n);
    else
        [nc,y] = wavework('cut',type,c,s);
    end
    

    wavepast.m

    function nc = wavepaste(type,c,s,n,x)
    % wavepaste 将type类型的系数矩阵x复制于小波分析结构的c中
    % 基于type和n,将x复制到c里面,然后返回一个新的分解结构
    % type  :近似系数'a'或细节系数'h','v,'d'
    % [c,s]  是小波分解结构
    % n是分解尺度(级数)(若type = 'a'则忽略)
    % x 是一个二维近似或细节系数矩阵,矩阵的维度适合于分解尺度(级别)n
    
    narginchk(5,5);
    nc = wavework('paste',type,c,s,n,x);
    

    再次说明:
    wavecopy:通过type和尺度n获取小波分解结构C中对应的分解系数
    wavecut:通过type和尺度n获得小波分解结构C中对应的分解系数y之后,将其置零,然后返回新的小波分解结构NC和y
    wavepaste:通过type和尺度n获得小波分解结构C中对应的分解系数,然后使用x对其进行替换,返回新的小波分解结构NC

    示例

    >> f = magic(8);
    >> [c1,s1] = wavedec2(f,3,'haar');  % 获得小波分解结构
    >> approx = wavecopy('a',c1,s1)  % 提取近似系数
    approx =
      260.0000 
    >> horizdet2 = wavecopy('h',c1,s1,2)  % 提取水平细节系数
    horizdet2 =
       1.0e-13 *
             0   -0.2842
             0         0
    >> [newc1,horizdet2] = wavecut('h',c1,s1,2);  % 将水平细节系数置零,并返回到新的小波分解结构newc1中
    >> horizdet2
    horizdet2 =
       1.0e-13 *
             0   -0.2842
             0         0
    >>  newhorizdet2 = wavecopy('h',newc1,s1,2)  % 提取newc1中的水平细节系数
    newhorizdet2 =
         0     0
         0     0
    

    3.5.3 显示小波分解系数(显示小波变换处理后的结果)

    wave2gray 用于显示小波分解系数。小波分解系数就是不同尺度n下的分解图像。
    wave2gray.m

    function w = wave2gray(c,s,scale,border)
    %wave2gray 显示小波分解系数
    %   显示并返回一个小波系数图像
    %	例子:
    %			wave2gray(c,s);						显示w或默认值
    %			foo = wave2gray(c,s);				显示并返回
    %			foo = wave2gray(c,s,4);				放大细节
    %			foo = wave2gray(c,s,-4);			放大的绝对值
    %			foo = wave2gray(c,s,1,'append');	存储边界值
    
    %	[c,s] 小波分解结构
    %
    %	scale			细节系数的比例
    %--------------------------------------------------------------------
    %	0 or 1         最大值范围 (默认)
    %	2, 3...        通过比例因子放大默认值
    %	-1, -2...      通过abs(scale)放大绝对值
    %
    %	border			小波分解间的边界
    %--------------------------------------------------------------------
    %	'absorb'		边界替代图像(默认)
    %	'append'		边界增加图像的宽度
    %
    %	Image w:	    ------- ------ -------------- -------------------
    %					|      |      |              |
    %					| a(n) | h(n) |              |
    %					|      |      |              |
    %					------- ------     h(n-1)    |
    %					|      |      |              |
    %					| v(n) | d(n) |              |      h(n-2)
    %					|      |      |              |
    %					------- ------ --------------
    %					|             |              |
    %					|    v(n-1)   |    d(n-1)    |
    %					|             |              |
    %					-------------- -------------- -------------------
    %					|                            |
    %					|           v(n-2)           |      d(n-2)
    %					|                            |
    
    % 检查输入参数的合理性
    narginchk(2,4);
    
    if (~ismatrix(c)) || (size(c, 1) ~= 1)
    	error('C必须是一个行向量.'); 
    end
    
    if (~ismatrix(s)) || ~isreal(s) || ~isnumeric(s) || (size(s, 2) ~= 2)
    	error('S必须是一个实数数值型的两列数组.');
    end
    
    elements = prod(s, 2);
    
    if (length(c) < elements(end)) ||   ...
    		~(elements(1) + 3 * sum(elements(2:end - 1)) >= elements(end))
    	error(['[C S] 必须是一个标准小波 ' ...
    				'分解结构.']);
    end
    
    if (nargin > 2) && (~isreal(scale) || ~isnumeric(scale))
    	error('SCALE必须是一个实数、数值型标量.');
    end
    
    if (nargin > 3) && (~ischar(border))
    	error('BORDER必须是字符串.');
    end
    
    if  nargin == 2
    	scale =1;    % 默认比例.
    end
    
    if nargin < 4
    	border = 'absorb'; %默认边界.
    end
    
    %比例系数并确定填充量.
    absflag = scale < 0;
    scale = abs(scale);
    if scale == 0
    	scale = 1;
    end
    
    [cd, w] = wavecut('a', c, s);  w = mat2gray(w); % mat2gray 将矩阵转换为灰度图像
    cdx = max(abs(cd(:))) / scale;
    if absflag
    	cd = mat2gray(abs(cd), [0, cdx]); fill = 0;
    else
    	cd = mat2gray(cd, [-cdx, cdx]); fill = 0.5;
    end
    
    % 一次建立一个灰度图像的分解
    for i = size(s, 1) - 2:-1:1
    	ws = size(w);
    	h = wavecopy('h', cd, s, i);
    	pad = ws - size(h);    frontporch = round(pad / 2);
    	h = padarray(h, frontporch, fill, 'pre');
    	h = padarray(h, pad - frontporch, fill, 'post');
    	v = wavecopy('v',   cd,   s,   i);
    	pad = ws - size(v);            frontporch = round(pad  /  2);
    	v = padarray(v,   frontporch,   fill,   'pre');
    	v = padarray(v,  pad - frontporch,  fill,   'post');
    	d = wavecopy('d',   cd,   s,   i);
    	pad = ws - size(d);            frontporch = round(pad  /  2);
    	d = padarray(d,   frontporch,   fill,   'pre');
    	d = padarray(d,   pad - frontporch,   fill,   'post');
       % 加一像素的白色边框
    	switch  lower(border)
    		case   'append'
    			w =  padarray(w,   [1   1],   1,   'post');
    			h = padarray(h,   [1   0],   1,   'post');
    			v = padarray(v,   [0  1],   1,   'post');
    		case   'absorb'
    			w(:,   end)   =  1;         w(end,   :)   =  1;
    			h(end,   :)   =  1;          v(:,   end)   =  1;
    		otherwise
    			error('无法识别边界参数.');
    		end
    	w = [w h; v d];                                          % 级联.
    end
    
    if nargout   == 0
    	imshow(w);                                               % 显示结果
    end
    

    示例:使用wave2gray显示变换系数

    img = imread('C:\Users\win\Desktop\sunflower.jpg');
    img = rgb2gray(img);  % wavefast用于处理灰度图像,故要进行图像转换
    [c,s] = wavefast(img,2,'haar');
    figure;wave2gray(c,s);title('小波分解系数(自动缩放)');
    figure;wave2gray(c,s,8);title('小波分解系数(细节8倍放大)');
    figure;wave2gray(c,s,-8);title('小波分解系数(细节8倍放大的绝对值)');
    

    结果:
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述

    3.6 快速小波反变换的实现

    将wfilters(小波工具箱自带函数)和wavefilter(非自带函数)中的参数type的值设置为“r”,即可得到重构滤波器(高通重构和低通重构)。
    在小波工具箱中,wavedec2表示二维尺度分解,得到小波分解结构[C,S]。其对应的waverec2二维小波重构,用来计算小波分解结构[C,S]的快速小波反变换。语法如下:

    g = waverec2(C,S,wname)   % g是重构后的二维图像(double类)
    g = waverec2(C,S,Lo_R,Hi_R)  
    

    当小波工具箱不可用时,可用waveback来实现二维小波重构。
    waveback.m

     function [varargout] = waveback(c, s, varargin)  
    %  waveback  多尺度二维快速小波反变换 
    %   [VARARGOUT] = WAVEBACK(C, S, VARARGIN) 计算分解结构的二维N尺度的部分或完全小波重构 
    %   
    %   Y = WAVEBACK(C, S, 'WNAME');     
    %   Y = WAVEBACK(C, S, LR, HR);   
    %   用低通和高通重构滤波器(LR和HR)或通过wavefilter的参数'wname'获得的滤波器,输出快速小波反变换矩阵Y
    %
    %   [NC, NS] = WAVEBACK(C, S, 'WNAME', N);     
    %   [NC, NS] = WAVEBACK(C, S, LR, HR, N);   
    %
    %   在N尺度重构后,输出新的小波分解结构
    %   See also WAVEFAST and WAVEFILTER.   
    
    %  检查输入输出参数的合理性 
    narginchk(3, 5);  
    narginchk(1, 2);  
    
    if (~ismatrix(c)) || (size(c, 1) ~= 1)  
      error('C必须是行向量.');     
    end  
    
    if (~ismatrix(s)) || ~isreal(s) || ~isnumeric(s) || (size(s,2) ~= 2)  
      error('S 必须是实数型的两列数组.');     
    end  
    
    elements = prod(s, 2);  
    if (length(c) < elements(end)) || ...  
          ~(elements(1) + 3 * sum(elements(2:end - 1)) >= elements(end))  
      error(['[C S] 必须是标准小波分解结构.']);   
    end  
    
    % [C, S]上的最大尺度.  
    nmax = size(s, 1) - 2;          
    % 获得第三个输入参数并初始化检查flags 
    wname = varargin{1};  filterchk = 0;   nchk = 0;      
    
    switch nargin  
    case 3  
       if ischar(wname)     
          [lp, hp] = wavefilter(wname, 'r');   n = nmax;  
       else   
          error('未定义滤波器.');    
       end  
       if nargout ~= 1   
          error('输出参数个数有误.');    
       end  
    case 4  
       if ischar(wname)  
          [lp, hp] = wavefilter(wname, 'r');     
          n = varargin{2};    nchk = 1;  
       else  
          lp = varargin{1};   hp = varargin{2};     
          filterchk = 1;   n = nmax;  
          if nargout ~= 1   
             error('输出参数个数有误.');    
          end  
       end  
    case 5  
        lp = varargin{1};   hp = varargin{2};   filterchk = 1;  
        n = varargin{3};    nchk = 1;  
    otherwise  
        error('输入参数个数不当.');       
    end  
    
    fl = length(lp);  
    if filterchk                                        % 检查滤波器.  
      if (~ismatrix(lp)) || ~isreal(lp) || ~isnumeric(lp) ...  
            || (~ismatrix(hp)) || ~isreal(hp) || ~isnumeric(hp) ...  
            || (fl ~= length(hp)) || rem(fl, 2) ~= 0  
          error('LP 和 HP必须是偶数并且是长度相等的实数滤波器向量.');   
      end  
    end  
    
    if nchk && (~isnumeric(n) || ~isreal(n))          % 检查 N 
        error('N 必须是实数.');   
    end  
    if (n > nmax) || (n < 1)  
       error('请求重建的数量(N)无效.');      
    end  
    if (n ~= nmax) && (nargout ~= 2)   
       error('输出参数不够.');   
    end  
    
    nc = c;    ns = s;    nnmax = nmax;             % 初始化分解 
    for i = 1:n  
       % 计算新的近似
       a = symconvup(wavecopy('a', nc, ns), lp, lp, fl, ns(3, :)) + ...  
           symconvup(wavecopy('h', nc, ns, nnmax), ...  
                     hp, lp, fl, ns(3, :)) + ...  
           symconvup(wavecopy('v', nc, ns, nnmax), ...  
                     lp, hp, fl, ns(3, :)) + ...  
           symconvup(wavecopy('d', nc, ns, nnmax), ...  
                     hp, hp, fl, ns(3, :));  
    
        % 更新分解 
        nc = nc(4 * prod(ns(1, :)) + 1:end);     nc = [a(:)' nc];  
        ns = ns(3:end, :);                       ns = [ns(1, :); ns];  
        nnmax = size(ns, 1) - 2;  
    end  
    
    % 对完全重建,将输出格式化为二维
    if nargout == 1  
        a = nc;   nc = zeros(ns(1, :));     nc(:) = a;      
    end  
    
    varargout{1} = nc;  
    if nargout == 2   
       varargout{2} = ns;    
    end  
    
    %-------------------------------------------------------------------%  
    function z = symconvup(x, f1, f2, fln, keep)  
    % 上采样行并用f1列卷积,上采样列并用f2行卷积,然后提取中心(假设对称扩展)
    
    y = zeros([2 1] .* size(x));      y(1:2:end, :) = x;  
    y = conv2(y, f1');  
    z = zeros([1 2] .* size(y));      z(:, 1:2:end) = y;  
    z = conv2(z, f2);  
    z = z(fln - 1:fln + keep(1) - 2, fln - 1:fln + keep(2) - 2);  
    

    示例:使用waveback进行图像的重构

    f = imread ( 'C:\Users\win\Desktop\sunflower.jpg' ) ;  
    f = rgb2gray(f);
    [C,S] = wavefast( f,4,'haar') ;  % 4层小波分解
    wave2gray(C,S,8); % 显示处理后的小波系数
    f = wavecopy ('a',C,S) ;  % f 是第4层分解后的近似系数
    figure; subplot(2,3,1);imshow(mat2gray(f)) ; title('4层分解后的近似图像');  % 显示第4层分解后的近似图像
    
    [C,S] = waveback(C,S,'haar',1) ;   % 一次重构
    f = wavecopy ('a',C,S) ;   %  f 是第3层分解后的近似图像
    subplot(2,3,2); imshow (mat2gray(f)) ; title('3层分解后的近似图像'); % 显示第3层分解后的近似图像
    
    [C,S] = waveback(C,S,'haar',1) ; % 二次重构
    f = wavecopy ('a',C,S) ;  %  f 是第2层分解后的近似系数
    subplot(2,3,3); imshow(mat2gray(f)) ; title('2层分解后的近似图像') % 显示第2层分解后的近似图像
    
    [C,S]= waveback(C, S,'haar',1);   %三次重构
    f = wavecopy('a',C,S);    %f是第1层分解后的近似函数
    subplot(2,3,4); imshow(mat2gray(f)); title('1层分解后的近似图像')%显示第1层分解后的近似图像
    
    [C, S]= waveback(C,S,'haar',1); %四次重构
    f = wavecopy('a',C,S);  %f是最终重构后的图像
    subplot(2,3,[5,6]); imshow(mat2gray(f)); title('重构后的最终图像')%显示重构后的最终图像
    

    结果
    在这里插入图片描述
    在这里插入图片描述
    wrcoef2 用二维小波系数重构某一层级的系数

    % 其实就是根据小波分解之后的系数c来重建其对应的图像。重建好的图像的尺度与原始图像一致。即无论你要重构哪个层的系数,最终它的维度都是和原始图像的尺度一致。
    X = wrcoef2('type',C,S,wname,N)  % 基于小波分解结构[C,S],计算N尺度上的重构系数
    X = wrcoef2('type',C,S,Lo_R,Hi_R,N)   % Lo_R,Hi_R可通过wfilters或wavefilter得到
    
    X = wrcoef2('type',C,S,wname)  % 未明示N的值时,此时计算最大尺度N上的重建系数
    X = wrcoef2('type',C,S,Lo_R,Hi_R)
    

    示例

    img  = imread('C:\Users\win\Desktop\sunflower.jpg');
    % 对图像进行3层的小波分解
    [c,s] = wavedec2(img,3,'db1');
    
    % 对各层的近似图像a进行重构
    a1 = wrcoef2('a',c,s,'db1',1);
    a2 = wrcoef2('a',c,s,'db1',2);
    a3 = wrcoef2('a',c,s,'db1',3);
    
    figure(1);
    subplot(2,2,1); imshow(img); title('原图');
    subplot(2,2,2); imshow(a1); title('1层重构近似系数');
    subplot(2,2,3); imshow(a2); title('2层重构近似系数');
    subplot(2,2,4); imshow(a3); title('3层重构近似系数');
    

    结果:
    在这里插入图片描述

    四、小波在图像处理中的应用

    基于小波的图像处理的基本方法如下:
    1、计算一幅图像的二维小波变换
    2、修改变换系数
    3、计算反变换

    小波的方向性和边缘检测

    f=imread('C:\Users\win\Desktop\sunflower.jpg');
    f = rgb2gray(f);
    figure;
    imshow(f);
    title('原图');
    [c,s]=wavefast(f,1,'sym4'); %快速小波变换
    figure;
    wave2gray(c,s,-6); 
    title('显示小波系数(细节放大6倍绝对值)');
    
    [nc,y]=wavecut('a',c,s);%近似系数置0
    figure;
    wave2gray(nc,s,-6);  
    title('显示小波系数(将近似系数置0,细节放大6倍绝对值)');
    edges=abs(waveback(nc,s,'sym4'));% 快速小波反变换,边缘重构
    
    figure;
    imshow(mat2gray(edges));
    title('图像的边缘')
    

    结果
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    基于小波的图像平滑
    为进行流线型的平滑处理,引入wavezero函数
    wavezero.m

    function [nc, g8] = wavezero(c, s, l, wname)  
    %  wavezero 将小波变换的细节系数置零  
    %   [NC, G8] = WAVEZERO(C, S, L, WNAME) 
    %   将L尺度上的小波分解结构[C,S]的细节系数置零,并计算WNAME小波的反变换
    
    [nc, foo] = wavecut('h', c, s, l);  
    [nc, foo] = wavecut('v', nc, s, l);  
    [nc, foo] = wavecut('d', nc, s, l);  
    i = waveback(nc, s, wname);  
    g8 = im2uint8(mat2gray(i));  
    figure; imshow(g8);  
    

    示例

    f=imread('C:\Users\win\Desktop\sunflower.jpg');
    f = rgb2gray(f);
    [c,s] = wavefast(f,3,'sym4');
    wave2gray(c,s,20);
    [c,g8] = wavezero(c,s,1,'sym4');  % 将第一级细节系数置为零
    [c,g8] = wavezero(c,s,2,'sym4');  % 将第二级细节系数置为零
    [c,g8] = wavezero(c,s,3,'sym4');  % 将第三级细节系数置为零
    

    结果
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    参考:小波的选择和小波的分类

    五、小波包

    由于正交小波变换只对信号的低频(近似)信息做进一步分解,而对高频(细节)信息不再继续分解,使得它的频率分辨率随频率升高而降低。所以小波变换能够很好地表征一大类以低频信息为主要成分的信号,但它不能很好地分解和表示包含大量细节信息(细小边缘或纹理)的信号,如非平稳机械振动信号、遥感图象、地震信号和生物医学信号等。与之不同的是,小波包变换可以对高频部分提供更精细的分解,而且这种分解既无冗余,也无疏漏,所以对包含大量中、高频信息的信号能够进行更好的时频局部化分析。
    一个三层的小波包分解如下:
    在这里插入图片描述
    式中 A 表示低频部分,D 表示高频部分,末尾的序号表示小波包分解的层数(即尺度数)。如果分解的是二维的图像,每一层的每个节点的分解系数就代表着一张图像,该图像亦或来自上一层近似系数的分解,亦或来自上一层细节系数的分解。

    小波包分解
    wpdec 一维小波包分解
    参考:wpdec 一维小波包分解

    T = wpdec(X,N,wname,E,P)   % 一维小波包分解函数,X为要分解的信号,N为分解级数(尺度数),wname为小波类型,E为熵的类型,默认的熵为‘shannon’,P为熵所对应的参数
    

    示例

    % 加载信号
    load noisdopp;  % noisdopp是小波工具箱自带的信号,直接加载即可
    x = noisdopp;
    wpt = wpdec(x,3,'db1','shannon');
    % 画对应的小波包树
    plot(wpt)
    

    结果
    在这里插入图片描述
    点击图中左侧的不同节点,右侧会显示对应节点的信号信息。

    wpdec2 二维小波包分解
    示例

    X = imread('C:\Users\win\Desktop\photo\peach.jpg');
    t = wpdec2(X,2,'db1'); 
    plot(t)
    

    结果
    在这里插入图片描述
    小波包重构
    wprec 一维重构
    示例

    load noisdopp;  % noisdopp是小波工具箱自带的信号,直接加载即可
    x = noisdopp;
    
    wpt = wpdec(x,3,'db1','shannon');
    % 一维重构
    y = wprec(wpt);
    plot(abs(x - y));
    title('重构后的信号和原信号的差值');
    

    结果
    在这里插入图片描述
    注意看纵坐标的单位是乘以10-15,所以差异还是比较小的。

    wprec2 二维重构
    示例

    X = imread('C:\Users\win\Desktop\photo\peach.jpg');
    t = wpdec2(X,2,'db1'); 
    Y = wprec2(t);
    figure;
    subplot(1,2,1); imshow(X); title('原图');
    subplot(1,2,2); imshow(Y); title('小波包重构后的图像');
    

    结果
    在这里插入图片描述
    获取小波包树上某个节点的小波包系数
    wpcoef 获取小波包树某节点的系数 (一维或二维)

    X = wpcoef(T,N)  % 返回小波包树某个节点N对应的小波包系数
    X = wpcoef(T)  % 相当于 X = wpcoef(T,0)  即获得原始信号节点(小波包树的根节点)的小波包系数 
    

    示例

    load noisdopp; 
    x = noisdopp;
    
    % 获得小波包树
    wpt = wpdec(x,3,'db1');
    
    % 获得[2,1]节点的系数
    cfs = wpcoef(wpt,[2 1]);
    
    % 显示
    figure(1); subplot(211); 
    plot(x); title('原始信号');
    
    subplot(212); 
    plot(cfs); title('节点[2,1]处的系数');
    

    在这里插入图片描述
    wprcoef 获取小波包树某节点的重构系数(一维或二维)

    X = wprcoef(T,N) %计算小波包树T的节点N的重构系数
    

    示例

    load noisdopp; x = noisdopp;
    t = wpdec(x,3,'db1','shannon');
    % 重构小波包结点(2,1)
    rcfs = wprcoef(t,[2 1]);
    figure(1); subplot(211);  plot(x); title('原始信号');
    subplot(212); plot(rcfs); title('重构小波包结点(2,1)');
    

    结果
    在这里插入图片描述
    将小波包树上的节点进一步分解
    wpsplt 分解小波包(一维或二维)

    T1 = wpsplt(T,N)   % T是小波包树,N为小波包树上的节点
    

    示例

    load noisdopp;
    x = noisdopp;
    
    wpt = wpdec(x,3,'db1');
    plot(wpt);
    wpt1 = wpsplt(wpt,9);  % 对小波包树wpt的第9个节点进行进一步分解,此处可以写9,也可以写[3,2]
    plot(wpt1);
    

    在这里插入图片描述
    在这里插入图片描述
    注: 小波包树的节点数是从0开始的,即树的根节点是第0个节点。

    将小波包树上的某些节点重组
    wpjoin 重组小波包某节点 (一维或二维)

    T = wpjoin(T,N)   % 将小波包树第N个节点以下的分支合并。N为0或者不给出N,都表示合并到根节点
    

    示例

    load noisdopp;
    x = noisdopp;
    
    wpt = wpdec(x,3,'db1');
    plot(wpt);
    wpt1 = wpjoin(wpt,[1,1]);
    plot(wpt1);
    

    在这里插入图片描述
    计算最优的小波包树
    besttree 获得最优小波包树(一维或二维)

    T = besttree(T)
    

    示例

    load noisdopp;
    x = noisdopp;
    
    wpt = wpdec(x,3,'db1');
    plot(wpt);
    wpt1 = besttree(wpt);
    plot(wpt1);
    

    结果
    在这里插入图片描述
    计算小波包树的最优层数
    beatlevt 计算最优层数(一维或二维)

    T = bestlevt(T)
    

    示例

    load noisdopp; 
    x = noisdopp;
    wpt = wpdec(x,3,'db1'); 
    
    wpt = wpsplt(wpt,[3 0]);
    plot(wpt);
    % 计算最优层数
    blt = bestlevt(wpt);
    plot(blt);
    

    结果
    在这里插入图片描述
    wpcutree 剪小波包树(一维或二维)

    T = wpcutree(T,L)  % 剪掉小波包树的第L层
    

    示例

    load noisdopp; 
    x = noisdopp;
    wpt = wpdec(x,3,'db1'); 
    plot(wpt);
    % 剪掉wpt的第二层
    newpt = wpcutree(wpt,2);
    plot(newpt);
    

    在这里插入图片描述
    wp2wtree 从小波包树提取小波树(一维或二维)

    T = wp2wtree(T)  % 计算对应于小波包分解树的修改后的小波包树T
    

    示例

    load noisdopp; x = noisdopp;
    wpt = wpdec(x,3,'db1');
    % 画出小波包树
    plot(wpt);
    % 提取小波树
    wt = wp2wtree(wpt);
    % 画出小波树
    plot(wt);
    

    结果
    在这里插入图片描述

    展开全文
  • 关于小波变换的一个入门级的介绍,对刚开始学习小波变换很有帮助。
  • 尽管时间和频率分辨率问题是物理现象(海森堡不确定性原理)的结果,并且无论使用哪种变换都存在,但是可以通过使用称为多分辨率分析(MRA)的替代方法来分析任何信号。 顾名思义,MRA可以分析具有不同分辨率的不同...

    原文

    1. 多分辨率分析

    尽管时间和频率分辨率问题是物理现象(海森堡不确定性原理)的结果,并且无论使用哪种变换都存在,但是可以通过使用称为多分辨率分析(MRA)的替代方法来分析任何信号。 顾名思义,MRA可以分析具有不同分辨率的不同频率的信号,但不能像STFT那样对每个频谱分量进行均等的解析。

    MRA被设计为在高频时具有良好的时间分辨率和较差的频率分辨率,在低频时具有良好的频率分辨率和较差的时间分辨率。 当手头信号在短时间内具有高频分量而在长时间内具有低频分量时,这种方法尤其有意义。 幸运的是,在实际应用中遇到的信号通常是这种类型的。 例如,下图显示了这种信号。 它在整个信号中具有相对较低的频率分量,并且在中间附近的短时间内具有较高的频率分量。

    在这里插入图片描述

    2. 连续小波变换

    连续小波变换被开发为短时傅立叶变换的替代方法,以克服分辨率问题。 小波分析以与STFT分析类似的方式进行,在某种意义上,信号乘以一个函数 t h e w a v e l e t {\it the wavelet} thewavelet,类似于STFT中的窗口函数,并且针对不同的信号分别计算变换时域信号的片段。 但是,STFT和CWT之间有两个主要区别:

    (1)没有对窗口信号进行傅立叶变换,因此将看到对应于正弦波的单个峰,即,不计算负频率。
    (2) 在为每个单个频谱分量计算变换时,窗口的宽度会发生变化,这可能是小波变换的最重要特征。

    连续小波变换定义如下:

    C W T x ψ ( τ , s ) = Ψ x ψ ( τ , s ) = 1 ∣ s ∣ ∫ x ( t ) ψ ∗ ( t − τ s ) d t CWT_x^\psi(\tau,s) = \Psi_x^\psi(\tau,s) = \frac{1}{\sqrt{|s|}} \int x(t) \psi^* \left( \frac{t - \tau}{s} \right) dt CWTxψ(τ,s)=Ψxψ(τ,s)=s 1x(t)ψ(stτ)dt

    如上式所示,变换后的信号是两个变量 τ \tau τ和s的函数,分别是" translation"和"scale"参数。 ψ ( t ) \psi(t) ψ(t)是变换函数,称为母小波( mother wavelet)。 术语“母小波”之所以得名,是因为小波分析的两个重要属性,如下所述:

    wavelet一词是指small wave。 small 是指此(窗口)函数具有有限长度(得到紧凑支持)的条件。wave是指该函数具有振荡性的条件。 术语“mother”表示在转换过程中使用的具有不同支持区域的函数是从一个主要函数或“母小波”派生的。 换句话说,母小波是用于生成其他窗口函数的原型。

    术语" translation"的含义与STFT中使用的含义相同; 它与窗口的位置有关,因为窗口在信号中移动。 显然,该术语对应于变换域中的时间信息。 但是,我们没有STFT之前的频率参数。 取而代之的是,我们将"scale"参数定义为 1 f r e q u e n c y \frac {1} {frequency} frequency1。 术语"frequency"是为STFT保留的。 下一部分将更详细地说明"scale"。

    3. The Scale

    小波分析中的参数“scale”似于地图中使用的“scale”。 与地图的情况一样,“high scale”对应于(信号的)非详细全局视图,而“low scale”对应于详细视图。 同样,就频率而言,低频(“high scale”)对应于信号的全局信息(通常跨越整个信号),而高频(“low scale”)对应于信号中隐藏模式的详细信息( 通常会持续相对较短的时间)。 下图给出了对应于各种刻度的余弦信号的示例。

    在这里插入图片描述
    幸运的是,在实际应用中,“ low scales”(高频)不会在信号的整个持续时间内持续存在,与图中所示不同,但它们通常会不时出现为短脉冲或尖峰。 “high scale”(低频)通常会持续整个信号持续时间。

    缩放(Scaling)作为一种数学运算,可以放大或压缩信号。 较大的"scale"对应于扩张(或扩展)的信号,较小的"scale"对应于压缩的信号。 图中给出的所有信号均来自相同的余弦信号,即它们是同一功能的扩张或压缩版本。 在上图中,s = 0.05是最小比例,而s = 1是最大比例。

    在数学函数方面,如果 f ( t ) f(t) f(t)是给定函数,则在s> 1的情况下, f ( s t ) f(st) f(st)对应于 f ( t ) f(t) f(t)的压缩(压缩)版本,而在s<1时,对应于f(t)的扩展(扩张)版本 。

    然而,在小波变换的定义中,在分母中使用了比例项,因此,上述陈述的反义成立,即,比例尺s> 1放大信号,而比例尺s <1压缩信号。 此处对“scale”的解释将在全文使用。

    4. CWT的计算

    在本节中将解释上述方程式。 令 x ( t ) x(t) x(t)是要分析的信号。 选择母小波作为该过程中所有窗口的原型。 使用的所有窗口都是母小波的放大(或压缩)版本和平移版本。 有许多用于此目的的功能。 Morlet小波和Mexican hat函数是两个候选者,它们用于示例的小波分析,这些示例将在本章后面介绍。

    一旦选择了母小波,计算就从s = 1开始,并对小于和大于“1”的所有s值计算连续小波变换。 但是,根据信号,通常不需要完全转换。对于所有实际目的,信号都是带宽受限的,因此,在有限的比例尺间隔内进行变换计算通常是足够的。 在本研究中,使用了s的一些有限值区间,这将在本章后面介绍。

    为了方便起见,该过程将从标度s = 1开始,并针对s的增加值继续进行,即,分析将从高频开始,朝低频进行。 s的第一个值将对应于压缩程度最高的小波。 随着s值的增加,小波将扩大。

    将小波放在信号的开始处与time= 0对应的点处。将小波函数以“1”的比例乘以信号,然后在所有时间进行积分。 然后,将积分结果乘以常数 1 s \frac {1} {\sqrt {s}} s 1。 该乘法是出于能量归一化的目的,因此变换后的信号在每个尺度上都将具有相同的能量。 最终结果是变换的值,即在时间零和标度s = 1处的连续小波变换的值。 换句话说,它是与时标平面中的 τ = 0 , s = 1 \boldsymbol\tau = 0,s = 1 τ=0s=1对应的值。

    然后将比例为s = 1的小波向右移 τ \tau τ到位置 t = τ t = \boldsymbol \tau t=τ,并计算上述等式以获得时间为 t = τ t=\tau t=τ,s = 1的变换值。 频率平面。

    重复该过程,直到小波到达信号的末端。 现在完成了时标平面上标度s = 1的一排点。

    然后,s增加一个小值。 请注意,这是一个连续的变换,因此 τ \boldsymbol \tau τ s s s都必须连续递增。 但是,如果此转换需要由计算机计算,则两个参数都将增加足够小的步长。 这对应于对时标平面进行采样。

    对s的每个值重复上述过程。 给定s值的每次计算都将填充时标平面的相应单行。 当针对s的所有期望值完成该过程时,已计算出信号的CWT。

    下图逐步说明了整个过程。

    在这里插入图片描述在图3.3中,显示了 τ \boldsymbol \tau τ四个不同值的信号和小波函数。 该信号是图3.1所示信号的截断版本。 “scale"值为1,对应于最低"scale"或最高频率。 请注意它有多紧凑(蓝色窗口)。 它应与信号中存在的最高频率分量一样窄。 小波函数的四个不同位置在图中显示为 t o = 2 , t o = 40 , t o = 90 和 t o = 140 \boldsymbol {t_o} = 2,\boldsymbol {t_o} = 40,\boldsymbol {t_o} = 90和\boldsymbol {t_o} = 140 to=2to=40to=90to=140。 在每个位置,它都乘以信号。 显然,乘积仅在信号落在小波支持区域内时才为非零,在其他地方为零。 通过在时间上移动小波,可以在时间上定位信号,通过更改s的值,可以在"sacle”(频率)上定位信号。

    如果信号的频谱分量对应于s的当前值(在这种情况下为1),则小波与该频谱分量存在的位置处的信号的乘积将给出一个较大的值。 如果信号中不存在与s的当前值相对应的频谱分量,则乘积值将相对较小或为零。 图3.3中的信号的频谱分量与t = 100 ms周围s = 1时窗口的宽度相当。

    图3.3中信号的连续小波变换将在100 ms左右的时间产生小范围的大数值,而在其他地方产生小数值。 另一方面,对于"high scale",连续小波变换将在信号的整个持续时间内给出较大的值,因为低频始终存在。

    在这里插入图片描述
    在这里插入图片描述
    图3.4和3.5分别说明了标度s = 5和s = 20的相同过程。 请注意窗口宽度如何随比例增加(频率减少)而变化。 随着窗口宽度的增加,变换开始拾取低频分量。

    结果,对于每个刻度和每个时间(间隔),将计算时标平面的一个点。 在一个尺度上的计算构成了时标平面的行,而在不同尺度上的计算构成了时标平面的列。

    现在,让我们看一个例子,看看小波变换的真实外观。 考虑图3.6中的非平稳信号。 除了在不同的频率外,这类似于为STFT给出的示例。 如图所示,信号由30 Hz,20 Hz,10 Hz和5 Hz的四个频率分量组成。
    在这里插入图片描述图3.7是该信号的连续小波变换(CWT)。 请注意,轴是“translation” 和“scale”,而不是时间和频率。 但是,“translation” 与时间严格相关,因为它指示了子小波所在的位置。 可以将母小波的“translation” 视为从t = 0开始经过的时间。 然而,“scale”却完全不同。 请记住,公式3.1中的“scale”参数s实际上是频率的倒数。 换句话说,无论我们说的关于频率分辨率的小波变换的性质如何,它的反面都会出现在表示时域信号WT的图中。

    在这里插入图片描述请注意,在图3.7中,较小的“scale”对应于较高的频率,即,频率随着“scale”的增加而减小,因此,图中“scale”在零附近的部分实际上对应于分析中的最高频率,而“scale”较高的那部分对应于最低的频率 频率。 请记住,信号首先具有30 Hz(最高频率)分量,然后以最低标度从0到30的转换出现。然后是20 Hz分量,第二最高频率,依此类推。 5 Hz分量出现在“translation” 轴的末端(如预期的那样),并且按预期再次出现在较高的音阶(较低的频率)上。

    在这里插入图片描述现在,回想一下这些分辨率属性:与STFT在所有时间和频率下均具有恒定的分辨率不同,WT在高频下具有良好的时间和较差的频率分辨率,而在低频下具有良好的频率和较差的时间分辨率。图3.8从另一个角度显示了图3.7中的相同WT,以更好地说明分辨率属性:在图3.8中,“low scale”(较高的频率)具有较好的scale分辨率(scale更窄,这意味着scale的确切值不那么模棱两可),这对应于较差的频率分辨率。同样,“high scale”具有scale频率分辨率,它对应于较低频率的更好的频率分辨率。

    图3.7和3.8中的轴已标准化,应进行相应评估。大致来说,平移轴上的100点对应于1000 ms,而刻度轴上的150点对应于40 Hz的频带(平移和刻度轴上的数字分别不对应于秒和Hz,它们只是计算中的样本数)。

    5. 时间和频率分辨率

    在本节中,我们将仔细研究小波变换的分辨率属性。 请记住,分辨率问题是我们从STFT切换到WT的主要原因。

    图3.9中的插图通常用于描述 如何解释时间和频率分辨率。 图3.9中的每个方框都对应于时频平面中的小波变换的值。 请注意,框具有某个非零区域,这意味着无法知道时频平面中特定点的值。 时频平面中落入一个框内的所有点都由WT的一个值表示。
    在这里插入图片描述
    让我们仔细看看图3.9:首先要注意的是,尽管盒子的宽度和高度会发生变化,但面积是恒定的。 也就是说,每个框代表时频平面的相等部分,但对时间和频率给出不同的比例。 请注意,在低频情况下,盒子的高度较短(对应于较好的频率分辨率,因为确切频率值的模棱两可性较小),但它们的宽度较长(这与较差的时间分辨率相对应,因为确切时间的值存在更多歧义)。 在较高的频率下,盒子的宽度减小,即时间分辨率变好,而盒子的高度增加,即频率分辨率变差。

    在结束本节之前,值得一提的是在STFT情况下分区是什么样子的。 回想一下,在STFT中,时间和频率分辨率由分析窗口的宽度决定,该宽度对于整个分析都是固定的,即时间和频率分辨率都是恒定的。 因此,在STFT情况下,时频平面由正方形组成。

    无论盒子的尺寸如何,STFT和WT中所有盒子的面积都是相同的,并由海森堡不等式确定。 总之,对于每个窗口函数(STFT)或母小波(CWT),框的面积是固定的,而不同的窗口或母小波会导致不同的面积。 但是,所有区域都由 1 4 π \boldsymbol {\frac {1} {4} \pi} 41π限制。 也就是说,由于海森堡的不确定性原理,我们无法尽可能地减少盒子的面积。 另一方面,对于给定的母小波,可以更改盒子的尺寸,同时保持面积不变。 这正是小波变换所做的。

    6. 小波理论:一种数学方法

    本节介绍了小波分析理论的主要思想,也可以将其视为大多数信号分析技术的基本概念。傅里叶定义的FT使用基本函数来分析和重构函数。向量空间中的每个向量都可以写成该向量空间中基本向量的线性组合,即,将向量乘以某个常数,然后取乘积之和。信号分析涉及对这些常数的估计(变换系数或傅立叶系数,小波系数等)。合成或重构对应于计算线性组合方程。

    与该主题相关的所有定义和定理都可以在Keiser的书《小波友好指南》中找到,但是对于基础函数如何工作的入门级知识对于理解小波理论的基本原理是必需的。因此,此信息将在本节中介绍。

    7. 基础向量

    向量空间V的基是一组线性独立的向量,因此V中的任何向量v都可以写成这些基向量的线性组合。 向量空间可能有多个基础。 但是,它们全部具有相同数量的向量,该数目称为向量空间的维数。 例如,在二维空间中,基础将具有两个向量。

    v = ∑ k ν k b k v = \sum\limits_{k} \nu^k b_k v=kνkbk·········································································(3.2)

    公式3.2展示了如何将任何向量v表示为基向量 b k \boldsymbol {b_k} bk和相应系数 ν k \boldsymbol {\nu ^ k} νk的线性组合。

    通过用基函数 ϕ k ( t ) \boldsymbol {\phi_k(t)} ϕk(t)替换基向量 b k \boldsymbol {b_k} bk,并用函数 f ( t ) f(t) f(t)替换向量 v v v,可以很容易地将这种就矢量而言的概念推广为函数 。 公式3.2变为

    f ( t ) = ∑ k μ k ϕ k ( t ) f(t) = \sum\limits_{k} \mu_k \phi_k (t) f(t)=kμkϕk(t)·······························································(3.2-a)

    复指数函数(正弦和余弦)是FT的基本函数。 此外,它们是正交函数,为重构提供了一些理想的属性。

    f ( t ) f(t) f(t) g ( t ) g(t) g(t) L 2 [ a , b ] L ^ 2 [a,b] L2[ab]中的两个函数。 ( L 2 [ a , b ] L ^ 2 [a,b] L2[ab]表示区间[a,b]中的平方可积函数的集合)。 两个函数的内积由公式3.3定义:

    < f ( t ) , g ( t ) > = ∫ a b f ( t ) ⋅ g ∗ ( t ) d t < f(t), g(t) > = \int_a^b f(t) \cdot g^*(t) dt <f(t),g(t)>=abf(t)g(t)dt········································(3.3)

    根据内积的上述定义,可以将CWT视为具有基函数 ψ ( τ , s ) ( t ) \psi_{(\tau ,s)}(t) ψ(τ,s)(t)的测试信号的内积:

    C W T x ψ ( τ , s ) = Ψ x ψ ( τ , s ) = ∫ x ( t ) ⋅ ψ τ , s ∗ ( t ) d t CWT_x^\psi(\tau, s) = \Psi_x^\psi(\tau, s) = \int x(t) \cdot \psi^*_{\tau, s}(t) dt CWTxψ(τ,s)=Ψxψ(τ,s)=x(t)ψτ,s(t)dt ·····················(3.4)

    其中:

    ψ τ , s = 1 s ψ ( t − τ s ) \psi_{\tau, s} = \frac{1}{\sqrt{s}} \psi \left( \frac{t - \tau}{s} \right) ψτ,s=s 1ψ(stτ)······························································(3.5)

    CWT的这一定义表明,小波分析是对基函数(小波)和信号本身之间相似度的度量。 在此,相似性是指相似的频率含量。 计算出的CWT系数是指在当前尺度下信号与小波的接近程度。

    这进一步澄清了先前关于信号与小波在一定规模上的相关性的讨论。 如果信号具有对应于当前标度的频率的主要分量,则当前标度的小波(基函数)将与该频率分量出现的特定位置处的信号相似或接近。 因此,在time-scale平面上此时计算的CWT系数将是一个相对较大的数字。

    8. 内积,正交性和正交性

    如果两个向量v,w的内积等于零,则称它们正交

    < v , w > = ∑ n v n w n ∗ = 0 < v, w > = \sum\limits_{n} v_n w^*_n = 0 <v,w>=nvnwn=0··························································· ······(3.6)

    类似地,如果两个函数f和g的内积为零,则称它们彼此正交:

    < f ( t ) , g ( t ) > = ∫ a b f ( t ) ⋅ g ∗ ( t ) ⋅ d t = 0 < f(t), g(t) > = \int_a^b f(t) \cdot g^*(t) \cdot dt = 0 <f(t),g(t)>=abf(t)g(t)dt=0············································(3.7)

    如果向量 v 1 , v 2 , . . . . , v n {\boldsymbol {v_1,v_2,....,v_n}} v1v2....vn的集合是成对正交的,则它们成对相互正交,并且所有长度均为“ 1”。 可以表示为:

    < v m , v n > = δ m n < v_m, v_n > = \delta_{mn} <vm,vn>=δmn············································································(3.8)

    同样,如果满足以下条件,则一组函数 ϕ k ( t ) , k = 1 , 2 , 3 , . . . {\phi_k(t)},k = 1,2,3,... ϕk(t)k=1,2,3...是正交的:

    ∫ a b ϕ k ( t ) ⋅ ϕ l ∗ ( t ) ⋅ d t = 0 , k ≠ l \int_a^b \phi_k(t) \cdot \phi^*_l(t) \cdot dt = 0,k \neq l abϕk(t)ϕl(t)dt=0,k=l(正交条件)···········································(3.9)

    ∫ a b { ∣ ϕ k ( t ) ∣ } 2 d x = 1 \int_a^b \{ | \phi_k(t) | \}^2 dx = 1 ab{ϕk(t)}2dx=1········································································(3.10)

    或同等:

    ∫ a b ϕ k ( t ) ⋅ ϕ l ∗ ( t ) ⋅ d t = δ k l \int_a^b \phi_k(t) \cdot \phi_l^*(t) \cdot dt = \delta_{kl} abϕk(t)ϕl(t)dt=δkl································································(3.11)

    其中, δ k l \delta_{kl} δkl是Kronecker delta函数,定义为:

    δ k l = { 1 , k = l 0 , k ≠ l \delta_{kl} = \left\{ \begin{array}{ll} 1, & k = l \\ 0, & k \neq l\\ \end{array} \right. δkl={1,0,k=lk=l······································································(3.12)

    如上所述,可能有不止一组基函数(或向量)。 其中,正交基函数(或向量)特别重要,因为它们在查找这些分析系数时提供了很好的特性。 正交正态基数允许使用正交正态性以非常简单明了的方式计算这些系数。

    对于正交基,系数 μ k \mu_k μk可以计算为

    μ k = < f , ϕ k > = ∫ f ( t ) ⋅ ϕ k ∗ ( t ) ⋅ d t \mu_k = < f, \phi_k > = \int f(t) \cdot \phi_k^*(t) \cdot dt μk=<f,ϕk>=f(t)ϕk(t)dt···············································(3.13)

    然后可以用公式3.2_a代入 μ k \mu_k μk系数来重建函数 f ( t ) f(t) f(t)。 这产生:

    f ( t ) = ∑ k μ k ϕ k ( t ) = ∑ k < f , ϕ k > ϕ k ( t ) f(t) = \sum\limits_{k} \mu_k \phi_k(t) = \sum\nolimits_{k} < f, \phi_k > \phi_k(t) f(t)=kμkϕk(t)=k<f,ϕk>ϕk(t)····································(3.14)

    正交基可能不适用于可以使用通用版本双正交基的每种类型的应用程序。 术语“双正交”是指彼此正交的两个不同的基,但是每个都不形成正交集。

    但是,在某些应用中,在这种情况下,可以使用双正交基。 框架是小波理论的重要组成部分,感兴趣的读者可以参考前面提到的Kaiser的书。

    遵循与第2章中STFT相同的顺序,接下来介绍连续小波变换的一些示例。 示例中给出的数字是由编写用于计算CWT的程序生成的。

    在结束本节之前,我想介绍两个在小波分析中常用的小波。 墨西哥帽小波定义为高斯函数的二阶导数:

    w ( t ) = 1 2 π ⋅ σ e − t 2 2 σ 2 w(t) = \frac{1}{\sqrt{2\pi} \cdot \sigma} e^{\frac{-t^2}{2 \sigma^2}} w(t)=2π σ1e2σ2t2········································································(3.15)

    ψ ( t ) = 1 2 π ⋅ σ 3 ( e − t 2 2 σ 2 ⋅ ( t 2 σ 2 − 1 ) ) \psi(t) = \frac{1}{\sqrt{2 \pi} \cdot \sigma^3} \left( e^{\frac{-t^2}{2 \sigma^2}} \cdot \left( \frac{t^2}{\sigma^2} - 1 \right) \right) ψ(t)=2π σ31(e2σ2t2(σ2t21))·················································(3.16)

    Morlet小波定义为:

    w ( t ) = e i a t ⋅ e − t 2 2 σ w(t) = e^{i a t} \cdot e^{-\frac{t^2}{2\sigma}} w(t)=eiate2σt2······································································(3.16-a)

    其中a是调制参数, σ \sigma σ是影响窗口宽度的缩放参数。

    9. 例子

    下面给出的所有示例均对应于现实生活中的非平稳信号。 这些信号来自数据库信号,其中包括正常人和患有阿尔茨海默氏病患者的事件相关电位。 由于这些不是简单的正弦波那样的测试信号,因此很难解释它们。 此处显示它们只是为了让您了解现实中的CWT到底是什么样子的。

    图3.11所示的以下信号属于正常人。

    在这里插入图片描述以下是其CWT。 轴上的数字对我们并不重要。 这些数字只是表明CWT是在translation-scale平面上的350个translation和60个scale位置计算的。 这里要注意的重要一点是,该计算不是真正的连续WT,这从有限数量的位置处的计算可以明显看出。 这只是CWT的离散版本,将在本页后面解释。 但是请注意,这不是离散小波变换(DWT),这是本教程第四部分的主题。

    在这里插入图片描述图3.13从不同角度绘制了相同的变换,以实现更好的可视化。

    在这里插入图片描述图3.14绘制了被诊断患有阿尔茨海默氏病的患者的事件相关电位

    在这里插入图片描述图3.15说明了它的CWT:
    在这里插入图片描述这是从另一个角度来看的另一种观点:
    在这里插入图片描述

    10.小波综合

    只要满足公式3.18,连续小波变换就是可逆变换。 幸运的是,这是一个非限制性的要求。 即使基本函数通常可能不是正交的,如果满足公式3.18,连续小波变换也是可逆的。 通过使用以下重建公式,可以进行重建:

    x ( t ) = 1 C ψ 2 ∫ s ∫ τ [ Ψ x ψ ( τ , s ) 1 s 2 ψ ( t − τ s ) ] d τ ⋅ d s x(t) = \frac{1}{C_\psi^2} \int_s \int_\tau \left[ \Psi^\psi_x(\tau, s) \frac{1}{s^2} \psi \left( \frac{t - \tau}{s} \right) \right] d\tau \cdot ds x(t)=Cψ21sτ[Ψxψ(τ,s)s21ψ(stτ)]dτds··················(3.17)

    其中 C ψ C_ \psi Cψ是一个常数,取决于所使用的小波。 重构的成功取决于此常数(可接纳性常数)是否满足以下可接纳性条件:

    C ψ = { 2 π ∫ − ∞ ∞ ∣ ψ ^ ( ξ ) ∣ 2 ∣ ζ ∣ d ξ } 1 2 < ∞ C_\psi = \left\{ 2 \pi \int_{-\infty}^{\infty} \frac{|\hat{\psi}(\xi)|^2}{|\zeta|} d\xi \right\} ^{\frac{1}{2}} < \infty Cψ={2πζψ^(ξ)2dξ}21<·····································(3.18)

    其中 ψ ^ ( ξ ) \hat{\psi}(\xi) ψ^(ξ) ψ ( t ) \psi(t) ψ(t)的FT。 公式3.18表示 ψ ^ ( 0 ) = 0 \hat {\psi}(0)= 0 ψ^(0)=0,这是

    ∫ ψ ( t ) ⋅ d t = 0 \int \psi(t) \cdot dt = 0 ψ(t)dt=0·······························································(3.19)

    如上所述,公式3.19并不是一个严格的要求,因为可以找到许多积分为零的小波函数。 为了满足公式3.19,小波必须是振荡的。

    11.连续小波变换的离散化:小波级数

    在当今世界,计算机用于执行大多数计算(嗯,…好吧…几乎所有计算)。 显然,FT,STFT或CWT都无法通过使用解析方程,积分等来实际计算。因此,有必要离散化变换。 就像在FT和STFT中一样,最直观的方法是简单地对time-frequency (scale)平面进行采样。 更直观地讲,以均匀的采样率对time-frequency (scale)平面采样听起来是最自然的选择。 但是,在WT的情况下,可以使用缩放比例降低采样率。

    根据Nyquist的规则,在"high scale"(较低的频率)下,可以降低采样率。 换句话说,如果需要"scale" s 1 \boldsymbol {s_1} s1 N 1 \boldsymbol {N_1} N1的采样率对time-scale平面进行采样,则可以在以下时间以 N 2 \boldsymbol {N_2} N2的采样率对同一平面进行采样 缩放 s 2 \boldsymbol {s_2} s2,其中 s 1 < s 2 \boldsymbol {s_1 <s_2} s1<s2(对应于频率 f 1 > f 2 \boldsymbol {f_1> f_2} f1>f2)和 N 2 < N 1 \boldsymbol {N_2 <N_1} N2<N1 N 1 \boldsymbol {N_1} N1 N 2 \boldsymbol {N_2} N2之间的实际关系为:

    N 2 = s 1 s 2 N 1 N_2 = \frac{s_1}{s_2} N_1 N2=s2s1N1···········································································(3.20)

    N 2 = f 2 f 1 N 1 N_2 = \frac{f_2}{f_1} N_1 N2=f1f2N1···········································································(3.21)

    换句话说,在较低的频率下,可以降低采样率,这将节省大量的计算时间。

    然而,此时应注意,就信号分析而言,离散化可以任何方式进行而没有任何限制。 如果不需要合成,则甚至不需要满足奈奎斯特标准。 当且仅当需要信号重构时,对离散化和采样率的限制才变得重要。 奈奎斯特的采样率是允许从其离散采样中重建原始连续时间信号的最小采样率。 由于这个原因,前面提到的基向量特别重要。

    如前所述,满足方程式3.18的小波 ψ ( τ , s ) \boldsymbol {\psi(\tau,s)} ψτs允许通过方程式3.17重构信号。 但是,对于连续变换而言确实如此。 问题是:如果离散化时间和比例参数,我们还能重构信号吗? 在某些条件下,答案是“是”(正如他们在商业广告中经常说的:有一定的限制!!!)。

    比例参数s首先在对数网格上离散。 然后将时间参数相对于比例参数离散化,即每个比例使用不同的采样率。 换句话说,采样是在图3.17所示的二元采样网格上完成的:

    在这里插入图片描述将轴所覆盖的区域视为整个时间刻度平面。 CWT为该平面上的连续点分配一个值。因此,存在无限数量的CWT系数。首先考虑scale轴的离散化。在对数无限的点中,使用对数规则仅取有限数。对数的底数取决于用户。由于其便利性,最常见的值为2。如果选择2,则仅刻度2、4、8、16、32、64等被计算。如果值为3,则刻度为3、9、27、81、243等可以被算出来。然后根据比例轴的离散化来离散时间轴。由于离散比例的变化是2倍,因此在每个比例上,时间轴的采样率都会降低2倍。

    注意,在最小刻度(s = 2)上,仅采样了时间轴的32个点(对于图3.17中给出的特殊情况)。在下一个标度值s = 4处,由于标度增加了2倍,因此时间轴的采样率降低了2倍,因此仅采集了16个样本。下一步,s = 8,并及时获取8个样本,依此类推。

    尽管它被称为time-scale平面,但将其称为translation-scale平面更为精确,因为变换域中的“时间”实际上对应于小波在时间上的移动。对于小波系列,实际时间仍然是连续的。

    与连续傅立叶变换,傅立叶级数和离散傅立叶变换之间的关系相似,存在连续小波变换,半离散小波变换(也称为小波级数)和离散小波变换。

    用数学术语表示上述离散化过程,scale离散化为 s = s 0 j \boldsymbol {s = s_0 ^ j} s=s0j,translation离散化为 τ = k ⋅ s 0 j ⋅ τ 0 \boldsymbol {\tau = k \cdot s_0 ^ j \cdot \tau_0} τ=ks0jτ0其中 s 0 > 1 \boldsymbol {s_0> 1} s0>1 τ 0 > 0 \boldsymbol {\tau_0> 0} τ0>0。 请注意,translation离散化是如何依赖于带有 s 0 \boldsymbol {s_0} s0的小数位数离散化的。

    连续小波函数:

    ψ τ , s = 1 s ψ ( t − τ s ) \psi_{\tau, s} = \frac{1}{\sqrt{s}} \psi \left( \frac{t - \tau}{s} \right) ψτ,s=s 1ψ(stτ)································································(3.22)

    ψ j , k ( t ) = s 0 − j 2 ψ ( s 0 − j − k τ 0 ) \psi_{j, k}(t) = s_0^{\frac{-j}{2}} \psi \left( s_0^{-j} - k \tau_0 \right) ψj,k(t)=s02jψ(s0jkτ0)·················································(3.23)

    通过插入 s = s 0   j \boldsymbol{s = s_0^{\, j}} s=s0j τ = k ⋅ s 0   j ⋅ τ 0 \boldsymbol{\tau = k \cdot s_0^{\, j} \cdot \tau_0} τ=ks0jτ0

    如果 { ψ ( j ,   k ) } \boldsymbol{ \left\{ \psi_{(j, \, k)} \right\} } {ψ(j,k)}构成正交基,则小波级数变换变为

    Ψ x ψ   j , k = ∫ x ( t )   ψ j ,   k ∗ ( t )   d t \Psi^{\psi_{\, j,k}}_x = \int x(t) \, \psi^*_{j, \, k}(t) \, dt Ψxψj,k=x(t)ψj,k(t)dt······················································(3.24)

    或者:

    x ( t ) = c ψ ∑ j ∑ k Ψ x ψ   j , k   ψ   j , k ( t ) x(t) = c_\psi \sum\limits_{j} \sum\limits_{k} \Psi^{\psi_{\, j,k}}_x \, \psi_{\, j,k} (t) x(t)=cψjkΨxψj,kψj,k(t)············································`···(3.25)

    小波级数要求 ψ ( j ,   k ) \boldsymbol{ {\psi_{(j, \, k)}} } ψ(j,k)可以是正交,双正交或框架。 如果 ψ ( j , k ) \boldsymbol{ {\psi_{(j,k)}} } ψ(j,k)不正交,则公式3.24变为:

    Ψ x ψ   j , k = ∫ x ( t )   ψ ( j ,   k ) ∗ ( t ) ^   d t \Psi^{\psi_{\, j,k}}_x = \int x(t) \, \hat{\psi^*_{(j, \, k)}(t)} \, dt Ψxψj,k=x(t)ψ(j,k)(t)^dt····················································(3.26)

    其中 ψ j ,   k ∗ ( t ) ^ \boldsymbol{\hat{\psi_{j, \, k}^*(t)}} ψj,k(t)^是双双正交基或双帧(请注意,*表示共轭)。

    如果 { ψ ( j ,   k ) } \boldsymbol{ \{ \psi_{(j, \, k)} \} } {ψ(j,k)}是正交或双正交的,则该变换将是非冗余的,就好像它们形成一个帧一样,该变换将是冗余的。 另一方面,查找框架比查找正交或双正交底物容易得多。

    以下类推可能会清除此概念。 将整个过程视为查看特定对象。 人眼首先确定粗略视图,该粗略视图取决于眼睛到物体的距离。 这对应于调整比例参数 s 0 − j \boldsymbol {s_0 ^ {-j}} s0j。 当观看非常接近的物体时,其细节非常丰富,j为负且较大(低比例,高频,分析信号中的细节)。 非常缓慢地移动头部(或眼睛),并以很小的增量(根据所观察的对象,其角度,距离)相应于 τ = k ⋅ s 0   j ⋅ τ 0 \boldsymbol{\tau = k \cdot s_0^{\, j} \cdot \tau_0} τ=ks0jτ0。 请注意,当j为负且较大时,它对应于时间 τ \boldsymbol {\tau} τ的小变化(高采样率)和 s 0   j \boldsymbol{s_0^{\, j}} s0j的大变化(小比例,高频 ,其中采样率很高)。 比例参数也可以认为是放大倍数。

    采样率有多低,仍然允许信号重构? 这是优化程序要回答的主要问题。 发现最方便的值(就编程而言)对于 s 0 s_0 s0为“ 2”,对于 τ \tau τ为“ 1”。 显然,当采样率被迫尽可能低时,可用正交小波的数量也会减少。

    本章中给出的连续小波变换示例实际上是给定信号的小波序列。 根据信号选择参数。 由于不需要重建,因此对于不同的示例,采样率有时远低于临界值,其中 s 0 s_0 s0从2变为10, τ 0 \tau_0 τ0从2变为8。

    到此结束本教程的第三部分。 希望您现在对小波变换的全部内容有一个基本的了解。 但是,还有一件事需要讨论。 即使离散小波变换可以在计算机上进行计算,根据信号大小和所需的分辨率,此计算可能需要几秒钟到几小时不等。 实际上,一种惊人的快速算法可用于计算信号的小波变换。 离散小波变换(DWT)在本教程的最后一章的第四部分中介绍。

    展开全文
  • 在双正交情况下,有两个标度函数\phi,\tilde\phi,可能产生不同的多分辨率分析,相应地有两个不同的小波函数\psi,\tilde\psi。 因此,缩放序列 a,\tilde a 中系数的数量 M 和 N 可能不同。 缩放序列必须满足以下双...
  • 很好的小波及多分辨率分析的入门和提高资料,深入浅出的全面分析了小波的各种基本理论,并且用恰当的实例和比喻对分辨率进行分析和解释,有利于初学者学习
  • 小波分析2006(第4讲:离散小波变换多分辨率分析)1.pdf
  • 小波变换多分辨率分析).ppt 电子科技大学 图像处理 课件 简单
  • 数字图像处理与Python实现笔记摘要绪论1 数字图像处理基础知识2 彩色图像处理初步3 空间滤波4 频域滤波5 图像特征提取6 图像压缩7 图像小波变换多分辨率 摘要 简要介绍数字图像处理涉及的一些基本概念、基本运算...

    摘要

    1. 简要介绍数字图像处理涉及的一些基本概念、基本运算、基本类型,以及如何通过Python对数字图像进行读取和简单操作。
    2. 以彩色图像为例对数字图像处理的基本操作进行介绍,熟悉数字图像处理的基本过程,主要包括颜色空间的基本概念、伪彩色图像处理操作,彩色图像处理简单操作。
    3. 瞄准在空间域中对图像进行增强,介绍空间滤波的机理、基本概念以及使用的基本技术。本章内容包括空间滤波基本概念、基于空间滤波的图像平滑处理、基于空间滤波的锐化操作以及混合空间增强。
    4. 从频域角度入手对图像处理及增强方法展开介绍。因为频域滤波所需的数学知识较多,所以本章采取由浅入深的策略,首先介绍一维傅里叶变换,其次介绍二维傅里叶变换和快速傅里叶变换,最后介绍图像频域滤波中出现的各种技术,其大体可分为低通滤波和高通滤波两大类。
    5. 从全局特征提取和局部特征提取两方面入手,分别介绍颜色特征、纹理特征、形状特征、边缘特征、点特征的提取方法。本章内容是目前机器视觉和图像处理领域的学者关注较多的内容,通过穿插较多的实例,帮助读者理解图像特征提取的基本技术。
    6. 瞄准如何减少图像传输及存储数据大小,介绍主要使用的压缩技术,包括有损压缩和无损压缩等,并使用JPEG压缩技术串讲全章知识点。
    7. 介绍图像的小波域表示及多分辨率表示。

    绪论

    • 人工智能是引领未来发展的战略性技术,是新一轮科技革命和产业变革的重要驱动力量,将深刻地改变人类社会生活。

    • 促进人工智能和实体经济的深度融合,构建数据驱动、人机协同、跨界融合、共创分享的智能经济形态,更是推动质量变革、效率变革、动力变革的重要途经。

    • 进年来,我国人工智能新技术、新产品、新业态持续涌现,与农业、制造业、服务业等行业的融合步伐明显加快,在技术创新、应用推广、产业发展等方面成效初显。

    • 人工智能技术并不是一个新生事物,它在最近几年引起全球性关注并得到飞速发展的主要原因,在于它的三个基本要素(算法、数据、算力)的迅猛发展,其中又以数据和算力的发展尤为重要。

    • 物联网技术的蓬勃发展使得数据累计的难度越来越低,而芯片算力的不断提升,使得过去只能通过云计算才能完成的人工智能运算,现在可以下沉到最普通的设备上完成。

    • 物联网技术为机器带来感知能力,而人工智能则通过计算算力为机器带来了决策能力,正如感知和大脑对自然生命进化所起到的必然性作用。

    1 数字图像处理基础知识

    https://hulin.blog.csdn.net/article/details/107570020

    2 彩色图像处理初步

    https://hulin.blog.csdn.net/article/details/107578369

    3 空间滤波

    https://hulin.blog.csdn.net/article/details/107589248

    4 频域滤波

    https://hulin.blog.csdn.net/article/details/107609844

    5 图像特征提取

    https://hulin.blog.csdn.net/article/details/107639032

    6 图像压缩

    https://hulin.blog.csdn.net/article/details/107693170

    7 图像小波变换与多分辨率

    • 小波变换是近年来图像处理中十分受重视的新技术,面向图像压缩、特征检测、纹理分析等提出了很多新方法,如多分辨率分析、时频域分析、金字塔算法等,都属于小波变换的范畴。
    • 信号分析是为了时间和频率之间的相互关系。傅里叶变换提供了有关频域的信息,但有关时间的局部化信息却基本丢失。与傅里叶变换不同,小波变换是通过缩放母小波(Mother Wavelet)的宽度获得信号的频域特征,通过平移母小波获得信号的时间信息。对母小波的缩放和平移是为了计算小波系数,这些小波系数反映了小波和局部信号之间的相关程度。
    • 像傅里叶分析一样,小波分析就是把一个信号分解为将母小波经过缩放和平移之后的一系列小波,因此小波是小波变换的基函数。小波变换可以理解为用经过缩放和平移的一系列小波函数,代替傅里叶变换的正弦波和余弦波进行傅里叶变换的结果。
    • 小波变换中的是指在时域具有紧支集或近似紧支集,是指具有正负交替的波动性,直流分量为0。小波本质上是定义在有限间隔而且其平均值为0的一种函数。
    • 与傅里叶变换相比,小波变换是空间(时间)和频率的局部变换,通过伸缩平移运算,对信号逐步进行多尺度细化,最终达到高频处时间细分,低频处频率细分,能自动适应时频信号分析的要求,从而可聚焦到信号的任意细节。
    • 小波变换是基于具有变化的频率和有限持续时间的小型波进行的。

    7.1 从傅里叶变换到小波变换

    7.1.1 小波

    1. 小波的概念

    • 小波是在有限时间范围内变化且其平均值为0的数学函数。具有两个特点。
      (1)具有有限的持续时间和突变的频率和振幅。
      (2)在有限的时间范围内,它的平均值为0
    • 小波变换的结果为各种小波系数,这些系数由尺度和位移函数组成。

    2. 小波变换

    • 通过小波对一个信号在空间和时间上进行局部化的一种数学变换,通过平移母小波,捕获到信号信息。通过缩放母小波的宽度捕获到信号的频域特性。对母小波的平移和缩放操作是为计算小波分量的系数,这些系数代表局部信号和小波之间的相互关系,这些参数反映了信号的时间属性和频率属性。

    7.1.2 感性认识小波变换

    • 傅里叶变换一直是信号处理领域应用最广泛、效果最好的一种分析手段,是时域到频域互相转化的工具。从物理意义上,傅里叶变换的实质是把对原函数的研究转化为对其傅里叶变换的研究。但是,傅里叶变换只能提供信号在整个时域上的频率,不能提供信号在某个局部时间段上的频率信息。
    • 傅里叶变换:在时域的常量函数,在频域将表现为冲击函数,表明具有很好的频域局部化性质。
      在这里插入图片描述
    import numpy as np
    import matplotlib.pyplot as plt
    from scipy.fftpack import fft
    
    
    # 中文显示工具函数
    def set_ch():
        from pylab import mpl
        mpl.rcParams['font.sans-serif'] = ['FangSong']
        mpl.rcParams['axes.unicode_minus'] = False
    
    
    set_ch()
    t = np.linspace(0, 1, 400, endpoint=False)
    cond = [t < 0.25, (t >= 0.25) & (t < 0.5), t >= 0.5]
    f1 = lambda t: np.cos(2 * np.pi * 10 * t)
    f2 = lambda t: np.cos(2 * np.pi * 50 * t)
    f3 = lambda t: np.cos(2 * np.pi * 100 * t)
    y1 = np.piecewise(t, cond, [f1, f2, f3])
    y2 = np.piecewise(t, cond, [f2, f1, f3])
    Y1 = abs(fft(y1))
    Y2 = abs(fft(y2))
    
    plt.figure(figsize=(12, 9))
    
    plt.subplot(221)
    plt.plot(t, y1)
    plt.title('信号1 时间域')
    plt.xlabel('时间/s')
    
    plt.subplot(222)
    plt.plot(range(400), Y1)
    plt.title('信号1 频率域')
    plt.xlabel('频率/Hz')
    
    plt.subplot(223)
    plt.plot(t, y2)
    plt.title('信号2 时间域')
    plt.xlabel('时间/s')
    
    plt.subplot(224)
    plt.plot(range(400), Y2)
    plt.title('信号2 频率域')
    plt.xlabel('频率/Hz')
    
    plt.show()
    
    
    • 从时域上看,相差很大的两个信号,在频域上却非常相近。一个很自然的方法是加窗(短时距傅里叶变换),将长时间信号分成数个较短的等长信号,然后再分别对每个窗进行傅里叶变换,从而得到频率随时间的变化,这就是短时距傅里叶变换
      在这里插入图片描述
    • 小波变换可以解决这个问题。
    import numpy as np
    import matplotlib.pyplot as plt
    import pywt
    
    
    # 中文显示工具函数
    def set_ch():
        from pylab import mpl
        mpl.rcParams['font.sans-serif'] = ['FangSong']
        mpl.rcParams['axes.unicode_minus'] = False
    
    
    set_ch()
    t = np.linspace(0, 1, 400, endpoint=False)
    cond = [t < 0.25, (t >= 0.25) & (t < 0.5), t >= 0.5]
    
    f1 = lambda t: np.cos(2 * np.pi * 10 * t)
    f2 = lambda t: np.cos(2 * np.pi * 50 * t)
    f3 = lambda t: np.cos(2 * np.pi * 100 * t)
    
    y1 = np.piecewise(t, cond, [f1, f2, f3])
    y2 = np.piecewise(t, cond, [f2, f1, f3])
    
    cwtmatr1, freqs1 = pywt.cwt(y1, np.arange(1, 200), 'cgau8', 1 / 400)
    cwtmatr2, freqs2 = pywt.cwt(y2, np.arange(1, 200), 'cgau8', 1 / 400)
    
    plt.figure(figsize=(12, 9))
    
    plt.subplot(221)
    plt.plot(t, y1)
    plt.title('信号1 时间域')
    plt.xlabel('时间/s')
    
    plt.subplot(222)
    plt.contourf(t, freqs1, abs(cwtmatr1))
    plt.title('信号1 时间频率关系')
    plt.xlabel('时间/s')
    plt.ylabel('频率/Hz')
    
    plt.subplot(223)
    plt.plot(t, y2)
    plt.title('信号2 时间域')
    plt.xlabel('时间/s')
    
    plt.subplot(224)
    plt.contourf(t, freqs2, abs(cwtmatr2))
    plt.title('信号2 时间频率关系')
    plt.xlabel('时间/s')
    plt.ylabel('频率/Hz')
    
    plt.tight_layout()
    plt.show()
    
    
    • 图中不仅可以看到信号中有哪些频率,还可以看到不同的频率成分在什么时间出现。傅里叶变换类似棱镜,可以将不同的信号分解。小波变换类似于显微镜,不仅知道信号中有哪些成分,还可以知道各种成分在什么位置出现。
      在这里插入图片描述
    • 经典的傅里叶变换把信号按正弦、余弦展开,将任意函数表示为具有不同频率的谐波函数的线性叠加,能较好的刻画信号的频率特性。但在时空域上无任何分辨,不能做局部分析。
    • 小波分析优于傅里叶分析之处在于:小波分析在时域和频域同时具有良好的局部化性质,因为小波函数是紧支集,而正弦、余弦的区间是无穷区间,所以小波变换可以对高频成分采用逐渐精细的时域或空域取代步长。,从而可以聚焦到对象的任意细节。

    7.2 简单小波示例

    7.2.1 哈尔小波构建

    • 哈尔基函数可以捕捉到信号的尺度信息,而哈尔小波函数信号的细节信息。哈尔小波变换要做的就是将信号分解到不同的基函数及小波函数上,并求每个函数对应分量的值。
    • 哈尔小波具有如下特点:
      (1)哈尔小波在时域是紧支撑的,即其非零区间为[0,1)
      (2)哈尔小波属于正交小波。
      (3)哈尔小波是对称的。系统的单位冲击响应若具有对称性,则该系统具有线性相位,这对于去除相位失真是非常有利的,哈尔小波是目前唯一一个具有对称性,又是有限支撑的正交小波。
      (4)哈尔小波仅取+1和-1,计算简单。
      (5)哈尔小波是不连续小波,在实际的信号分析与处理中受到了限制。

    7.3 图像多分辨率

    7.3.1 小波多分辨率

    • 多分辨分析是小波分析中最重要的概念之一,将一个函数表示为一个低频成分与不同分辨率下的高频成分,并且多分辨分析能提供一种构造小波的统一框架,提供函数分解与重构的快速算法。
    • 单调性。
    • 伸缩性。
    • 平移不变性。
    • Riesz基。

    7.3.2 图像金字塔

    • 当观察图象时,通常看到是相连接的纹理与灰度级相似的区域,它们相互结合形成物体。如果物体的尺寸较小或者对比度不高,通常采用较高的分辨率观察;如果物体的尺寸很大或者对比度很强,只需要较低的分辨率。
    • 如果物体的尺寸有大有小,或者对比度有强有弱的情况同时发生,那么,以若干个分辨率对他们进行研究将具有优势。
    • 以多分辨率解释图像的一种有效但概念简单的结构是图像金字塔,图像金字塔最初用于机器视觉或者图像压缩,将图像表示为一系列分辨率逐渐降低的集合。
    • 多分辨率分析:将信号或图像用多种分辨率层次表示。在某个分辨率层次上难以检测到的特征,可能很容易在其他分辨率特征上检测出来,每层中包含一个近似图像和一个残差图像,多种分辨率层次联合起来可以称为图像金字塔。
      在这里插入图片描述

    7.3.3 图像子带编码

    • 在子带编码中,一幅图像被分解成一系列频带受限的分量,称为子带。子带可以重组在一起,无失真的重建原始图像。每个子带通过对输入图像进行带通滤波而得到。

    7.4 图像小波变换

    7.4.1 二维小波变换基础

    • 通过小波变换得到N个小波系数,而二维离散小波变换输入的是二维矩阵,每个步骤输出的是近似图像,水平细节,垂直细节和对角细节
    import numpy as np
    import matplotlib.pyplot as plt
    import pywt.data
    
    
    # 中文显示工具函数
    def set_ch():
        from pylab import mpl
        mpl.rcParams['font.sans-serif'] = ['FangSong']
        mpl.rcParams['axes.unicode_minus'] = False
    
    
    set_ch()
    original = pywt.data.camera()
    # Wavelet transform of image, and plot approximation and details
    titles = ['近似图像', '水平细节', '垂直细节', '对角线细节']
    coeffs2 = pywt.dwt2(original, 'haar')
    LL, (LH, HL, HH) = coeffs2
    fig = plt.figure(figsize=(12, 3))
    for i, a in enumerate([LL, LH, HL, HH]):
        ax = fig.add_subplot(1, 4, i + 1)
        ax.imshow(a, interpolation="nearest", cmap=plt.cm.gray)
        ax.set_title(titles[i], fontsize=10)
        ax.set_xticks([])
        ax.set_yticks([])
    fig.tight_layout()
    plt.show()
    
    
    • 基于哈尔小波对图像进行二维小波变换结果
      在这里插入图片描述

    7.4.2 小波变换在图像处理中的应用

    • 小波变换在图像处理中的应用与傅里叶变换类似,基本方法是:
      (1)计算一幅图像的二维小波变换,并得到小波系数
      (2)对小波系数进行修改,保留有效成分,过滤不必要成分
      (3)使用修改后的小波系数进行图像重建

    • 基于小波变换的图像去噪步骤
      (1)图像小波变换。选择一个小波,计算噪声图像的小波系数。
      (2)对细节系数通过阈值进行过滤。选择一个细节系数阈值,并对所有细节系数进行阈值化操作。
      (3)基于阈值化过滤后的细节系数及原始近似系数,使用小波变换对图像进行重建。

    import numpy as np
    import matplotlib.pyplot as plt
    import pywt
    from skimage.restoration import denoise_wavelet
    
    
    # 中文显示工具函数
    def set_ch():
        from pylab import mpl
        mpl.rcParams['font.sans-serif'] = ['FangSong']
        mpl.rcParams['axes.unicode_minus'] = False
    
    
    set_ch()
    x = pywt.data.ecg().astype(float) / 256
    
    sigma = .05
    x_noisy = x + sigma * np.random.randn(x.size)
    x_denoise = denoise_wavelet(x_noisy, sigma=sigma, wavelet='sym4', multichannel=False)
    
    plt.subplot(311)
    plt.title('原始信号')
    plt.plot(x)
    
    plt.subplot(312)
    plt.title('加噪图像')
    plt.plot(x_noisy)
    
    plt.subplot(313)
    plt.title('去噪图像')
    plt.plot(x_denoise)
    
    plt.show()
    
    • 小波对一维信号去噪
      在这里插入图片描述
    import matplotlib.pyplot as plt
    from skimage.restoration import (denoise_wavelet, estimate_sigma)
    from skimage import data, img_as_float
    from skimage.util import random_noise
    
    
    # 中文显示工具函数
    def set_ch():
        from pylab import mpl
        mpl.rcParams['font.sans-serif'] = ['FangSong']
        mpl.rcParams['axes.unicode_minus'] = False
    
    
    set_ch()
    original = img_as_float(data.coffee())
    
    plt.subplot(221)
    plt.axis('off')
    plt.title('原始图像')
    plt.imshow(original)
    
    sigma = 0.2
    noisy = random_noise(original, var=sigma ** 2)
    
    plt.subplot(222)
    plt.axis('off')
    plt.title('加噪图像')
    plt.imshow(noisy)
    
    im_haar = denoise_wavelet(noisy, wavelet='db2', multichannel=True, convert2ycbcr=True)
    plt.subplot(223)
    plt.axis('off')
    plt.title('使用haar去噪后')
    plt.imshow(im_haar)
    
    # 不同颜色通道的噪声平均标准差
    sigma_est = estimate_sigma(noisy, multichannel=True, average_sigmas=True)
    im_haar_sigma = denoise_wavelet(noisy, wavelet='db2', multichannel=True, convert2ycbcr=True, sigma=sigma_est)
    plt.subplot(224)
    plt.axis('off')
    plt.title('使用haar with sigma去噪后')
    plt.imshow(im_haar_sigma)
    plt.show()
    
    
    • 二维小波图像去噪
      在这里插入图片描述

    7.5 小结

    • 本章主要介绍小波相关内容。首先介绍了小波变换和傅里叶变换的区别,并从感性认识介绍了小波的优点,其次介绍了小波信号处理过程,包括小波构建、小波变换、小波逆变换,然后介绍了图像多分辨率,最后介绍了小波变换在图像处理中的应用。

    参考资料

    1. 岳亚伟《数字图像处理与Python实现》人民邮电出版社
    展开全文
  • Radon变换历史:Radon变换是J. Radon于1917年提出的。在Fourier变换及它们对应的卷积可以快速计算之前,Radon变换的计算几乎没有引起人们的兴趣。现在Radon变换已经成为医学成像和许多遥感成像等的主要工具而受到...
  • 该程序是小波变换、s变换、傅里叶变换对雷克子波的时频分析
  • 多分辨率分析小波变换--简介

    万次阅读 2013-03-03 09:50:50
    傅里叶变换的基础函数的正弦函数,而小波变换基于一些小型波,称为“小波”,它具有变化的频率和有限的持续时间。这就允许它们对图像提供一张等效的乐谱,不光阐明了要演奏的音符(频率),而且阐明了要何时演奏。...
  • Matlab讲义之胡广书现代信号处理教程-第十章 离散小波变换多分辨率分析.rar 喜欢的下
  • 序言 什么是小波 “小波”(wavelet)就是一种“尺度”很小的波动,并具有时间和频率...■小波变换基于一些小型波,称为小波,具有变化的频率和有限的持续时间。 ◆傅里叶变换反映的是图像的整体特征,其...
  • 在这项工作中,我们研究了使用小波变换的概念验证多分辨率分析,这是一种在信号处理和表示中使用的流行数学和分析框架,并且我们研究了其在无线传感器网络中压缩图像数据领域的应用。 。 所提出的方法包括小波变换...
  •   虽然自20世纪50 年代末起傅里叶变换就一直是基于变换的图像处理的基石,但近年来一种新的称为小波变换的变换使得压缩、传输和...1987年,首次证明小波是―种全新且有效的信号处理和分析方法,称为多分辨率理论的基
  • 小波变换(1)多分辨率分析:尺度函数和小波函数

    万次阅读 多人点赞 2014-04-23 00:40:32
    在小波分析中,常出现空间V和W,以及Scaling
  • 该方法充分利用了小波变换多分辨率分解的特性,结合分数傅里叶变换,不仅有效完成了图加密,而且实现了各图像独立的处理,加密后,每个图像都会有自己独立的密钥。分析小波变换类型和变换级次对加、解密效果的...
  • matlab小波变换的代码

    2020-04-11 10:43:27
    %% 3.MALLET分解算法(圆周卷积的快速傅里叶变换实现) sig1=ifft(fft(y).*fft(h)); % 低通(低频分量) sig2=ifft(fft(y).*fft(g)); % 高通(高频分量) figure(5); % 信号图 subplot(2,1,1) plot(real(sig1)); ...
  • 本文基于小波分析提出由高分辨率图像可得到各种较低分辨率图像的分辨方法符合人的视觉系统特性。为了提高图像的空间分辨率,本文根据小波子带间的相关性等特点,提出由低分辨率的高频细节预测高分辨率的高频细节方法...
  • 第四章多分辨率分析与正交小波变换 .ppt
  • 第四章多分辨率分析与正交小波变换.ppt
  • 第三章 多分辨率分析与离散序列的小波变换 3.1 概述 3.2 分辨率信号分解与重建的基本概念 3.3 尺度函数和小波函数的一些重要性质 3.4 由多分辨率分析引出采样率滤波器组 3.5 Mallat 算法实现中的一些问题 ...
  • 1.如何理解时间分辨率和频率分辨率? In this section we will take a closer look at the resolution properties of the wavelet transform. Remember that the resolution problem was the main reason why ...
  • 浅谈傅里叶变换、小波变换、HHT变换

    千次阅读 多人点赞 2019-08-29 12:18:42
    浅谈傅里叶变换、小波变换、HHT变换一、傅里叶变换1.1傅里叶变换介绍二、小波变换2.1小波变换正反变换公式2.2小波变换适应场景及其优缺点2.3小波变换的应用三、HHT变换3.1HHT产生的背景3.1 HHT变换介绍3.2 HHT对信号...
  • 第四章-多分辨率分析与正交小波变换.ppt
  • 利用小波变换多分辨率特性分析降水时间序列,其小波分解级数应是确定的。为解决此问题,提出一种新的框架:首先利用à trous小波变换对降水时间序列进行分解,然后对分解的序列进行尺度熵(MSE)分析以获取原始...
  • 前者用多分辨率分析的多孔算法计算变形条纹解析信号的离散小波变换, 寻找每一位置的小波系数在尺度方向上的模极大值点, 通过提取该点的相位就可以得到原条纹的相位值; 后者是在前者的基础上, 通过插值改变条纹信号的...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 3,247
精华内容 1,298
关键字:

小波变换多分辨率分析