Graph Signal Processing – Part II: Processing and Analyzing Signals on Graphs(文献翻译)
目录
摘要
1.引言
2. 问题陈述:一个说明性的例子
3.图上的信号和系统
3.1 邻接矩阵与图信号移位
3.3 图离散傅里叶变换(GDFT)
3.4 在GDFT域图上的系统
3.5 邻接矩阵频谱域的图信号滤波
3.5.1 邻接矩阵的归一化
3.5.2 邻接矩阵特征向量的谱排序
3.5.3 谱域滤波器设计
3.5.4 图传递函数上系统的多项式(切比雪夫)近似
3.5.5 图上的逆系统
3.6 基于拉普拉斯函数的图傅里叶变换
3.7 拉普拉斯谱域的排列组合和滤波
3.8 使用图拉普拉斯定义图系统
摘要
图上的数据分析是对在非规则但结构化的图域上获取的数据进行信息处理。本专题的第一部分的重点是基本的和高阶的图属性、图拓扑和图的谱表示。第一部分还为顶点聚类和图分割建立了严格的框架，并说明了图在各种数据关联任务中的强大功能。第二部分从这些概念着手，以解决围绕图上的数据/信号处理的算法和实际问题，也就是说，重点是对图上的确定性和随机数据的分析和估计。通过对多传感器温度场估计的一个简单直观、说明性强、通用性强的案例研究，介绍了与图信号相关的基本思想。利用图信号移位算子定义了图上系统的概念，推广了传统学习系统的相应原理。图信号和系统的谱域表示的核心是图离散傅里叶变换(GDFT)，它是基于邻接矩阵和图拉普拉斯矩阵的特征分解而定义的。然后用谱域表示作为基础来介绍图信号滤波的概念和解决它们的设计，包括切比雪夫多项式近似级数。本文提出了与图信号采样相关的思想，特别是通过图子采样进行数据降维这一具有挑战性的课题，并进一步将其与压缩感知联系起来。接下来复习图上时变信号的原理和与随机图信号相关的基本定义。顶点-频谱联合域的局部图信号分析称为顶点-频率分析，它可以看作是经典时频分析在信号图域的扩展。介绍了与局部图傅里叶变换(LGFT)相关的重要课题，以及它的各种形式，包括图谱和顶点域窗以及反演条件和反演关系。建立了带谱变窗的LGFT与谱图小波变换(SGWT)之间的联系。进一步讨论了利用多项式(Chebyshev)逼近实现LGFT和SGWT的方法，并给出了实例支持。最后，介绍了顶点-频率表示的能量形式，以及它们与经典时频分析的关系，包括满足边缘性质的顶点-频率分布。这篇文章有许多的例子来支持。
1.引言
图是一种不规则的结构，它自然地解释了数据的完整性，然而，传统的方法已经建立在机器学习和信号处理之外，主要集中于分析底层的图，而不是处理图上的信号。另一方面鉴于多传感器、多节点的迅速增加,可能记录在不规则或特定的网格,受益于图表可以说明空间感知意识,物理直觉和传感器重要性的能力,结合固有的局部和全局传感器联合,这将非常有利于分析等结构化数据信号图。因此我们的第二部分目的是建立一个图形信号之间的共同语言是观察不规则信号域,和一些最基本的学习范式——学习系统、信号处理和数据分析，如光谱分析、系统传递函数，数字滤波器的设计,参数估计和最优化去噪。在经典的数据分析和信号处理中，信号域由等距的瞬时时间或均匀网格上的一组空间传感点决定。然而,越来越多的实际数据传感领域甚至可能不是相关的物理维度的时间和/或空间,通常会表现出各种形式的不规则性,,例如,在社会或与web相关的网络,传感点及其连接关系属于特定对象/节点和特殊拓扑结构的链接。应该注意的是,即使数据获得良好的时空域定义,通过图引入信号样本之间的新关系,可产生对分析新的见解和提供增强的数据处理(例如,基于邻域局部相似性)。因此，图相对于经典数据域的优势是，图在问题定义中自然而全面地考虑了不规则数据关系，以及分析中相应的数据连接。
为了在图上的信号/数据的基本思想上建立直观感受，第2节首先考虑一个简单而通用的多传感器温度估计的例子。第3节介绍了图上信号和系统的基本概念，包括基本定义、运算和变换，推广了传统信号处理的基础。本文从全面介绍和介绍一种新颖的等度量图信号移位算子开始，对图上的系统进行了解释。然后，基于邻接矩阵和图拉普拉斯函数定义了图傅立叶变换，并以此为基础引入了图信号滤波的概念。在第4节中，我们将回顾与图信号采样有关的各种思想，特别是它们的子采样这一具有挑战性的主题。第5节和第6节介绍了图上时变信号的概念，并介绍了与随机图信号相关的基本定义。第7节讨论的局部图信号可以同时在顶点-频域进行表征。本节还涵盖了局部图傅立叶变换的重要主题，其变换的各种形式，与框架的关系和与图小波变换的联系。顶点频率表示的能量也被考虑，以及它们与经典时频分析的关系。
2. 问题陈述:一个说明性的例子
考虑使用多传感器来测量感兴趣区域内的温度场。根据特定地理区域对当地用户的重要程度来选择温度感知位置，共N = 16个感知点，如图1(a)所示。温度场用{x(n)}表示，n为传感器指标，其值的快照见图1(b)。每个测量的传感器信号可以用数学表示为
其中s(n)为理想测量条件下的真实温度，ε(n)包括局部环境对传感器读数的不利影响或传感器故障(被称为“噪音)。为了说明，在我们的研究中ε(n)被建模为白色、零均值、高斯分布，标准差σε = 2，即:
。将其添加到信号s(n)中，得到x(n)中的信噪比为 14.2 dB。
备注1:经典数据分析需要将图1(a)中典型的不规则空间温度传感排列重新排列成图1(b)中所示的线性结构。显然，这种“字典式”不适合利用与实际传感器位置相关、由地形决定的信息。这使得多传感器温度场的经典分析不适用(或是次优)，因为性能严重依赖于选择的传感器排序方案。这就证明了即使是最常规的多传感器测量设置需要一个比对应于经典信号处理框架的标准线性估计结构更复杂的估计结构。
图1
为了引入一种针对图1所示温度场的情境感知降噪方案，我们从一个局部信号平均算子开始，探索这个问题的图论框架。在经典的分析中，这可以通过移动平均算子来实现，例如，通过对邻近数据样本进行平均，或图1(b)中线性数据设置等效邻近传感器，以及对每个传感点进行平均。从物理上来说，这种局部邻域应该包括邻近的传感点，但只能是那些具有相似的气象属性的传感点，这些气象属性由传感器距离、高差和其他地形特定属性所定义。
换句话说，由于图1中的传感器网络测量了一组间隔不规则的传感器的相关温度，所以有效的估计策略应该包括用标准方法(线性路径图)无法实现的领域知识。
为了说明基于局部信息(基于邻域)的方法的优点，考虑遥感点n = 3(低地)，n = 6的(山脉)和n = 8(海岸)，如图2(a)所示，然后给出了各测点的累积温度(m在n附近)
使一个感测点n的局部平均温度可由累计温度y(n)除以所包含感测点的个数(局部邻域大小)得到。例如，对于图2(a)所示的感知点n = 3和n = 6，“域信息感知”的局部估计采用如下形式
为方便起见，传感点之间的全连接关系现在可以排列成矩阵形式：
y = x + Ax (4)
式(5)中邻接矩阵A表示感知位置的连通性结构;每个y(n)的计算都应该涉及到这个局部连通性结构。
这个简单的真实世界的例子可以在图形信号处理框架中解释:
- 信号被测量的传感点被指定为图节点(vertices)，如图1所示。
- 如图2(a)所示，表示感知点之间物理上的连接的顶点到顶点的线看作图边(edges)，
- 顶点和边构成一个图，如图2(b)，一个新的结构丰富的信号域，
- 这个图形，不是一个标准的传感点矢量，被用来分析和处理数据，因为它展示了空间和物理领域的感知能力，
- 测得的温度现在被解释为图上的信号样本，如图3所示。
- 与传统的信号处理类似，这种新的图信号可能在同一图上有多种形式，并且可能包含噪声，
- 因此，通过关系(4)，我们引入了一个简单的图形系统，用于物理和空间感知信号平均值(图形上的线性一阶系统)。
为了强调我们对特定传感器的依赖程度(即对传感器相关性建模)，可以提出一个权重方案，形式如下
其中wnm是加权矩阵W的元素。图边及其权值Wnm的定义有三类方法:
•已经有了物理意义上良好定义的边和权重，
•基于顶点位置的几何学定义边和权值，
•基于数据相似度的方法来学习底层的图拓扑。
定义边权值的三种方法在本专题的part3中都有详细介绍。
由于在我们的地理温度测量案例中，图权值不属于明显的和物理上定义良好的边和权值的类别，我们将使用基于“顶点位置几何学”的方法来定义边和权值。以这种方式，根据顶点水平距离rmn和高度差hmn，计算相邻顶点的权值元素Wnm：
其中α和β是合适的常数。由此得到的权矩阵W如(6)所示。
图2
图3 图形信号值用两种方式表示:顶图用长度与信号值成比例的垂直线表示，底图用热度图指定顶点上的信号值。
基于(4)，得到了累积温度的加权图信号估计：
为了产生无偏估计，代替(4)和(8)中的累积和，每个y(n)的估计内的加权系数应该总和为一。这可以通过归一化(10)实现，由
其中对角归一化矩阵D的元素等于度矩阵的元素，是一个随机漫步(扩散)移位算子
既然我们已经定义了图的顶点和边的权值，我们可以使用第一部分第4.3节中给出的不可知数据聚类方法，根据图的拓扑结构来聚类这个图的顶点。图4为基于图拉普拉斯矩阵的三个最光滑(λ最小对应的列特征向量)特征向量u1、u2、u3(不含常数特征向量u0)得到的聚类结果
图4 观察该图的正确聚类，海边地区(蓝色)，低山(红色)，低地(黄色)，和高山(绿色)。
图拉普拉斯矩阵L = D−W，其值如(7)所示。
求取L的特征值和特征向量：
L=[3.77 0 0 -0.97 -0.91 0 0 0 0 0 0 0 0 -0.05 -0.90 -0.94;0 1.19 -0.03 0 0 0 0 0 0 -0.37 0 0 0 -0.78 0 0;0 -0.03 2.93 0 0 -0.96 0 -0.95 -0.98 0 0 0 0 0 0 0;-0.97 0 0 2.81 0 0 0 0 0 0 0 0 0 0 -0.88 -0.96;-0.91 0 0 0 1.94 0 0 0 0 0 -0.06 -0.01 -0.01 0 -0.94 0;0 0 -0.96 0 0 1.94 0 0 -0.97 -0.01 0 0 0 0 0 0;0 0 0 0 0 0 1.02 0 0 -0.62 -0.40 0 0 0 0 0;0 0 -0.95 0 0 0 0 1.89 -0.94 0 0 0 0 0 0 0;0 0 -0.98 0 0 -0.97 0 -0.94 2.89 0 0 0 0 0 0 0;0 -0.37 0 0 0 -0.01 -0.62 0 0 1.24 -0.25 0 0 0 0 0;0 0 0 0 -0.06 0 -0.40 0 0 -0.25 1.56 0 0 -0.85 0 0;0 0 0 0 -0.01 0 0 0 0 0 0 0.93 -0.92 0 0 0;0 0 0 0 -0.01 0 0 0 0 0 0 -0.92 0.94 0 -0.01 0;-0.05 -0.78 0 0 0 0 0 0 0 0 -0.85 0 0 1.68 0 0;-0.90 0 0 -0.88 -0.94 0 0 0 0 0 0 0 -0.01 0 2.74 0;-0.94 0 0 -0.96 0 0 0 0 0 0 0 0 0 0 0 1.91;][U,lamda]=eig(L) %求矩阵L的特征值和特征向量,其中lamda的对角线元素是特征值,其中U的X的列是相应的特征向量,--------------------结果----------------------U =0.2299 0.1076 -0.1769 0.3208 -0.0036 -0.0069 0.0007 0.0013 0.0032 -0.0005 0.0171 0.0054 0.0001 -0.0002 -0.0626 -0.89280.2555 -0.0313 -0.2073 -0.2865 0.5029 -0.5442 -0.0511 -0.0011 0.3221 -0.0454 -0.3836 -0.0625 -0.0047 0.0105 0.0004 -0.00420.2527 -0.3456 0.2418 0.0682 -0.0002 -0.0003 -0.0001 -0.0000 0.0032 -0.0061 -0.0072 -0.0013 0.3261 -0.8044 -0.0001 0.00010.2295 0.1090 -0.1777 0.3325 -0.0106 -0.0424 0.2667 -0.0014 -0.0038 0.0005 0.0737 -0.4801 -0.0001 0.0002 -0.6410 0.27090.2300 0.1087 -0.1714 0.3056 -0.0059 0.0724 -0.6459 0.0070 0.0159 -0.0020 -0.0880 0.5096 -0.0000 0.0000 -0.2593 0.23470.2534 -0.3478 0.2442 0.0696 -0.0070 0.0016 -0.0003 0.0003 -0.1067 -0.6915 0.0031 0.0006 -0.4985 0.1110 0.0000 -0.00000.2596 -0.0246 -0.2249 -0.3144 -0.6273 0.1044 0.0354 -0.0027 0.6031 -0.0894 0.0705 0.0159 -0.0005 0.0002 -0.0008 0.00090.2534 -0.3494 0.2466 0.0717 -0.0041 0.0089 0.0013 -0.0002 0.1039 0.7067 0.0047 0.0008 -0.4739 0.1128 0.0000 -0.00000.2534 -0.3490 0.2460 0.0712 -0.0050 0.0060 0.0006 0.0000 -0.0026 -0.0127 0.0029 0.0005 0.6485 0.5725 0.0000 -0.00000.2604 -0.0283 -0.2204 -0.3115 -0.3533 -0.4513 -0.0648 0.0030 -0.6469 0.1012 0.1538 0.0256 0.0028 -0.0021 -0.0005 0.00100.2569 -0.0180 -0.2200 -0.2767 0.0733 0.6300 0.0574 0.0010 -0.3005 0.0432 -0.5460 -0.1129 -0.0010 0.0018 0.0068 -0.00970.2724 0.4772 0.4314 -0.1211 0.0002 -0.0008 0.0097 0.7051 0.0028 -0.0002 0.0007 -0.0034 0.0000 -0.0000 0.0013 -0.00040.2722 0.4752 0.4281 -0.1186 0.0001 -0.0005 0.0009 -0.7090 -0.0031 0.0002 -0.0004 0.0016 0.0000 -0.0000 -0.0016 -0.00100.2558 -0.0206 -0.2155 -0.2722 0.4722 0.2813 0.0516 -0.0005 0.0200 -0.0056 0.7022 0.1221 0.0021 -0.0044 -0.0012 0.01860.2293 0.1098 -0.1737 0.3229 -0.0090 0.0168 -0.3016 -0.0006 0.0170 -0.0023 0.0956 -0.4986 -0.0000 0.0000 0.6486 0.17570.2288 0.1084 -0.1783 0.3344 -0.0114 -0.0757 0.6377 -0.0021 -0.0225 0.0031 -0.1015 0.4782 -0.0000 0.0000 0.3118 0.2076lamda =0.0023 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0.0116 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0.0210 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0.0539 0 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0.7176 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 1.2864 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 1.5075 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 1.8550 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 1.8844 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 1.9151 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 2.7657 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 2.8632 0 0 0 00 0 0 0 0 0 0 0 0 0 0 0 3.8299 0 0 00 0 0 0 0 0 0 0 0 0 0 0 0 3.8936 0 00 0 0 0 0 0 0 0 0 0 0 0 0 0 4.0725 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4.7003
即使是这样一个简单的图聚类方案也能够识别不同的有物理意义的地理区域。这也意味着温度估计可以粗略地在每个聚类内执行，甚至可以作为一个独立的图(见第一部分的图切割，第4节)，而不是在整个传感器网络。
上面介绍的图数据估计框架是非常通用的，可以应用到许多不同的场景中，在确定了一个合适的图拓扑结构后，我们希望对这些图上获得的数据进行估计，这也是本专题的主题。
3.图上的信号和系统
在经典数据分析中，信号采样是在连续的、等间距瞬时进行的。这就决定了信号样本的顺序，x(n)前面是x(n−1)，后面是x(n+1)。因此，数据样本之间的“时间距离”是标准数据处理算法的一个固有参数。采样时刻之间的关系也可以用图形形式表示，信号被采样时对应的顶点和相应的边定义了线性采样(顶点)顺序。在经典应用中，瞬时采样的等距性质可以用所有边的权值来表示(例如，归一化为1)，如图5所示。
图5:定义在等距离散时间网格上的经典时域信号的有向路径图表示。
在离散时间中定义的算法(例如基于DFT或其他类似数据变换的算法)，通常假定所分析的信号具有周期性，即采样x(N−1)与采样x(0)连续，为一个循环序列。注意，这种情况对应于图6所示的循环图，它允许我们在许多标准数据转换使用该模型，例如DFT、DCT、小波，并基于这些转换定义对应图的其他处理算法。
图6 有向循环图
一般无向图(包括循环形)上的信号是通过将实(或复数)数据值x(n)与每个顶点相关联来定义的，如图7和图8所示。这样的信号值可以以矢量形式排列：
因此，一个图可以看作是一个广义的信号域。
图7
图8
一般来说，对于在顶点n观察到的图信号的任何线性处理方案，可以被定义为这个顶点的信号值x(n)与邻近顶点的信号样本x(m)的线性组合，即
其中，Vn为节点n附近的顶点集，h(m,n)为尺度系数。
注2:(12)中的估计形式高度依赖于顶点;它仅在正则图这样一个非常特定的情况下是顶点不变的，其中Vn是顶点n的k邻域，h(n,m) = h(n−m)。
现在我们继续使用图形上的移位来定义各种形式的顶点不变(vertex-invariant)滤波函数。然后，这些将被用来介绍有效的图信号处理方法
3.1 邻接矩阵与图信号移位
考虑一个图信号x，它的`x(n)是顶点n的观测样本。图上的信号移动可以定义为信号样本x(n)从原顶点n沿长度为1的所有路径移动，即K = 1，从顶点n开始。如果这样移动的信号用x1表示，则可以用图邻接矩阵A来定义其值
例1:作为图信号及其移位效果的说明，考虑图6(a)中圆形图上的信号。
原始信号x如图9(a)所示，其移位版次x1如图9(b)所示。图8 (a)中无向图上的另一个简单信号如图10(a)所示，其移位版次x1 = Ax如图10(b)所示
信号经过两次图移位后，再经过一次移位x1 = Ax得到。结果，两次移位，图形信号为：
因此，一般情况下，图上的m次平移信号为：
注3:与标准移位算子一样，图信号的二阶移位是通过对已经移位过的信号进行移位而得到的。移位算子的作用由邻接矩阵A来假设。
3.2 基于图移信号的系统
与基于标准线性移位的系统非常相似，图上的系统可以表示为图信号x及其图移位效果, m = 1,2，…, M−1的线性组合。系统在图上的输出信号可以写成：
其中，根据定义，h0, h1，…， hM−1为系统系数。对于循环(经典线性系统)图，这种关系可简化为著名的的有限关系脉冲响应(FIR)滤波器，即
矩阵Am在图中描述了长度为K = m的游走遍历(参见第1部分m2的性质)，输出图信号y(n)为输入图信号值和所考虑顶点n的(M−1)邻域的顶点上观察到的信号值的线性组合。
注4：当最小多项式和特征多项式的次数相同时，一个有物理意义的系统序(M−1)应小于顶点数N，即M≤N，在经典信号分析中对应的条件为:
(15)中的系统脉冲响应系数hM的数量M应该小于或等于信号采样总数N(对于图9中的循环图，表示均值图信号移位为m = 0,1,2，…，N−1，因为m = N的移位减少为m = 0的移位，m = N + 1的移位等价于m = 1的移位，以此类推)因此，一般情况下，系统序(M−1)应小于邻接矩阵a的最小多项式的次数Nm，更多细节见第一部分，3.1节
注5:任何M - 1阶≥Nm的系统都可以简化为Nm - 1阶的系统。
注6:如果系统阶大于或等于最小多项式的次，即M−1≥Nm，则存在多个系统对给定的输入信号产生相同的输出信号。图上所有这样的系统称为等价的。
最后三个注释中的说明及证明将在第3.5.3节中更详细地说明
例2:图11(a)为给出的信号，在图上有一个线性系统，由系数h0 = 1, h1 = 0.5定义。这个系统在图上对应于一个经典的一阶加权移动系统。然后输出的图信号表示顶点n处的信号值和K = 1附近的信号值的加权平均值。输出图形信号如图11(b)所示， y = x + 0.5 Ax.
图形上的一般系统。图上的一个系统在顶点域中可以定义为：
其中H(A)为顶点域系统(filter)函数。如果图上的系统满足以下性质，则该系统是线性的且移位不变的:
1. 线性：
2. 移不变性：
输入序列的移位仅引起输出序列做相同移位的系统为移不变系统
注7:图上的系统由是线性的且移位不变的，因为
3.3 图离散傅里叶变换(GDFT)
经典的探索性数据分析通常采用频谱域(傅里叶)信号估计;这就产生了许多简单而有效的算法。虽然标准的频谱分析在时间和频率上都采用等距网格，遵循图上系统的思想，我们接下来表明图信号的谱域表示是基于邻接矩阵或图拉普拉斯矩阵的谱分解。
信号x的图形傅里叶变换被定义为：
其中X为GDFT系数的向量，U是一个矩阵，其列表示邻接矩阵a的特征向量。用X(k)表示向量X的元素，k = 0,1，…，N−1，由前面可知，对于无向图，邻接矩阵是对称的，即A^T = A，对称矩阵的特征矩阵满足这个性质:
因此，图傅里叶变换向量X的元素X(k)表示的是考虑的图信号x(n)在A(一个基函数)的第k个特征向量上的投影：
这样，图形离散傅里叶变换作为标准正交基函数，可以被解释为一组投影(信号分解)到u0,u1，…，uN−1特征向量集合上。
然后，(18)直接得到反离散傅里叶变换图：
即
例如，从图6一个循环图，(19)和(21)中的离散傅里叶变换将图形简化为标准离散傅里叶变换(DFT)对。因此，(19)中的变换及其在(21)中的逆称变换为图离散
傅里叶变换(GDFT)与图反离散傅里叶变换(IGDFT)。
3.4 在GDFT域图上的系统
考虑(17)中所定义的一个图上的一般系统，
利用邻接矩阵的谱表示，我们得到：
用系统上的一个图传递函数：λ是A的特征值矩阵。
先把这个关系式(23)与U^−1相乘，得到:
由(18)可知，U^−1*y和U^−1*x分别为输出图信号y和输入图信号x的GDFT，使图关系上的谱域系统可写为：
在顶点域的输出图形信号可以计算为：
(26)中图系统的元素形式为：
其中λk表示邻接矩阵A的第k个特征值。由(24)和上式，我们现在可以在图上定义系统的传递函数，形式为:
注8:(15)中的经典线性系统可以直接从(28)中对应的循环图中得到。这是因为有向循环图的邻接矩阵的特征值(有关有向循环图的详细信息，请参阅第一部分3.2.1节)与经典DFT中单位圆上采样相同。
与经典信号处理中的z变换类似，对于图上的系统，我们也可以在z域中引入系统传递函数。
系统在图上的z域传递函数定义为：
n = 0,1，…, M−1。显然，从(28)，可以得知：
然而，满足关系的任意图信号(x(n)和y(n))的z-变换的定义并不简单，这限制了z-变换在图上的应用。这将在3.10节中进行更详细的讨论。
3.5 邻接矩阵频谱域的图信号滤波
图移信号的能量为
向量范数
但是，如图10所示，一般情况下，移位信号的能量与原始信号的能量是不相等的
另一方面，在图信号处理中，通常希望图移不增加信号能量。下面介绍了一个这样的图移位算子。
注9:利用矩阵二范数，直接证明了图移信号Ax与原图信号x的能量之比满足这个关系
其中，
3.5.1 邻接矩阵的归一化
我们可以使用归一化邻接矩阵，定义为
作为一个图上任何系统中的一个图移位算子。虽然这种归一化仍然不能在图上进行等距移位，但这样移位的信号的能量保证不大于原图信号的能量，因为
等号成立只适用于与λ max对应的特征向量成比例的特定信号
图上的基本移位、图上的系统和图的谱域表示都可以用(31)中的归一化邻接矩阵来实现，方法与原始邻接矩阵相同。一个不适用于标准邻接矩阵的重要性质是邻接矩阵的归一化产生一个更简单的特征向量和特征值排序方案，如下文。
3.5.2 邻接矩阵特征向量的谱排序
为了在图上进行有物理意义的低通和高通滤波，我们需要建立谱阶的概念。这反过来需要一个准则来将特征向量(对应于GDFT基函数)分类为慢变化和快变化的特征向量。
注10:在经典的傅里叶分析中，基函数是按其频率排序的，例如低通(慢变)基函数是低频率的谐波函数。另一方面，作为信号分解基础的图邻接矩阵特征向量的频率的概念没有定义，我们必须找到另一个准则对特征向量进行分类或排序。再一次，我们从经典的傅立叶分析中得到启示，即“信号变化”的能量。
信号能量的变化。对于图信号，第一个图差值可以定义为原始图信号与其图差的差值，即
与经典分析类似，信号变化的能量可以定义为一个图形信号x的一阶差分的能量，并采取如下形式:
当图信号具有邻接矩阵A的特征向量x = u的特定形式时，该特征向量变化的能量等于
其中用归一化邻接矩阵Anorm作为范数来限定移位的图信号的能量。在推导中我们也使用了
现在，E∆u的值越小，说明u是缓慢变化的，E∆u = 0表明信号是恒定的，而E∆u的值越大，说明u在时间上的变化越快。(32)中的形式也被称为基函数/特征向量的二范数总变差。
因此，如果基函数u的变化具有很大的能量，那么特征向量u就可以被认为对图信号的高谱含量来说是长的。
在数学中，基函数是函数空间中特定基底的元素。函数空间中的每个连续函数可以表示为基函数的线性组合，就像向量空间中的每个向量可以表示为基向量的线性组合一样。
备注11:从(32)开始，当λ = λ max时，图信号的变化率的能量是最小的，并且随着λ的减小而增大(见第1部分图10)。
我们已经建立了基于相应的特征值的特征向量一个判据，我们将继续在图上定义一个理想低通滤波器。低通滤波在图像域的作用是这样一个过滤器应该通过变化率比截止特征值λc(参见截止频率)所定义的变化率要慢的信号分量(A的特征向量)。而所有表现出比截止特征值λ c的变化更快的信号分量(特征向量)应该被抑制。因此，图域中的理想低通滤波器定义为：
例3在图12(a)中，我们得到一个图形信号，为这个图邻接矩阵的两个特征向量的线性组合为 x=(图的邻接矩阵的特征向量在第一部分图10中给出),信号受到信噪比为2.7dB的高斯白噪声 ε的干扰，噪声图信号如图12(b)所示，接下来用截止特征值为λc = 1的图滤波器对这个有噪声的信号进行滤波。输出信号xf如图12(c)所示，在信噪比out = 18.8dB的情况下，这种滤波方法可以使信号质量提高16.1dB。
备注12:当时，特征向量的变化率能量与经典的基于DFT的滤波一致。
------------------------------------------------------------------------------------
代码如下(上matlab代码：用于生成针状图；下python ：用于SVD滤波) 因高斯噪声为随机噪声，所以下面代码所生成的图与上方图12有差别A=[0,1,1,1,0,0,0,0;1,0,1,0,1,0,0,0;1,1,0,1,1,0,0,0;1,0,1,0,0,0,1,0;0,1,1,0,0,1,0,1;0,0,0,0,1,0,0,1;0,0,0,1,0,0,0,1;0,0,0,0,1,1,1,0]; [U,lamda]=eig(A); %求矩阵A的特征值和特征向量,其中lamda的对角线元素是特征值,其中U的X的列是相应的特征向量,u6=U(:,7); u7=U(:,8); x=2*u6+3.2*u7%绘制系统输出x的谱图 figure y=0:1:7; stem(y,x)----------------结果------------------ x =0.47040.98141.04470.52802.06781.68340.90801.9368import numpy as np import random import matplotlib.pyplot as pltdef gen_gaussian_noise(signal,SNR):""":param signal: 原始信号:param SNR: 添加噪声的信噪比:return: 生成的噪声"""noise=np.random.randn(*signal.shape) # *signal.shape 获取样本序列的尺寸noise=noise-np.mean(noise)signal_power=(1/signal.shape[0])*np.sum(np.power(signal,2))noise_variance=signal_power/np.power(10,(SNR/10))noise=(np.sqrt(noise_variance)/np.std(noise))*noisereturn noise## 1.待处理信号(9个采样点) t = np.arange(0,9,1) t1 = np.arange(0,9,1)u6=[-0.3811,-0.1634,-0.2586,-0.2672,0.3441,0.4944,0.1563,0.5500,0] u7=[0.3852, 0.4089,0.4881, 0.3320, 0.4311,0.2171,0.1860,0.2615,0] x=2*np.array(u6)+3.2*np.array(u7) x=np.array(x) r= gen_gaussian_noise(x,2.7) xr = x+r## 2.一维数组转换为二维矩阵 x2list = [] for i in range(3):x2list.append(xr[i*3:i*3+3]) x2array = np.array(x2list)## 3.奇异值分解 U,S,V = np.linalg.svd(x2array) S_list = list(S) ## 奇异值求和 S_sum = sum(S) ##奇异值序列归一化 S_normalization_list = [x/S_sum for x in S_list]## 4.画图 X = [] for i in range(len(S_normalization_list)):X.append(i+1)fig1 = plt.figure().add_subplot(111) fig1.plot(X,S_normalization_list) fig1.set_xticks(X) fig1.set_xlabel('Rank',size = 15) fig1.set_ylabel('Normalize singular values',size = 15) plt.show()## 5.数据重构 K = 1 ## 保留的奇异值阶数 for i in range(len(S_list)):# S_list[i+K] = 0.0if S_list[i]<1:S_list[i]=0S_new = np.mat(np.diag(S_list)) reduceNoiseMat = np.array(U * S_new * V) reduceNoiseList = [] for i in range(len(x2array)):for j in range(len(x2array)):reduceNoiseList.append(reduceNoiseMat[i][j])## 6.去燥效果展示 fig2 = plt.figure().add_subplot(111) fig2.plot(t,list(x),'y',label = 'Original data') fig2.plot(t,list(xr),'b',label = 'noisy data') fig2.plot(t1,reduceNoiseList,'r-',label = 'Processed data') fig2.legend() fig2.set_title('Rank is 1') fig2.set_xlabel('Sampling point',size = 15) fig2.set_ylabel('Value of data',size = 15) plt.show()
------------------------------------------------------------------------------------
3.5.3 谱域滤波器设计
我们将用表示定义在图上的系统的期望图传递函数。然后，具有该传递函数的系统既可以在谱域实现，也可以在顶点域实现。
在频谱领域，实现是直接的，可以通过以下三个步骤来实现:
1. 计算输入图像信号的GDFT：X=
2. 将输入图信号的GDFT乘以图传递函数，得到输出谱形式Y = X，和
3.将输出图信号Y进行反离散傅里叶变换(IGDFT)，即
对于大的图形来说，这个过程在计算上要求很高，在大的图形中，直接在顶点域实现所需的滤波器(或其近似值)可能更方便。
对于顶点域的实现，任务是找到(14)中的系数（参照标准脉冲响应），使它们的谱表示H(λ)等于(或近似等于)期望的G(λ)。这是按照以下方式执行的。
顶点域系统的传递函数如(28)所示，并且应该等于期望的传递函数，，由这个条件可以得出了一个线性方程组(33):
（28）
这个方程组的矩阵形式是:
其中Vλ是特征值λk的Vandermonde（范德蒙德）矩阵形式，即：
(36)是系统系数的向量,需要计算得到:
对(33)解的注解：
1. 考虑有N个顶点且具有所有邻接矩阵特征值的情况，最小多项式阶Nm等于顶点数N(换句话说，最小多项式等于特征多项式，Pmin (λ) = P(λ))。
(a)如果滤波器序列数M使M = N，则(33)的解是唯一的，因为Vandermonde矩阵的行列式总是不为零。
(b)如果滤波器序M使M < N，则(33)中的系统是多因素决定的(过定的)。因此，(33)的解只能在最小二乘意义下得到(如本节后面所述)
2. 如果某个特征值的次数大于1(最小多项式阶Nm低于顶点数N)，(33)中的系统简化为Nm个线性方程的系统(通过去除多个重复特征值的方程)。
(a)如果滤波器序M满足Nm<M≤N，(33)中的系统是欠确定的。在这种情况下(M - Nm) 个滤波器系数是自由变量，系统有无限个数的解，而所有这样得到的滤波器都是等价的。
(b)如果滤波器序M =Nm，则(33)中系统的解是唯一的。
(c)如果滤波顺序使M < Nm，则(33)中的系统是超定的，在最小二乘意义下得到解。
3.任何阶为M > Nm的滤波器都有一个唯一的阶为Nm的等价滤波器。该滤波器将自由变量设为0，即i = Nm,Nm + 1,...,N−1时,hi = 0。
最小二乘估计法是对过度确定系统，即其中存在比未知数更多的方程组，以回归分析求得近似解的标准方法。在这整个解决方案中，最小二乘法演算为每一方程式的结果中，将残差平方和的总和最小化。
https://baike.baidu.com/item/%E6%9C%80%E5%B0%8F%E4%BA%8C%E4%B9%98%E4%BC%B0%E8%AE%A1%E6%B3%95/8019852?fr=aladdin
求系统系数
精确解:对于M = N = Nm，即当滤波器的阶数等于顶点数和最小多项式的阶数时，(33)或(34)中系统的解是唯一的，即
最小二乘解:对于超定情况，当M <Nm时，h的均方近似通过最小化平方误差得到：
根据(令误差最小)，可以得出：
由于M < Nm，得到的解是的最小二乘近似。考虑到这个解可能不满足，设计的系数向量（其谱为),服从
通常不同于期望的系统系数g(谱G(Λ))。
pinv：伪逆矩阵是逆矩阵的广义形式。由于奇异矩阵或非方阵的矩阵不存在逆矩阵，但在matlab里可以用函数pinv(A)求其伪逆矩阵。函数返回一个与A的转置矩阵A' 同型的矩阵X，并且满足：AXA=A,XAX=X.此时，称矩阵X为矩阵A的伪逆，也称为广义逆矩阵。
例4:考虑图8(a)中的未加权图，设计合成频率响应为的期望滤波器，该滤波器考虑（33）中M = 1、2、4、6的各种滤波阶数，结果如图13所示。为简洁起见，考虑M = 4时顶点域滤波器为：
然而，精确的频率响应只有在M = N = 8时才能得到。
图13 在频谱域(参考标准频率响应)中具有指定传递函数的图滤波器的设计。期望的谱响应用蓝色的圆表示。红色星号表示例4中设计的滤波器的谱响应，且根据顶点域M个滤波器系数h0, h1，…， hm-1。
clc clear all %%参数设置 A=[0,1,1,1,0,0,0,0;1,0,1,0,1,0,0,0;1,1,0,1,1,0,0,0;1,0,1,0,0,0,1,0;0,1,1,0,0,1,0,1;0,0,0,0,1,0,0,1;0,0,0,1,0,0,0,1;0,0,0,0,1,1,1,0]; %邻接矩阵A [U,lamda]=eig(A); %%求取A的特征值和特征向量 M=4; %滤波器阶数%%查看特征值范德蒙德行列式 x=diag(lamda); %定义5维列向量x for i=1:1:8 %行控制变量i从1~5，步长为1 for j=1:1:M-1 %列控制变量j从1~5，步长为1 lamda_1(i,j)=x(i)^(j-1); %对矩阵元素A(i,j)赋值 end end%% x=x'; y=[0,0,0,0,0,0.5,1,1]; coefficient=polyfit(x,y,M-1); %用一次函数拟合曲线，想用几次函数拟合，就把M设成那个数 y1=polyval(coefficient,x); plot(x,y,'o',x,y1,'-')
3.5.4 图传递函数上系统的多项式(切比雪夫)近似
在不失一般性的情况下，可以认为是期望的传递函数由一个λ连续函数在λmin≤λ≤λmax最大区间内采样组成,所要求的传递函数，期望传递函数G(λ)的变量λ是连续的，并且图上的系统只使用离散点的值，所要求的传递函数，G(λ)的变量λ是连续的，但图上的系统只使用离散点的值。因此，对于期望传递函数G(λ)的多项式近似P(λ)，重要的是在考虑区间λmin≤λ≤λmax内的点的误差有界且足够小。
这个问题在代数上被称为最小-最大逼近问题，它的目标是找到一个与期望函数值有最小的最大绝对误差的逼近多项式。最小-最大多项式可以用截断的Chebyshev多项式P(λ)来近似，得到的函数近似具有几乎最小-最大性质。因此，可以用切比雪夫多项式来近似传递函数G(λ)：
其中，Tm(z)为切比雪夫多项式，定义为:
以变量λ为中心且归一化为：
其中−1≤z≤1(切比雪夫多项式定义要求)。从z到λ的逆映射为：
由于切比雪夫多项式是正交的，测度为，将期望函数G(z)展开成多项式级数的切比雪夫系数cm由下式可得：
例5:考虑图8(a)中的未加权图和期望的传递函数为
G(λ)在离散点上采样,对应于例4图13中输出值。因为最小特征值和最大特征值是λ min =−2和λ max = 3.19,这就得到了在标准化区间−1≤z≤1内变量z的期望传递函数:
其中，z5由（40）可得：
(M−1)= 3时的切比雪夫级数为
通过变量的变换，z→λ，我们得到了这种形式：
图信号滤波现在可以在顶点域使用
其中，
利用上面的顶点域滤波器可以去除图12的噪声，其效果如图15所示，信噪比提高了16.76dB
图14 用带有M项的(M−1)阶切比雪夫多项式逼近，，期望的谱响应G(λ)用蓝色虚线和蓝色点表示。红线表示被设计滤波器的谱响应切比雪夫逼近。
图15
A=[0,1,1,1,0,0,0,0;1,0,1,0,1,0,0,0;1,1,0,1,1,0,0,0;1,0,1,0,0,0,1,0;0,1,1,0,0,1,0,1;0,0,0,0,1,0,0,1;0,0,0,1,0,0,0,1;0,0,0,0,1,1,1,0] [U,lamda]=eig(A) %求矩阵A的特征值和特征向量,其中lamda的对角线元素是特征值,其中U的X的列是相应的特征向量,u6=U(:,7) u7=U(:,8) x=2*u6+3.2*u7 e=[1.65186931,-1.14711105, 0.37797999, 0.81724406, -0.9909782,-0.14525697, -1.18463069, 0.62088355]' x_e=x+e %P=0.07+0.36*A+0.11*A^2-0.04*A^3 %切比雪夫M=4 %P=0.1734+0.3532*A+0.08*A^2-0.0336*A^3 %二乘法M=4 P=0.0478883180835174+0.170777630157730*A+0.176308509911496*A^2+0.0502774460233508*A^3-0.0116121317299550*A^4-0.00549674172663644*A^5 %二乘法M=6 y=P*x_e; %绘制系统输出x的谱图 figure k=0:1:7; stem(k,y)
在标准滤波器设计理论中，可以对近似传递函数G(λ)的过渡区间进行适当的平滑处理，以改进近似。
一般来说，(40)中从λ到z的映射可以写成z = aλ + b，其中a = 2/(λ max - λ min)和b =−(λ max - λ min)/(λ max - λ min)。λ的切比雪夫多项式级数是这样的:
同样的关系也存在：
这种变量的变化允许进行如(39)的递归计算
3.5.5 图上的逆系统
利用G(λ)给出的图上系统H(λ)表示图上系统的逆，可以得到它们的一般关系：
根据(37)，这反过来意味着，如果所有且P(λ) = P min (λ)，则。
3.6 基于拉普拉斯函数的图傅里叶变换
类似于基于邻接矩阵的离散傅里叶变换图的频谱表示,一个图信号的频谱表示可以基于图拉普拉斯特征值分解，即
虽然基于邻接矩阵和图拉普拉斯矩阵的谱分解分析可以统一进行，但由于它们的性质和适用范围不同，因此将分别进行分析。
给出了图傅里叶变换的信号x，它利用图拉普拉斯特征值分解来定义其基函数
其中矩阵U在其列中包含了图拉普拉斯矩阵的特征向量。反图傅里叶变换是这样的形式：
在无向圆无权图的情况下，如图7(a)中的图，这种基于拉普拉斯的频谱分析简化为标准傅里叶变换，但使用实值基函数(特征向量)，如part1 式(32)所示。
3.7 拉普拉斯谱域的排列组合和滤波
如3.5.2节所示，图的移位和邻接矩阵与顶点域的特征向量的第一个有限差有关，而特征向量的变化率与其对应的特征值有关(参考标准频率)。类似的方法可用于基于拉普拉斯的特征分解。我们已经看到，对于标准时域信号，循环图的拉普拉斯图表示信号u(n)的二阶有限差分y(n)，即拉普拉斯矩阵中所有元素的紧凑表达式可以写成矩阵形式y = Lu。现在很明显，特征向量，u表现出小的变化，也应该有一个小的二阶差分的累积能量：
回想一下，这个表达式对应于特征向量u的二次形式，定义为。上述拉普拉斯二次型的推理也适用于图信号。拉普拉斯分析默认情况我们将考虑无向加权图，其中根据定义
这意味着一个特征向量uk的二次形式等于它对应的特征值。
(45)
显然。一个小意味着小的变化，因此，对应于小λk的特征向量对应于图信号的低通部分。也就是说，平滑指数越小，特征向量uk越平滑。
一个理想的拉普拉斯谱低通滤波器，其截止特征值为λc，因此可以定义为:
例6:考虑如图8(a)的无向图上的一个信号，该信号如图16(a)所示。这个图形信号是由两个拉普拉斯特征向量(对应于慢变信号部分)的线性组合产生的，使x = 2u0 +1.5u1。考虑的图的拉普拉斯特征向量在第一部分图13中给出。将原信号x在信噪比为−1.76 dB的情况下被高斯白噪声干扰，如图16(b)所示。接下来，使用理想谱域图滤波器对这个有噪声的图信号进行滤波，截止特征值λc= 2，得到滤波后的信号xf，如图16(c)所示。实现的输出信噪比为信噪比out = 21.29 dB，也就是说，与图16(b)中的有噪信号相比，图滤波器虽然简单，但其信噪比增益为23.05 dB。
为了进一步说明图滤波的原理，我们对图3中的噪声信号使用了一个谱截止值为λc = 0.25的滤波器进行滤波，结果如图17所示。同样的信号也被滤波使用多项式近似低通系统，如图18所示。
图13 选择λk最小的特征向量
图16
正则图的拉普拉斯与基于邻接矩阵的GDFT。对于J-正则无权图，可以建立基于邻接矩阵法和基于拉普拉斯法的谱分解之间的直接关系(见part 1 式(13)）：L = JI − A
可以得到：它们有相同的特征向量。
备注13:特征向量uk从低通到高通的排序，这是基于各自的特征值，基于邻接矩阵法和基于拉普拉斯法图谱分解的顺序完全相反，例如，为得到最平滑的特征向量，要求
或是
3.8 使用图拉普拉斯定义图系统
在第3.2节和公式(14)的讨论之后，图上的系统，使用图拉普拉斯的形式：
对于无权图，图上系统的这个定义可以与相应的邻接矩阵形式L = D−A相联系。然后通过拉普拉斯特征值分解得到系统在图上的谱域描述
我们利用矩阵多项式的特征分解的性质，
因此，可以得出：，以元素形式
在顶点域中，输出信号的第n个元素为：
其中，传递函数定义为：
图脉冲响应是：
注14:(50)中y(n)的表达式可以解释为图上的广义卷积，它是利用顶点域的脉冲响应hn(i)的广义图位移来实现的。接下来，我们通过对单位脉冲的响应来描述图上的广义卷积。
为了说明，考虑位于图顶点m的脉冲函数，由
对应的GDFT为：
它是基于图拉普拉斯特征向量定义的。
观察到，与标准时域相似，任何图信号都可以写成图顶点上的函数之和，也就是：
或是向量的形式：
其中，δi为元素δ(n−i)的向量，则系统输出y可以是以下形式：
其元素为:
式(52)中，hn(i)与H(λk)有关。
:根据(47)，(50)中给出的顶点n的图卷积算子的形式定域在顶点n的(M−1)邻域内。这种局部化性质对于大型图更为重要。
接下来将讨论两个任意图信号的广义卷积。
3.9 图上信号的卷积
考虑两个图信号，x(n)和h(n)。利用这两个信号的图拉普拉斯谱[18]，基于图上卷积谱的假设，定义了图上这两个信号的广义卷积算子:
等于对应的图形信号x(n)和h(n)频谱的乘积，以元素形式即：
以元素形。因此，广义图卷积运算的输出，x(n)∗h(n)，等于(55)中的谱乘积Y (k)的反离散傅里叶变换，即
注意在(56)中的H(k)和(51)中的H(λk)的定义之间的区别。下一节将更详细地讨论这两种形式。
在图上移位——另一种定义
上述广义图卷积框架也可作为方便定义图上位移的基础。考虑图像信号h(n)和δ函数位于一个顶点m。在这里,我们将使用hm(n)表示的指向顶点m的图信号h (n)移动版本。这种转移信号将根据经典信号处理中的推理，移位的信号是原始信号和适当移位的脉冲函数的卷积。因此，这里通过广义图卷积定义一个图移信号，根据(54)和(55)，其GDFT等于。那么图移信号就是的IGDFT，即
在计算X(k)H(k)的反离散傅里叶变换时，也会出现同样的关系：
上式(59)是图移信号的另一种形式。由于H(k)定义为信号H(n)的GDFT与(52)的定义不同，产生了不同的移位操作，分别表示为
备注16:注意这两种移位操作，(52)或(59)满足这样一个性质即移动0等于原始信号，即
例7:考虑图8(a)中的图上的一个信号，它由其图的拉普拉斯GDFT定义，给出 τ = 0.1573.用(59)中的移位算子得到的所有移位信号如图19所示。
图19 图移位算子的一个例子。顶图为由其拉普拉斯GDFT定义的图信号，即图信号hm(n)为在m = 0时“移位”到m = 7的移位信号，用(59)计算。
clc clear all h=zeros(8,8); A=[0,1,1,1,0,0,0,0;1,0,1,0,1,0,0,0;1,1,0,1,1,0,0,0;1,0,1,0,0,0,1,0;0,1,1,0,0,1,0,1;0,0,0,0,1,0,0,1;0,0,0,1,0,0,0,1;0,0,0,0,1,1,1,0]; D=[3,0,0,0,0,0,0,0;0,3,0,0,0,0,0,0;0,0,4,0,0,0,0,0;0,0,0,3,0,0,0,0;0,0,0,0,4,0,0,0;0,0,0,0,0,2,0,0;0,0,0,0,0,0,2,0;0,0,0,0,0,0,0,3]; L=D-A; [U,lamda]=eig(L); LAMDA=diag(lamda); H=exp(-2*LAMDA*0.1573); for m=1:1:8 %行控制变量m从1~8，步长为1 for n=1:1:8 %列控制变量n从1~8，步长为1 for k=1:1:8 %列控制变量k从1~8，步长为1 h(m,n)=h(m,n)+H(k,1)*U(m,k)*U(n,k); %对矩阵元素h(m,n)赋值 end end end for i=1:1:8subplot(4,2,i)k=0:1:7;stem(k,h(i,:),'r')axis([-0.5 7.5 0 0.6]);title(['m=' num2str(i-1)]) box off end
3.10 信号在图上的z变换
图信号移动算子的关系，可以用和的定义来建立。考虑(51)是信号h(n)的离散傅里叶变换。即
当系统系数hn, n = 0,1，…，M−1与H(λ k)的关系为(51)，即
对于M = N，后两个关系的向量形式为：
所以信号h(n)和系数hn的关系，可以用表示
注17:在经典的DFT中(有向圆图及其邻接矩阵的情况，用代替)，信号样本h(n),可由H(λk)的逆DFT得到或是系统系数hn得到，因为特征值等于相应的谱域移位算子，即
因此，对于经典的DFT分析，以下关系成立
这个关系在(35)和中是明显的，用来定义一个图信号的z变换
图信号的z变换。对于给定的图信号，根据(60)中的解释，系统的系数，对应于具有与图信号本身相同的GDFT的系统传递函数：
因此信号x的图像z变换等于系数的经典z变换：
所以下面的结论成立：
输出信号y(n)现在可以由下得到：
其中输出图信号y(n)是系数yn的反z变换的结果，即
在无向图邻接矩阵的分解中，当特征值为复值时，在复值z域的z变换表示可能是有意义的。以第一部分图1(b)的图及其邻接矩阵为例，特征值如图20所示。
定义:将解析图信号和图Hilbert变换在谱域内定义为：
式中为λk的虚部。如果将这些关系式应用到的标准DFT中，就可以得到相应的经典信号处理定义。
3.11 频谱域的移位算子
频谱域的移位操作可以用顶点域的移位相同的方式定义。考虑定义在无向图上的两个图信号x(n)y(n)的乘积。这个积的GDFT是:
可以认为是Y(k)对i个谱指数的移位。
备注18:按期想中，在谱域中i = 0的移动为原始值，直到常数因子。这个关系对于顶点域的移动算子不成立。
3.12 图上的Parseval定理
考虑在无向图上观测到的两个图信号x(n)和y(n)及其谱X(k)和Y(k)，则Parseval定理具有如下形式
它对任意两个图信号都成立。
要在图上证明Parseval定理，需考虑
在原始域和谱域的能量之间产生Parseval等价。假设图是无向的，因此U−1 = U T成立。这个定理是相当普遍的，适用于基于无向图分解的拉普拉斯矩阵与邻接矩阵。
3.13 最优的去噪
考虑一个测量值x，由一个慢变化的图形信号s和一个快速变化的扰动ε组成：x = s + ε
目的是设计一种抑制干扰的滤波器(去噪)，其输出表示为y = H(x)。
最优去噪任务可以定义为目标函数的最小化：
物理上，第一项的最小值迫使输出信号y尽可能接近可用的观测值x，根据它们的欧式距离的能量(最小误差能量)，而第二项表示y的平滑度度量(见第3.7节)，这在物理上也有意义，因为原始输入的s是低通的，比扰动ε更平滑。参数α提供了输出y与x的接近度和输出平滑度准则之间的平衡度。
为了解决这个最小化问题，我们求导：
可得到：，这个关系的谱域形式由L=可以得到
然后，上述谱输入/输出关系的元素传递函数的形式为：
注19:对于小的α，我们有H(λk)≈1，即(65)的全通行为，不进行信号平滑，得到y≈x。H(λk)≈δ(k)。结果y≈const。是最大平滑的(一个恒定的输出，没有任何变化)。
例8:用(65)中α = 1的最优滤波器对图3中的含噪信号进行滤波，结果如图21所示。获得的信噪比为19.16 dB。
其他成本函数。在许多可能的替代方案中，我们将引入两个额外的代价函数来进行图信号去噪，利用不同的约束强加在解决方案上。
为了避免使输出信号平滑，我们可能希望它与线性形式(即满足Ly = 0的信号y)的偏差尽可能小，这可以通过给定的代价函数来实现
产生了一个封闭形式的去噪解决方案与相应元素的谱域关系
例如，(64)和(66)中的两种成本函数形式的组合可以在滤波器传递函数的设计中提供额外的灵活性
产生传递函数为：
该传递函数形式可以通过参数α和β的选择进一步微调。例如，如果我们希望的分量不衰减，我们可以使用α+βλ1 = 0。这种成本函数可以直接推广为M个未衰减分量的传递函数。
稀疏促进解决方案。一些应用要求提高输出图信号的稀疏性，而不是平滑性。这样的解决方案自然依赖于压缩感知理论，该理论要求用促进稀疏性的规范取代先前成本函数中的两范数。这类成本函数的两个例子是
备注20:当p = 0时，零范数在促进稀疏性方面是最好的，因为当p = 0时，(67)中的成本函数的第二项计数(并最小化)了Ly中非零元素的数量。最小化Ly的稀疏性促进y的常数(或线性)解，具有最小数量的不连续点(向量Ly的非零元素)。在(68)中的第二个代价函数中，零范数促进了非零元素的最小可能数目。这也被称为全变差(TV)方法 the total variation (TV) approach。然而，这些目标函数的最小化不能通过分析的方式实现，例如在p = 2的标准MSE情况下。
另一方面，p = 1的一范数使上述代价函数为凸函数，允许使用梯度下降法求解，同时在一些缓和的条件下产生与p = 0时相同的解。“1范式”可以作为“0 范式”的分析替换
3.14 用随机漫步拉普拉斯算子定义的图上的系统
4 分采样、压缩感知和重构(Subsampling, Compressed Sensing, and Reconstruction)
图可以由大量的顶点组成，其数量级可达数百万甚至更高。相关的计算和存储问题带来了对图上定义的分采样和压缩感知的潜在优势的考虑。这里，我们提出了几种基本的子采样方法，以及它们与经典信号处理的关系。
4.1 低通图信号的分采样
为了方便起见，我们将从最简单的情况开始即考虑的图信号是低通性质的。
这样的信号可以表示为具有最低平滑指数的图拉普拉斯的K < N个特征向量线性组合：
该(k -稀疏)信号在GDFT域中的GDFT域系数为如下形式：
回顾可知，如果k<<N,图信号在GDFT域是稀疏的
因此，恢复稀疏信号所需的最小图信号样本数M的要求是M=K<N。为了重构的稳定性，通常采用K≤M < N个图信号样本。可用的图信号样本向量称为测量向量，用y表示，而图信号样本可用的顶点集合(V ={0,1,2，…，N−1}的随机子集)用表示
度量矩阵现在可以使用IGDFT中的x = UX，其元素形式由(76)给出。(76)式中对应顶点n∈M = {n1,n2，…，nM}，然后定义系统
其中其中是度量矩阵，是测量向量并由可用的图形信号样本组成。一般来说，由于M < N，这个方程组是欠定的，在没有附加约束的情况下，不能唯一地求解X。
假设信号的频谱表示仅包含K≤M个慢变特征向量的线性组合，我们可以排除(77)中的GDFT系数X(K)，X(K + 1)，…，X(N−1)，因为它们是零值的，对图信号样本的生成没有贡献。考虑到这一点，式(78)中的M × N方程组简化为如下M × K方程组:
若是M = K个独立的测量值，该系统可以被唯一地求解，而对于M > K，该系统是典型的超确定的，在最小二乘(LS)意义上找到解，如：
计算后，所有的GDFT值直接按代入，零值写入X(K)， X(K + 1)，…X (N−1)。然后使用x = UX在所有顶点恢复图形信号。
恢复状态。(80)中的信号重构是可能的，如果逆存在，意味着就矩阵条件数而言，这个要求等价于
就矩阵条件数而言，这个要求等价于
备注21:对于图信号的噪声测量，重构的GDFT系数中的噪声与输入噪声和矩阵条件数直接相关。
如果我们能够选择可用的信号样本位置(顶点)，那么采样策略将是找到一组测量值，使它们产生尽可能接近统一的条件数(为了稳定和减少噪音的影响)。
例10:为了演示由简化的图信号样本集重构的原理，考虑一个图信号在M = 3个顶点处的值，给定如图23上图所示。假设图信号是低通类型，当K = 2时，最小非零GDFT系数为X(0)和X(1)。然后可以从(82)重构这个图信号的GDFT系数:
根据(76)中的定义，假设在n = 0, n = 2，和n = 6的三个顶点，对于两个非零系数X(0)和X(1)，可用信号样本x(n)为:
矩阵A32的秩是2,对应的矩阵条件数为而重构的GDFT非零值为X(0) = 2和X(1) = 1，得到重构图信号x = UX(其中),如图23底图所示。
图23:低通图信号的下采样Subsampling图示。上图:在顶点1、3、4、5和7处缺少样本的图形信号；下图:重构图信号。
clc clear all h=zeros(8,8); A=[0,1,1,1,0,0,0,0;1,0,1,0,1,0,0,0;1,1,0,1,1,0,0,0;1,0,1,0,0,0,1,0;0,1,1,0,0,1,0,1;0,0,0,0,1,0,0,1;0,0,0,1,0,0,0,1;0,0,0,0,1,1,1,0]; D=[3,0,0,0,0,0,0,0;0,3,0,0,0,0,0,0;0,0,4,0,0,0,0,0;0,0,0,3,0,0,0,0;0,0,0,0,4,0,0,0;0,0,0,0,0,2,0,0;0,0,0,0,0,0,2,0;0,0,0,0,0,0,0,3]; L=D-A; [U,lamda]=eig(L); LAMDA=diag(lamda); y=[1.140,0.996,0.563]' k1=[1,3,7] ; k2=[1,2] ; %对应节点0，2，6 for m=1:1:3 %行控制变量m从1~3，步长为1 for n=1:1:2 %列控制变量n从1~2，步长为1 Amn(m,n)=U(k1(m),k2(n)); %对矩阵元素A(m,n)赋值,取对应节点0，2，6的特征向量u0.u1 end end r=rank(Amn); %矩阵Amn的秩 cond_Amn=cond(Amn'*Amn); %矩阵Amn'*Amn的条件数 Xk=((Amn'*Amn)^(-1))*Amn'*y; %Xk=pinv(Amn'*Amn)*y X=zeros(8,1) for m=1:1:2 %行控制变量m从1~2，步长为1 for n=1:1:1 %列控制变量n从1~1，步长为1 X(m,n)=Xk(m,n); %将求得的Xk写入X对应的位置，其余位置为0 end end x=U*X %重构信号%绘图 k=0:1:7; stem(k,x,'r') title('重构信号' )
注22:对于特征向量的有向循环图，上述采样和插值的关系与经典信号处理中相同。
4.2 稀疏图信号的下采样Subsampling
其次，在非零GDFT系数位置已知和未知的情况下，将考虑在GDFT域中稀疏的图信号的下采样。
4.2.1 系数在GDFT中的位置已知
前面的分析不仅适用于低通类型的图信号x及其对应的GDFT X，而且适用于任意已知谱位置的K个非零值的GDFT X，即:
与(78)类似，对应的方程组为：
矩阵形式为用于求解非零谱值X(k)，，与4.1节给出的低通信号相同
4.2.2 支持矩阵，下采样和上采样(Support Matrices, Subsampling and Upsampling)
在信号处理的过程中，下采样问题通常被称为支持矩阵。假设一个图信号x的下采样方式是在顶点的子集上可用，而不是全部的顶点集合。对于这个下采样信号，我们可以定义它的上采样版本Xs，在信号不可用的节点处加0。利用数学形式，原始信号x的下采样和上采样版本，xs为：
这里的支持矩阵B是一个N × N的对角矩阵，矩阵中1在对角的位置对应
和0在矩阵其他地方。得到信号x的下采样和上采样版本xs是这样一种方式:将信号x在一个简化的顶点集上进行下采样，然后在未定义下采样信号的原始信号位置加零进行上采样。
回顾一下，在没有附加约束的情况下，一般情况下，一个具有N个独立值的信号x不能从它在xs中的M < N个非零值重构出来。然而，对于在GDFT域中也是稀疏的图信号，额外的约束是信号x只有在GDFT域中的K≤M个非零系数，X=,,所以拥有以下关系：
X = CX
其中，支持矩阵C是一个N × N的对角矩阵，其对角位置为1对应于和其他位置为0。注意，在这个方程的两边都存在GDFT X，这与(84)中的x不同。重构公式如下：
因为，所以可得：如果BUC的秩是K(如果有K个线性无关的方程)，可能有K个CX的非零系数，也就是说
这个条件等价于(81)，因为矩阵BUC的非零部分等于(83)中的
4.2.3 系数的位置未知
当非零谱系数的位置是未知 时，重构问题更为复杂;这种情况已经在标准压缩感知理论中得到解决，可以表述为
其中是X的非零元素的个数（零范数）
虽然解决这个最小化问题的方法是多种多样的，我们在这里采用一个简单的两步方法:
1）使用M > K个信号样本的估计非零系数位置
2）在估计位置K重建X的非零系数，以及在所有顶点上的信号x，使用已知非零系数位置的重建方法，如第4.1和4.2.1节。位置K的非零系数计算为。
第1步中GDFT的非零位置可以通过测量值(可用信号样本)y在测量矩阵上的投影来估计
X0中K个最大值的位置被用作非零位置的估计。这个过程也可以通过迭代的方式实现：
(i)在第一次迭代中，我们假设K = 1，然后继续估计图信号中的最大谱分量。在确定其位置为之后，非零位置的初始空集合变为。然后将重构向量，(其中)，从测量值y中移除。在这种情况下，矩阵AM1是由索引k1定义的矩阵AMN的一列，e = y−y1的差值作为下一步的测量向量
(ii)通过求解k2 = 来估计图信号中第二大谱分量的位置。非零位置的集合现在变成。图信号的第一个和第二个分量现在估计为，其中矩阵AM2是测量矩阵AMN的子矩阵，AMN由指标k1和k2定义的列组成。重构向量y2 = A2X2，从测量值y中移除，误差e = y−y2，现在作为新的测量向量
(iii)迭代重复K次，直到e中剩余的测量值可以忽略。在稀疏性K未知的情况下，迭代过程直到，其中是预定义的精度。
例11:考虑一个稀疏图信号，稀疏度K = 2，在顶点n = 2,3,4,5,7上测量，取值如图24(上)所示。我们的任务是重建完整的信号，即寻找缺失的样本x(0),x(1), x(6)。
根据(86)对给定的测量值y进行计算非零元素在GDFT中的估计位置以及初始估计值X0，因为K = 2，两个非零系数的位置估计为X0中两个最大的值的位置。在考虑的情况下，，如图24(下所示。然后在稀疏度K = 2的情况下重构GDFT系数，如X2 =，得到X(0) = 2, X(3) = 1.2，如图24(右下)所示。最后，所有顶点的重构图信号x = UX显示在的图24中间。
clc clear all h=zeros(8,8); A=[0,1,1,1,0,0,0,0;1,0,1,0,1,0,0,0;1,1,0,1,1,0,0,0;1,0,1,0,0,0,1,0;0,1,1,0,0,1,0,1;0,0,0,0,1,0,0,1;0,0,0,1,0,0,0,1;0,0,0,0,1,1,1,0]; D=[3,0,0,0,0,0,0,0;0,3,0,0,0,0,0,0;0,0,4,0,0,0,0,0;0,0,0,3,0,0,0,0;0,0,0,0,4,0,0,0;0,0,0,0,0,2,0,0;0,0,0,0,0,0,2,0;0,0,0,0,0,0,0,3]; L=D-A; [U,lamda]=eig(L); LAMDA=diag(lamda); y=[0.707,1.307,0.407,1.307,0.407]' k1=[3,4,5,6,8] ; k2=[1,2,3,4,5,6,7,8] ; %对应节点2,3,4,5,7 for m=1:1:5 %行控制变量m从1~3，步长为1 for n=1:1:8 %列控制变量n从1~2，步长为1 Amn(m,n)=U(k1(m),k2(n)); %对矩阵元素A(m,n)赋值,取对应节点0，2，6的特征向量u0.u1 end endX0=Amn'*y X=zeros(8,1) k1=[1,4] ; k2=[1,2,3,4,5] ; %对应节点2,3,4,5,7 for m=1:1:5 %行控制变量m从1~3，步长为1 for n=1:1:2 %列控制变量n从1~2，步长为1 Am2(m,n)=Amn(k2(m),k1(n)); %对矩阵元素A(m,n)赋值,取对应节点0，2，6的特征向量u0.u1 end end X2=pinv(Am2)*y X(1,1)=X2(1); %将求得的Xk写入X对应的位置，其余位置为0 X(4,1)=X2(2); %将求得的Xk写入X对应的位置，其余位置为0 x=U*X %重构信号k=0:1:7; title('重构信号' ) stem(k,x,'r')
4.2.4 独特的重建条件
与标准压缩感知问题一样，当，初始GDFT估计X0将产生非零元素X(k)的正确位置，重构将是唯一的，其中µ等于测量矩阵AMN任意两列间内积的最大值，
(µ被称为一致性指数 )
为了说明重构的唯一性，回顾一下k -稀疏信号可以写成
其中(86)的初始估计等于X0 =，其元素为：
Graph Signal Processing – Part II: Processing and Analyzing Signals on Graphs(文献翻译)相关推荐
- GNN笔记：图信号处理（Graph Signal Processing）
1 图信号处理定义 图信号处理(Graph Signal Processing,以下简称 GSP)用来处理那些定义在图上的非规则域的信号.换句话说,就是处理图上定义的信号,但信号所在域是非规则的. 2 ...
- 图信号处理基础------Foundations in Graph Signal Processing
本篇博文意在用例子的形式解释图信号处理基础,目的是记录总结自己的学习过程,有兴趣的读者也可也一起交流. 前言 图是一种包含多种数据属性的不规则结构.然而,传统的图论处理方法都注重分析底层的图结构而不是 ...
- Processing入门教程-Processing的“前世今生”
很早以前大概13.14年就通过清华大学付志勇教授了解到了Processing这个工具,起初只是初步了解并没有下定决心学习(当初资料太少了).由于当时只是初步的看了看,所以很多内容和知识点都是一知半解的 ...
- MULTI-CHANNEL SPEECH ENHANCEMENT USING GRAPH NEURAL NETWORKS 文献翻译
MULTI-CHANNEL SPEECH ENHANCEMENT USING GRAPH NEURAL NETWORKS 文献翻译 来自于脸书实验室的一篇文章,将图神经网络用在了多通道语音增强上面,思 ...
- Graph Signal Processing——Part I: Graphs, Graph Spectra, and Spectral Clustering （文献翻译）
目录 目录 0.Abstract 1.Introduction 2.图形定义和属性 2.1 基本定义 2.2 一些常用的图拓扑 2.3 图及其相关矩阵的性质 3.图矩阵的谱分解(特征分解) 3.1 邻 ...
- 【文献翻译】Evaluating five different adaptive decomposition methods for EEG signal seizure detection
五种不同的自适应分解方法在脑电信号癫痫检测和分类中的应用 文章目录 摘要 1 引言 2 - 方法 2.1 数据集 2.1.1. 波恩大学(UOB)数据集 2.1.2. 新德里神经病学和睡眠中心(NSC ...
- 【多目标轨迹预测】Path-Aware Graph Attention for HD Maps in Motion Prediction（ICRL 2022）翻译
参考: (31条消息) [文献阅读]路径感知的图注意力做运动预测(Fang Da等人,ArXiv,ICRA 2022)_全部梭哈迟早暴富的博客-CSDN博客_lanegcn lanercnn ICRA ...
- （3）小波变换原理及应用
Part1-Introduction To The Wavelet Transform(简介) 1.Origin of the wavelet transform: The theories of W ...
- Signal Processing Toolbox
MathWorks官网是个宝,没事就去逛逛兴许就能学到东西.MATLAB用户很多很多,但会去逛MathWorks官网的人并不多.官网有很多好东西,只是上面的内容实在太多了,在不熟悉官网架构的情况下要找 ...
最新文章
- java 16 -12 静态导入
- Console命令详解,让调试js代码变得更简单
- 小Z的房间[HEOI2015] （matrix-tree定理）
- svn文件大小类型限制，提交必须加多少字的说明
- 安卓手机老是自动保存图片_Redmi K30 Pro自动亮度调节和iPhone基本一致，安卓手机的大进步...
- 线性模型第3讲：Lasso方法
- 修改pip install镜像源
- Android蓝牙开锁讲解
- Activity启动模式singleTask模式
- linux打开txt文件命令_linux系统文件及常用命令
- 嗖嗖移动业务大厅项目_会员合作项目：10086移动外呼业务
- 火车票查询软件测试自学,火车票订购系统的测试报告.doc
- orbslam 2 运行 tum 数据集中的 walking xyz 序列
- macd的顶背离和底背离
- c# forbidden.html,C＃Web API方法返回403 Forbidden
- oracle order siblings by,sql中ORDER SIBLINGS BY排序的含义
- python用matplotlib画五角星_3.用Python画五角星
- 如何使用ansible管理多台远程服务器
- java中 ^ 是什么意思
- 一个能永久存储网页快照的网站