第六章功率谱估计电子教案.doc

第六章 功率谱估计的经典方法6.1 引言 从第二章的讨论中，我们已经知道一个随机信号在各时间点上的值是不能先验确定的，它的每个实现(样本)往往是不同的，因此无法象确定信号那样可以用数学表达式或图表精确地表示它，而只能用它的各种统计平均量来表征它。其中，自相关量作为时移的函数是最能较完整地表征它的特定统计平均量值。而一个随机信号的功率谱密度(函数)，则是自相关函数的傅立叶变换。对于一个随机信号来讲，其能量通常为无限大，它本身的傅立叶变换是不存在的，常常需要研究其功率在频域上的分布。因此，功率谱密度是一个随机信号的一种最重要的表征形式。我们要在统计意义下了解一个随机信号，就要求知道(或估计)它的功率谱密度。 如果我们用表示随机信号的自相关函数，表示它的功率谱密度(以下简写成PSD)，则有： (6.1)而其中(6.2)即为滞后积的数学期望。 根据各态历经假设，零均值广义平稳随机过程的集合的平均可以用一个样本序列的时间的平均代替，于是上式可写成 (6.3)实际上，首先不可能获得样本序列的所有数据，即无数个x(n)，其次通常检测到的数据都是含有噪声或者干扰。因此，只能根据有限个含有噪声的检测数据估计随机信号的自相关序列，进而估计功率谱。将式(6.3)代入式(6.1)得 令，则，上式可写成 (6.4)式(6.4)在的极限情况下是不可能收敛的，这是因为对于无限时域的随机信号，它的傅氏变换是不存在的。实际上只有将式(6.4)求平均，成为 (6.5)才有意义。以后我们还将会看到，只有将式(6.4)经过求平均(或平滑)，即只有式(6.5)才能满足一个正确的估计必须满足的一致估计的条件。 由于实际得到的随机信号只能是它的一个实现或一个样本序列的片段，因此问题是如何根据它的有限个样本序列来估计信号的自相关函数或功率谱密度。这是本章要讨论的中心内容。当我们用一个样本的记录的有限个数据来估计自相关函数和功率谱密度时，有(6.6)(6.7)这里(6.8)wN为矩形函数，或按式(6.4)(6.9)这里是有限长序列的傅氏变换。 一个好的估计应该是无偏估计，最小方差估计。如果我们用表示某个随机变量的真值，表示它的估计值，则希望满足： (1) 无偏估计 无偏估计，即的偏差(Bias)为零，所谓偏差(用B表示)定义为无偏估计即B0，的估计。图6.1中的估计1和估计2都属于无偏估计。 (2) 最小方差估计 最小方差估计，即方差为最小的估计。图6.1中，估计2较之估计1方差小。图6.1 二种估计的概率密度分布 但是常常发生这种情况，一种估计的偏倚较小，而方差较大；另一种估计偏倚较大而方差较小；此时很难定哪一种估计好。因此也常常用均方误差的大小来衡量估计的优劣。在第二章中我们已经讨论到均方误差定义为不难证明均方误差为与偏差和方差均有关，要最小就要求B2与之和最小。 由于，当时式(6.6)就成为式(6.3)。因此应有这就是说，当观察到的样本的数据有无限多个时，按照无穷多个这样的样本数据估计到的自相关函数应该就是自相关函数的真值(各态历经假设)。换句话说，一个正确的估计应有及(6.10)满足式(6.10)的估计称为一致估计。一个正确的估计应该满足一致估计的条件(这是正确估计的必要条件，不是充分条件)。反之，如果某种估计方法不能满足式(6.10)一致估计的条件，则这种估计方法一定是不正确的。下面我们在讨论各种估计方法时，常常以此作为估计正确与否的主要准则之一。 功率谱估计有着极其广泛的应用，不仅在认识一个随机信号时，需要估计它的功率谱。它还被广泛地应用于各种信号处理中。下面我们举三个应用的例子。 在信号处理的许多场所，要求预先知道信号的功率谱密度(或自相关函数)。例如，在最佳线性过滤问题中，要设计一个维纳滤波器就首先要求知道(或估计出)信号与噪声的功率谱密度(或自相关函数)。根据信号与噪声的功率谱(或)才能设计出能够尽量不失真的重现信号，而把噪声最大限度抑制的维纳滤波器(见第二章)。 常常利用功率谱估计来得到线性系统的参数估计。例如，当我们要了解某一系统的幅频特性时，可用一白色噪声通过该系统。再从该系统的输出样本y(n)估计功率谱密度。由于白色噪声的PSD(用表示)为一常数即，于是有 故通过估计输出信号的PSD，可以估计出系统的频率特性(模特性)。 在第七章将要讨论到用自回归模型法估计PSD的一节中，我们将要具体讨论系统参数估计与PSD估计间的关系。 从宽带噪声中检测窄带信号。这是功率谱估计在信号处理中的一个重要用途。但是这要求功率谱估计有足够好的频率的分辨率，否则就不一定能够清楚地检测出来。所谓谱估计的分辨率可以粗略地定义为能够分辨出的二个分立的谱分量间的最小频率间隙(距)。 谱估计问题无论从认识一个随机信号或从其它应用方面来讲都是重要的。因此对谱估计方法的研究引起了国内外学者的广泛注意与重视。它是当前在信号处理中的一个十分活跃的课题。功率谱估计总的来讲可以分为经典谱估计方法和现代谱估计方法。经典谱估计为线性、非参数化方法，需采用经典的傅里叶变换及窗口截断。经典谱估计方法包括周期图法，相关图法等，对长序列有良好估计。现代谱估计为非线性、参数化方法，包括最大似然估计，最大熵法，AR模型法，预测滤波器法，ARMA模型等。对短序列的估计精度高，与经典法相互补充。现代谱估计是融合经典变换理论、统计估计理论、系统辨识、信息论、时间序列分析及计算方法等理论与技术的一门新学科。目前应用广泛，发展迅速。功率谱估计的经典方法以傅立叶变换为基础，很早就被提出，在FFT算法出现后得到了广泛的应用。本章内容为经典谱估计方法。6.2 自相关序列的估计功率谱的估计要求计算自相关序列，下面讨论自相关序列的估计方法。6.2.1 自相关序列的无偏估计设观察到N个样本序列的值：。现在要由此N个数据来估计自相关函数。由于只能观察到的N个值，而时的值是不知道的，那么滞后积序列是一个长为N|m|的序列，因此，式(6.6)成为(6.11)式中m取绝对值是因为，m为负值时上式仍适用。式(6.11)规定的求和上下限的原则是保持充分利用全部(N个)数据。这种估计方法的效果如何，我们首先需要看它的偏差与方差是否满足一致估计的条件。 由式(6.11)，得 上面Rxs(m)是自相关函数的真实值。所以故这种估计，当时，属于无偏估计。 现在来求，按定义 (6.12)又按式(6.11)(6.13)当随机序列是零均值，实、高斯序列时，有所以 代入式(6.13)，得代入式(6.12)，得令，显然的最小值为，最大值为，且的情况将出现次，的情况将出现次以此类推，对不同值的情况，出现的次数将为，于是上式可写成 (6.14)当Nm时，上式以1/ N趋于零，即。故的方差满足一致估计的条件。如果不是高斯过程，在上式中需要再加一项，但此项往往是可以忽略的，因此，式(6.14)仍近似适用。式(6.11)这种估计自相关函数的方法，虽然当m很小于N时能得到一致估计，但当N一定，|m|接近于N时，即时，就变得很大，因而不能得到有用的估计。6.2.1 自相关序列的有偏估计因此很多学者如Jenkins-Watts和Parzen等都主张按下式估计，我们用表示 (6.15)同时这相当于它的均值等于真值用三角窗函数加权。故是有偏的，其偏差为因此， 是的渐进无偏估计。同时 （） 事实上，将用三角窗函数加权后，不仅使方差减小，而且有利于钝化的截断边界，从而改进对的估计。 讨论：虽然的Bias和Var都不等于零，但当时，说明是的渐进无偏估计和有效估计，因此，是的一致估计。且，同时可以证明的均方误差小于的，以及可改进对的估计，所以今后我们还是用作为自相关函数的估计，即(6.16) 此外，m越大，分辨率越高，但是自相关序列估计的偏差也相应增大，通常取，效果较好。通过将自相关函数的估计进行傅氏变换求得功率谱估计的方法即为BT PSD法。6.3 周期图作为功率谱的估计功率谱估计的经典法实质上就是传统的傅立叶分析法，它又可分成二种。一种是间接法，它先通过式(6.6)对自相关函数进行估计(一般都需要窗函数将自相关值加权，以减小自相关序列截断的影响)，然后再通过式(6.7)作傅立叶变换得功率谱估计值。这种方法是1958年由Blackman与Tukey提出的，简称BT PSD估计法，也称为相关图法。另一种是直接法，它是将观察到的有限个样本数据利用FFT算法作傅立叶变换直接按式(6.9)进行功率谱估计(不通过自相关函数的估计)，这种方法称为周期图法。本章的经典法中主要讨论周期图法。 用周期图(包括平滑后的周期图)作为功率谱估计的方法可利用FFT进行计算，因而有计算效率高的优点，在谱分辨率要求不高的地方常用这种周期图法进行谱估计。它的一个主要缺点是频率分辨率低。这是由于周期图法在计算中把观察到的有限长的N个数据以外的数据认为是零。这显然与事实不符合。把观察不到的值认为是零，相当于将x(n)在时域里乘上了一个矩形窗口函数，在频域里相当于与一个sinc函数卷积，由于sinc函数与函数比较有二方面的差别，其一是其主瓣不是无限窄，其二是有旁瓣，因此卷积的结果必然造成失真。由于主瓣不是无限窄，如果原来真实的功率谱是窄的，与主瓣卷积后将使功率向附近频域扩散，使信号模糊，降低了分辨率，主瓣愈宽分辨率愈差。由于存在旁瓣，又将产生两个后果，一是PSD主瓣内的能量“泄漏”到旁瓣将使谱估计的方差增大，二是与旁瓣卷积得到的信号功率谱完全属于干扰。在严重的情况下，强信号与旁瓣的卷积可以大于弱信号与主瓣的卷积，使弱信号淹没在干扰的强信号中，而无法检测出来。这正是周期图作为功率谱估计的二个主要缺点。 对于BT PSD估计法，由于它是按式(6.7)将功率谱用有限个自相关函数值的傅氏变换代表按式(6.1)用无限个自相关函数值的傅氏变换求得的真实功率谱。这相当于将序列乘了一个矩形窗函数，因此也同样存在上述二个缺点。许多学者想尽办法，企图选择适当窗口函数的形状来提高经典法的谱分辨率，但是发现，所有能降低旁瓣的窗口函数都是以主瓣的增宽为代价的；反之亦然。这两个缺点只能互换而无法同时改善。因此，近几年来，在提高功率谱估计的分辨率方面提出了很多新的方法，其中包括以1967年Burg提出的最大熵谱分析法为代表的现代谱估计法，现代谱估计将在下一章阐述。 6.3.1 周期图 将式(6.16)代入(6.7)得 这里的下标N表示它们是有限长(长为N)的序列，令l=n+m，有 其中的为即是有限长序列x(n)的傅氏变换。显然是周期性的。直接将的模的平方除以N求得的功率谱的估计称为周期图，并用表示。于是有(6.17)如果我们观察到x(n)的N个值：，可以通过FFT直接求得频率离散化的。然后按照式(6.17)直接求得(不必先通过估计自相关函数)。这种将周期图作为功率谱估计的方法的主要优点是计算方便，它可以直接用FFT算法从x(n)得到。从而得到。正是由于它的这个优点，使这种方法成为一种十分通用的方法。6.3.2 周期图的估计性能 为了了解周期图作为功率谱估计的估计效果，让我们来讨论它的偏差和方差。为此首先需要求周期图的期望值。 按式(6.16)，得 (6.18)这里wN代表矩形序列。令(6.19)由于是二个矩形函数的卷积，因此它必是一个三角函数，常常称它为Bartlett窗函数，用表示，不难证明而它的傅氏变换为(6.20.0)图6.1所示为Bartlett函数及其傅立叶变换示意图。图6.2 Bartlett函数及其傅立叶变换示意图又自相关函数的真值将上式与式(6.19)代入式(6.18)，并求其傅氏变换，得 (6.20)由式(6.20)可见，除非为函数，将不等于，故周期图作为的估计是有偏的：当时故因此，周期图作为功率谱的估计，当时是无偏的。是渐进无偏估计。 其次求周期图的方差。为了得到周期图的方差，首先假设序列是一个实，白色，零均值过程的样本，具有高斯概率分布函数。按方差的定义应有(6.21)按式(6.17)，周期图可表示为 将 和代入上式，并求期望值：(6.22)其中 由于我们假设是白色的，所以代入上式得 为了求得，按式(6.21)除了需要求得以外，还需要求得，为此我们先求在频率为的协方差，参照式(6.22)可得 (6.23)对于白色高斯过程可以证明矩分解定理考虑到对于白色过程有所以代入式(6.23)得 (6.24)上式当n=k及p=q时只有第一项存在(其它二项为0)；当n=p及k=q时只有第二项存在；当n=q，k=p时，只有第三项存在。又因为故式(6.24)成为(6.25)而其中同理可得代入式(6.25)得当时得 因为上式方括号中只有一项不为0，故有 由公式可见，当，。这说明周期图不满足一致估计的条件。不论N取多长，都有的量级。因此周期图不是对功率谱的好的估计。实际上当时式(6.17)成为式(6.4)。正像我们在讨论式(6.4)时所指出的，对于无限能量的随机序列，它的傅氏变换是不存在的，因此式(6.4)在的极限情况下是不可能收敛的。故也就不能期望当时会等于它的真值，而满足一致估计的条件。为了使周期图作为功率谱估计满足一致估计的条件，我们必须将周期图进行平滑(或平均)处理。 如果我们需要求周期图的协方差，则有如果以及(k和l为整数)则上式成为上式在以及不是N的整倍数时等于零。因此以的整倍数为频率间距的周期图的值是不相关的，当N增大时，协方差为零的功率谱样本之间的间距减小。因此可以预料到N增加将会使周期图的起伏加快。图6.3说明了周期图的这种特性。图(a)、(b)、(c)分别表示N为16，32和64时的高斯白噪声序列的周期图。图6.3 高斯白噪声序列的周期图由于 故有 虽然本节推导的结果是以假设高斯概率密度为根据的，但其定性结果在一个相当宽广的范围内成立。总之，周期图法不是功率谱的一致估计，而且周期图都在真实功率谱附近随机起伏。这主要是由于时域序列加窗后，频域发生泄漏现象，包括主瓣平滑作用和旁瓣泄漏作用。所以人们作了许多改进，加各种形式的窗，还有分段求估计再求平均的平均周期图法，平滑周期图法等。然而，由于将测得数据以外的数据均看作零（实质相当于乘上了一个矩形窗口），因此频率分辨率低，频谱能量向旁瓣泄漏，且不是功率谱的一致估计。即使进行改进，也不能从根本上改善周期图的性能，只适用于数据较长和频率分辨率要求不高的场合同时，如果增加信号的数据点数，数据量增大，相当于滞后窗加宽，相应的傅立叶变换的主瓣变窄，即功率扩散的频率范围变窄。同时，方差不会降低，随机起伏现象仍然存在。且数据点越多，随机起伏越密集。经典谱估计的优点是计算简单。 6.4 平滑后的周期图作为PSD的估计 上面我们已证明周期图作为功率谱的估计不满足一致估计的条件。因为当，方差，因此必须作一些修正。 本节将表明如果将周期图进行平滑(平均是一种主要的平滑方法)将会使方差减小，得到一致的谱估计。平滑的主要方法有二种。在FFT出现并广泛应用以前主要用窗口处理法进行平滑，即选择适当的窗口函数作为加权函数进行加权平均来加快收敛速度。另一种平滑方法是平均周期图的方法，即先将数据分段，再求各段周期图的平均值。后一种方法又称为Bartlett方法，是当前用得最多的一种平滑方法。Welch又对Bartlett方法进行了改进并提出了用FFT计算的具体方法。 6.4.1 巴特利特(Bartlett)平均周期图的方法 首先让我们来看一下为什么周期图经过某种平均(或平滑)后会使它的方差当时趋于零，达到一致估计的目的。 如果是不相关的随机变量，每一个具有期望值，方差，则可以证明它们的数学平均的期望值等于，数学平均的方差等于，现证明如下： 所以 (6.26)由(6.26)可见，L个平均的方差比每个随机变量的单独方差小L倍。当，可达到一致谱估计的目的。因而降低估计量的方差的一种有效方法是将若干个独立估计值进行平均。把这种方法应用于谱估计应归功于Bartlett。 Bartlett平均周期图的方法是将序列分段求周期图再平均。设将分成L段，每段有M个样本，因而，第i段样本序列可写成第i段的周期图为假定各段的周期图是互相独立的，则谱估计可定义为L段周期图的平均，即(6.27) 于是它的期望值为将式(6.20.0)与式(6.20)代入上式得 (6.28)这里，因此Bartlett估计的期望值是真实谱与三角窗函数的卷积。由于三角窗函数不等于函数，所以Bartlett估计也是有偏估计即，但当时，。 由于我们假定各段周期图是相互独立的，所以可按式(6.26)得到下式： (6.29)由此可见，随着L的增加是下降的，当时，。因此Bartlett估计是一致估计。 M和N如何选择？ 比较式(6.28)的与式(6.20)的可见在二种情况的估计量的期望值都是真值与窗口函数的卷积形式，但后者将前者WB中N改为M，。因而使主瓣的宽度增大。由于主瓣的宽度愈窄，愈接近函数，偏差愈小。式(6.28)中的主瓣宽度大于式(6.20)中的主瓣宽度，因而，而主瓣愈宽显然使分辨率愈差。因此Bias可用来说明谱的分辨率，Bias愈大说明谱分辨率愈差。一个固定的记录长度N，周期图分段的数目L愈大将使方差愈小，但M也愈小，因而使Bias愈大，谱分辨率变得愈差。因此Bartlett方法中Bias或谱分辨率和估计量的方差间是有相互制约关系的。M和N的选择一般是要根据所研究的信号由用户选定。例如，如果我们知道谱有一个窄峰，同时如果分辨出这个峰是重要的，那么我们必须选择M足够大。又从方差的表达式我们可以确定谱估计的可接受的方差所要求的记录长度。 由此可见Bartlett法使谱估计的方差减小是以增加Bias以及降低谱分辨率为代价的。实际上，当N一定时，Bias与Var的相互制约性是谱估计的一个固有特性。 例如 为了说明经平均后的周期图作为功率谱估计的实际效果，设有一零均值高斯分布的随机过程，其功率谱密度为这一功率谱密度是由一零值高斯分布单位方差的噪声序列通过一个其的滤波器形成的。为了简便，设选用的矩形窗函数。其周期图的期望值(用表示)与真值均表示在图6.4中。它说明的周期图可以得到的Bias的情况。图6.5表示分4段与16段二种情况平均后的周期图。显然L=4的方差比L=16的大。将L=16的曲线与图6.4的曲线比较可见在这种估计中误差的大部分起因于Bias而不是方差(因为二种情况均为，Bias相仿，误差也相仿)。M=16，L=2及8的周期图表示在图6.6。图6.6中L不同造成的影响也是明显的。但是这二个曲线的起伏都很大，因此有理由认为误差主要起因于方差。比较与M=16的周期图可见，在6.3中得到的M愈大(即6.3中的N愈大)，将使周期图的起伏愈增快的结论，在这里也同样成立。 比较图6.5与图6.6发现在这个例子中最好的选择是应用L=16，的估计而不是L=8、M=16的估计，即宁可减小方差，牺牲Bias。在实际中，当然功率谱密度的真值是不知道的。但是谱的窗口函数以及关于功率谱密度的某些信息往往是预先知道的。通过改变M和L以及利用预先知道的情况，通常可以很好地进行选择。平均周期图的方法特别适合于应用FFT算法。因此在FFT出现以后这个方法比下面将要讨论的利用窗口函数的处理法用得更多。而在FFT出现以前主要是用窗口函数处理法平滑周期图。 图6.4 与的特性 图6.5 平滑后的周期图（每段取8个数据）图6.6 平均后的周期图(每段取16个数据) 6.4.2 加窗平滑周期图 滞后窗的加窗效应源自无限长序列取样时加的数据窗。周期图的期望值被主瓣平滑的程度和功率向旁瓣泄露的多少，取决于加在信号x上的数据窗w(n)的类型。当w(n)为矩形窗，其傅里叶变换的主瓣比较窄，故周期图被平滑的程度最小，但是由于矩形窗的傅里叶变换的旁瓣比其他窗的旁瓣要高，因而旁瓣泄漏也最严重。如果可以选择合适的窗，比如海明窗，以减弱旁瓣泄漏，使弱窄带信号也显现出来，从而提高频率分辨率，但是旁瓣的下降是以主瓣的变宽为代价的。因此，平滑周期图的另一种常用的方法，下面讨论加窗平滑后周期图的偏差。平滑周期图，实际上就是用一个合适的窗口谱函数与周期图卷积，即(6.30)或(6.31)这里与分别是与的傅氏反变换，并设序列长2M-1。在第一章中我们已讲到是一个实，偶，非负函数，应是一个偶序列，并且满足条件：(6.32)例如前面用到的三角或Bartlett窗函数是满足此条件的，而海宁(Hanning)与汉明(Hamming)窗函数就不满足式(6.32)的条件虽然这后二种窗函数能够提供较好的频率分辨率以及较低的旁瓣，但会产生负的功率谱估计。 如果我们将时域与频域对换，由式(6.30)的卷积形式，可以看成是通过一个具有“单位样本响应”的滤波器产生的。的“频率响应”具有低通特性，因此可以平滑，对中快变化成分有滤除作用，因此我们称为平滑后的周期图。 因式(6.30)可得又按式(6.20)，有因此，有这里现在如果M比N小得多，则，于是式(5.48)可近似写成(6.33)将式(6.33)与式(6.20)比较可见，如将式(6.20)右边的换成，则式(6.20)就全同于式(6.33)。事实上，由于周期图法是将观察到的x(n)的N个以外的值均认为是零，这实际上已把x(n)乘上了一个矩形窗口(即使我们并没有作窗口处理)。式(6.20)中的三角窗函数正是这个矩形窗口造成的，加窗口处理后将式(6.20)改变成式(6.33)。等于将式(6.20)中的窗函数改变成。如果的主瓣宽度大于的主瓣宽度，则可进一步平滑谱估计，减小谱估计的方差，但是同样带来Bias的增加和谱分辨率的变差。至于选用怎样的窗函数才能有较少的Bais的增加与分辨率的降低来换取较多的谱估计的方差的改善以及窗口处理对谱估计方差的影响的具体研究，这里不再讨论。常用窗函数例如：汉宁窗函数：海明窗函数：例题1：为高斯白噪声，用平均周期图法和Bartlett法求信号的功率谱密度估计。N=1024;Ns=256;pxxo=10*log10(abs(fft(xn).2)/(N-1);%周期图无重叠分段平均法pxx1=abs(fft(xn(1:256),Ns).2)/Ns;pxx2=abs(fft(xn(257:512),Ns).2)/Ns;pxx3=abs(fft(xn(513:768),Ns).2)/Ns;pxx4=abs(fft(xn(769:1024),Ns).2)/Ns;pxxb=10*log10(pxx1+pxx2+pxx3+pxx4)/4);结果分析：从频谱看出，与周期图谱估计比较，无重叠分段平均周期图法Bartlett法估计曲线均较为平坦，但是波峰加宽，频率分辨率下降。 6.4.3 Welch法Welch提出对Bartlett法的修正，使其更适合于用FFT进行计算。他主要提出二方面的修正，一方面是在分段时，可使各子序列数据段之间有重迭，这样将会使方差减小(当N与M一定时)。他建议以重迭5070%为宜；另一方面是是选择适当的窗函数，并在周期图计算前直接加进每个子序列。设子序列为xi(n) 长度为L,相邻两个子序列有LD点重叠，于是第i个子序列为：K个子序列总长度为N=L+D(K-1)。若D=L，则N=KL，即为Bartlett法。这样得到的每一段的周期图为(6.34)这里为归一化因子，而Bartlett法每段的周期图为(6.35)这样，Welch周期图表示式为：或 (6.36)这样加窗函数的优点是无论什么样的窗函数均可使谱估计非负。 例题2: 对例题1的信号分别分析分段加窗平滑周期图和重叠分段加窗平滑周期图。pxxo=10*log10(abs(fft(xn).2)/(N-1);%周期图无重叠分段加窗平均周期图w=hanning(256);pxx1=abs(fft(w.*xn(1:256),Ns).2)/norm(w)2;pxx2=abs(fft(w.*xn(257:512),Ns).2)/norm(w)2;pxx3=abs(fft(w.*xn(513:768),Ns).2)/norm(w)2;pxx4=abs(fft(w.*xn(769:1024),Ns).2)/norm(w)2;px