CN1933434A - 频谱估计的方法和装置 - Google Patents

Abstract

估计100、200、300、430来自所采样的数据的数据流的样本集合的频谱，直到实现110、220、320、340目标信噪比为止。

Description

频谱估计的方法和装置

技术领域

本发明涉及频谱估计的方法和装置。

背景技术

数据收集网络通常在第一位置处对数据进行采样，然后将所采样的数据发送到另一个位置以便处理和分析。在某些数据收集网络中，可能会以时间函数的方式随机地或无规则地对数据进行采样。具体而言，个体样本之间的时间间隔本质上可能会以时间函数的方式而随机地变化。

这种数据收集网络的示例包括但不限于具有随机采样的宽带测试系统和时间同步的低功率传感器网络。当伴随着样本的精确时间戳(例如时间同步)时，具有随机采样的宽带测试系统促进了利用远低于信号的传统奈奎斯特采样率的平均采样率来进行宽带信号表征。在另一种宽带信号测试情形中，特定测试常要求有限频谱范围上的精确数据(例如射频设备的单音调或双音调测试)。在这种情形下，随机化的数据采样可能会使得执行测试所需的数据总量最小化。对于低功率联网传感器，每个传感器的功率消耗通常与传感器的采样率直接相关。在许多情形下，通过采用随机化的采样来降低数据速率促进了低功率操作。此外，网络所施加的约束(例如网络协议和相关联的定时)常常会对采样间隔产生实际限制，从而导致间距不均匀或不规则的样本。这里通过引用结合进来的Barford的美国专利No.6,735,539B2教导了这样一种系统，它使用了具有间距不均匀的具有时间戳的样本的联网的传感器。

已知多种用于从间距不均匀的、不规则的和/或随机间距的带时间戳的数据处理频谱估计的技术。在已知的技术中，有自回归滑动平均(ARMA)方法、离散傅立叶变换(DFT)方法和矩阵因子分解(MF)方法。ARMA方法将带时间戳的输入数据集合拟合到ARMA模型。ARMA方法要求过高的计算时间，从而使得这种方法对于许多实时应用都不合适。DFT方法一般要求极大的输入数据集合以及相对较小的样本间时间间距，以使DFT的近似傅立叶积分的求和的误差最小化。MF方法包括但不限于尝试将最小二乘频谱拟合到所采样的数据的方法。对于这种MF方法，不存在确定足够数目的样本的实际可行的手段，从而导致需要过长的测量和计算时间，尤其是在考虑实时操作(例如在制造环境中)时更是如此。

因此，需要有一种提供间距规则的和间距随机或不规则的数据样本之一或二者的频谱估计的方法，它在诸如制造和测试这样的环境中提供相对精确的频谱估计。这种方法将会解决频谱估计领域长时间以来存在的需求。

发明内容

在本发明的某些实施例中，提供了对所采样的数据的数据流D进行频谱估计的方法。该频谱估计方法包括从数据流D中选择第一样本集合。该方法还包括选择其他样本集合以逐步与所述第一样本集合相组合以便创建组合样本集合。其他样本集合被逐步组合，直到所述组合样本集合的所估计的频谱的信噪比的量度指示实现了目标信噪比为止。

在本发明的另一个实施例中，提供了一种频谱估计器。该频谱估计包括计算机处理器、存储器和存储在存储器中并被处理器执行的计算机程序。计算机程序包括这样的指令，这种指令在被处理器执行时，实现从带时间戳的样本的数据流中选择第一样本集合，并且还实现随后选择带时间戳的样本的其他集合，以逐步与所述第一集合相组合以便创建组合样本集合。样本被逐步组合，直到所述组合样本集合的所估计的频谱的信噪比的量度指示实现了目标信噪比为止。

在本发明的另一个实施例中，提供了采用频谱估计的测试系统。该测试系统包括连接到被测设备的输入的激励源和连接到所述被测设备的输出的采样器。该系统还包括连接到采样器的输出的频谱估计器。采样器将来自被测设备的所采样的数据提供给频谱估计器。频谱估计器实现从所采样的数据中选择带时间戳的数据以产生样本集合，实现计算所述样本集合的所估计的频谱的信噪比的量度，实现逐步选择附加的带时间戳样本以与所述样本集合相组合，并且实现重新计算组合样本集合的所估计的频谱的信噪比的量度。逐步选择和重新计算量度被重复，直到量度小于预定限度λ为止。当量度小于预定限度λ时，实现了目标信噪比。

本发明的某些实施例除了上述特征外还具有其他特征。本发明的这些和其他特征在下文中参考附图详细描述。

附图说明

在结合附图参考以下详细描述的情况下，可以更容易理解本发明的实施例的各种特征，附图中类似的标号指代类似的结构元件，其中：

图1示出根据本发明实施例的频谱估计方法的流程图；

图2示出根据本发明实施例的构造最小二乘(L-S)矩阵A的流程图；

图3示出根据本发明实施例的用于所采样的数据的数据流D的频谱估计方法的流程图；

图4示出根据本发明实施例的在感兴趣的频率集合下估计所采样的数据的频谱的方法的流程图；

图5示出根据本发明实施例的采用频谱估计的测试系统的框图；

图6示出根据本发明实施例的频谱估计器的框图。

具体实施方式

本发明的实施例促进了利用来自以时间函数的方式采样的信号或过程的数据来估计信号或过程的频谱。被时间采样的并且被估计频谱的信号或过程本质上可以是任何信号或过程，包括但不限于被联网的传感器远程监视的信号或过程，和与测试系统相关联的信号。例如，信号可以是来自被测设备(DUT)的信号，其中测试或测量系统对信号采样。在另一个示例中，传感器可以以时间函数的方式对物理过程(例如温度、辐射等)进行采样，并且将所采样的数据传递到处理器，以便进行频谱估计。本发明的某些实施例实质上适用于对所采样的数据执行频谱估计的任何情形。

根据本发明的各种实施例，频谱估计可以采用间距规则的、时间采样的数据，或者间距不规则的(即随机的)、时间采样的数据。具体而言，间距规则的时间采样数据包括以恒定时间间隔采样的数据点集合。例如，利用周期性采样时间间隔采样的数据可以生成间距规则的时间采样数据。间距不规则或随机的时间采样数据包括这样的数据点集合，其中相邻数据点对之间的采样时间间隔或时间间距以时间函数的方式而不同或变化。例如，间距规则或随机的时间采样数据是用作为随机时间函数的采样时间间隔来生成的。在另一个示例中，间距不规则的时间采样数据是在采样时间间隔本质上具有规则的间距时生成的，只不过数据中存在缝隙。数据中的缝隙例如可能是由一个或多个连续数据点的传输期间的丢失产生的。

在某些实施例中，所估计的频谱以向量形式提供。具体而言，在这种实施例中，所估计的频谱是作为幅度和相位值的向量提供的。在其他实施例中，所估计的频谱是作为包括实值和虚值的向量提供的。此外，某些实施例在任意的预定频率下估计频谱。

根据本发明的某些实施例的频谱估计采用了最小二乘方法，该方法将“最小二乘”频谱拟合到可用的所采样的数据。在某些实施例中，由频谱估计产生的估计的频谱使用接近最小数目的样本。具体而言，根据这些实施例，接近最小数目的样本是利用仅略多于求解最小二乘拟合以估计频谱所需的多个浮点操作来识别的。在某些实施例中，可以通过控制频谱估计操作来进行频谱计算时间和频谱精度之间的折衷。

一般来说，数据的每个数据点与指示数据点的采样时间的时间戳相关联。时间戳可以是由采样器生成并与数据点一起被发送到处理器的实际或显式时间戳。例如，采样器可以记录一对值(t_n，y_n)，其中“t_n”代表实际采样时间，“y_n”表示实际样本值。或者，时间戳可以是虚拟时间戳，其中只有样本值被采样器记录，而采样时间是隐式地提供的。例如，采样器可以采用具有先验已知的采样时间序列的采样时钟或采样定时器。先验已知(即隐式或虚拟)采样时间彼此之间可以有规则的或不规则的(例如随机的)的间距。这样，虚拟时间戳通常与其中采样时间是从采样器生成样本的方式推断出来的数据相关联。

此外，时间戳可以是相对时间戳(例如指示从上次样本开始经过的时间)或者绝对时间戳(即基于固定参考或开始时间)。这样，实际时间戳可以是相对时间戳或绝对时间戳。类似地，虚拟时间戳可以是相对的或绝对的。例如，可以采用绝对实际时间戳、相对实际时间戳、绝对虚拟时间戳和相对虚拟时间戳中的一个或多个。为了论述简便起见，并且不是以限制性的方式，以下采用绝对实际(显式)时间戳。本领域的技术人员可以很容易地采用上述其他时间戳中的任何一种，而仍处于本发明实施例的范围内。

这里所采用的“随机采样的”数据是其中采样时间没有以时间函数的方式具有规则间距的所采样的数据。例如，由利用在平均采样间隔周围随机变化的采样时钟对模拟信号采样的测量设备(例如传感器)生成的所采样的数据就是随机采样的数据。在另一个示例中，随机采样的数据代表这样的数据样本：这些数据样本本质上是规则采样的数据，但是在数据样本序列中却具有偶发的缝隙或丢失的样本。在随机采样的数据的另一个示例中，采样时间间隔以时间函数的方式(例如时间的正弦函数)变化。为了以下论述简便起见，将假定随机采样的数据。但是，本领域的技术人员可以很容易地采用符合奈奎斯特采样标准的间距规则的数据，而仍处于本发明实施例的范围内。

例如，带时间戳的数据可以在相对于数据的奈奎斯特采样率较慢的平均采样率下被取得或采样。具体而言，平均采样频率可能低于要估计的频谱中的任何感兴趣的频率。此外，数据中相对长的缝隙是可以容忍的。这样，可以采用具有相对于被测量的信号的时间比例有较长的恢复时间的测量设备来取得样本。允许较长的恢复时间和/或相对较低的采样频率通过使测量设备的实时方面集中在生成时间戳上而远离数据递送，从而促进了测量设备设计。

在某些实施例中，带时间戳的数据由数据流表示。这里所使用的“数据流”是指在一段时间中到达处理位置的固定大小的数据点集合(例如固定长度的带时间戳数据点列表)或带时间戳的数据点的序列。具体而言，当采用带时间戳数据点序列时，某些实施例可以简单地在继续处理数据流之前暂停或等待下一个或附加的新的带时间戳数据点到达。从而，某些实施例可以用来生成带时间戳数据的“完整”集合的频谱估计，其中完整集合代表所有可获得或将会获得的数据。其他实施例可以生成初步估计，该初步估计随后在数据流中的新数据到达并被处理时被更新。

一般而言，数据流D可以被表示成值对的集合(t_n，y_n)，其中每对包括时间戳t_n值和相应的数据点y_n值，其中n是每对的索引，并且术语“索引”是指该对在数据流内的位置。具体而言，在某些实施例中，数据流D是值对的集合(t_n，y_n)，其中该集合具有基数或大小|T|(即，n＝1，2，…|T|)。因此，数据流D可以由方程(1)定义：

         D＝{(t₁，y₁)，(t₂，y₂)，…，(t_T，y_|T|)}        (1)

由于大小|T|是可变的，因此数据流D的这种定义既适应固定大小集合也适应顺序到达的集合，而没有含糊性。

本发明的某些实施例在以m为索引的任意预定感兴趣频率f_m下估计频谱。在某些实施例中，大小为|F|的感兴趣的频率的有序集合F(即，m＝1，2，…|F|)由方程(2)定义：

                    F＝{f₁，f₂，…，f_|F|}               (2)

其中0＜f₁＜f₂＜…＜f_|F|。在某些实施例中，数据流D的频谱是在f₀＝0Hz下估计的。

由本发明的某些实施例产生的估计的频谱是一个向量S，该向量S包括在感兴趣的频率的有序集合F中的每个频率下的多对频谱系数，例如实频谱幅度和虚频谱幅度(c_m，s_m)，或者相应的幅度和相位(A_m，φ_m)。在某些实施例中，所估计的频谱S由方程(3)给出：

            S＝{(c₀，s₀)，(c₁，s₁)，…，(c_|F|，s_|F|)}

                          或                            (3)

            S＝{(A₀，φ₀)，(A₁，φ₁)，…，(A_|F|，φ_|F|)}

其中m＝1，2，…，|F|。

频谱估计确定可用数据到频谱S的最小二乘拟合。一般来说，所确定的最小二乘拟合本质上代表由方程组所例示的超定问题的近似解，其中一般有比解答中的未知数更多的方程(即，超定最小二乘(L-S)拟合问题)。在估计的频谱S的情况下，最小二乘拟合找出最佳地近似真实频谱的频谱系数(例如(c_m，s_m)或者等同的(A_m，φ_m))。在矩阵表示法中，超定L-S拟合问题由Ax＝b表示，其中x和b是向量，A是矩阵。一般来说，求解任何超定问题的目标是找出使得Ax合理地“接近”b的解答x。就最小二乘拟合来说，目标是找出使得差Ax-b在“最小二乘”意义上最小(即，使‖Ax-b‖₂最小，其中‖·‖2是矩阵2范数)的解答向量x。

一般来说，采用各种分解或因子分解，其中包括但不限于QR因子分解和奇异值分解(SVD)，来执行最小二乘(L-S)拟合。QR因子分解确定一对因子矩阵Q和R，其积等于矩阵A。一般来说，矩阵A是实数的m乘n矩阵(即

)，其中m和n是整数，并且m≥n。因子矩阵R是m乘z上三角矩阵(即， )，它定义为这样一个矩阵，其中主对角线下的所有非对角线元素都等于零，而因子矩阵Q是正交m乘m矩阵(即 )。从而，矩阵A的QR因子分解可以由方程(4)给出：

A＝QR       其中

其中，因子矩阵Q和其转置的积等于单位矩阵(即，Q^TQ＝I)。

一般来说，矩阵A的QR因子分解可以以多种不同方式来实现，其中包括但不限于Householder变换、Givens变换、分块Householder变换、以及快速Givens变换。以上列出的每种变换以及其他QR因子分解的细节可以在本领域的技术人员已知的关于线性代数和矩阵计算的多本教科书中找到，其中包括但不限于Golub和Van Loan所著的Matrix Computations，2^ndEd.，The Johns Hopkins University Press，Baltimore，MD，1989，这里通过引用将其完全结合进来。在实践中，QR因子分解一般可以以与计算机程序和相关数值计算环境相关联的函数或子例程的形式获得，所述计算机程序和相关数值计算环境包括但不限于注册给Mathworks，Inc.，Natack，MA并由其发表的Matlab，以及由http://www.netlib.org分发和维护的LAPACK。本领域的技术人员将会很容易认识到，可以用其他的求解L-S方程的分解和方法替代QR分解，而不会脱离本发明的各种实施例的范围。

图1示出了根据本发明实施例的频谱估计方法100的流程图。频谱估计方法100包括构造110最小二乘(L-S)问题的最小二乘(L-S)矩阵A，其将所采样的数据拟合到所估计的频谱S。具体而言，L-S矩阵A包括所采样的数据的傅立叶变换的基集合。L-S拟合问题的解答向量x包括傅立叶变换的系数，其中系数本质上代表所估计的频谱S。虽然这里为了论述是就矩阵来表达的，但是本领域的技术人员可以在不进行过度的实验的情况下很容易地将下面的内容延伸到其他等同表示方式。

在某些实施例中，L-S矩阵A在第一列中包括比例因子，即[1/]。L-S矩阵A还包括第一列之后的第一子矩阵和第二子矩阵。第一子矩阵包括作为任意预定频率F和数据流D的所选样本的时间戳t_n的正弦函数的元素。第二子矩阵包括作为任意预定频率F和数据流D的所选样本的时间戳t_n的余弦函数的元素。在某些实施例中，L-S矩阵A由方程(5)给出：

$A=\left[\begin{array}{ccc}\frac{1}{\sqrt{2}}& \sin \left(2{\mathrm{\πt}}_{1}f\right)& \cos \left(2\mathrm{\π}{t}_{1}f\right)\\ \frac{1}{\sqrt{2}}& \sin \left(2\mathrm{\π}{t}_{2}f\right)& \cos \left(2\mathrm{\π}{t}_{2}f\right)\\ \frac{1}{\sqrt{2}}& \sin \left(2\mathrm{\π}{t}_{3}f\right)& \cos \left(2\mathrm{\π}{t}_{3}f\right)\\ \·& \·& \·\\ \·& \·& \·\\ \·& \·& \·\\ \frac{1}{\sqrt{2}}& \sin \left(2\mathrm{\π}{t}_{\left|T\right|}f\right)& \cos \left(2\mathrm{\π}{t}_{\left|T\right|}f\right)\end{array}\right]...\left(5\right)$

其中行是以时间戳t_n为索引的，索引n的范围是1到|T|(即，n＝1，2，…|T|)。此外，函数sin(2πt_nf)和cos(2πt_nf)分别是以频率的有序集合F中的频率f_m为索引的行向量，并且由方程(6)和(7)给出：

sin(2πt_nf)＝[sin(2πt_nf₁)，sin(2πt_nf₂)，…，sin(2πt_nf_|F|)]    (6)

cos(2πt_nf)＝[cos(2πt_nf₁)，cos(2πt_nf₂)，…，cos(2πt_nf_|F|)]    (7)

其中方程(5)、(6)和(7)中的“f”代表包含频率集合F中以m为索引的频率的向量，m的范围是1到|F|(即，m＝1，2，…|F|)。为了进行论述，方程(6)和(7)中的索引“n＝1”代表第一所选样本，索引“n＝2”代表第二所选样本，以此类推，直到最后的所选样本由索引“n＝|T|”表示。在某些实施例中，第一所选样本不需要是数据流D中的第一样本，第二所选样本不需要是数据流D中的第二样本，以此类推。在某些实施例中，最后的所选样本是具有小于或等于数据流D的最后一个样本n＝|T|(即(t_|T|，y_|T|))的时间戳的样本。

本领域的技术人员将会容易认识到，可能存在用于确定估计的频谱S的其他形式的矩阵A，从而可以用来代替以上所描述的那些，而不会脱离本发明的范围。例如，另一种已知的傅立叶变换的形式采用指数基集合，而不是等同的正弦/余弦基集合。这种傅立叶变换将会产生包括e^j2πtf形式的项的矩阵A。矩阵A的这种替换形式很明显是在本发明的范围内的。

一般来说，通过从数据流D中选择一定数目的样本来构造110L-S矩阵A，其中该数目足以实现由L-S拟合产生的所估计频谱S的目标信噪比(SNR)。具体而言，从数据流D中选择样本集合。然后，构造110L-S矩阵A，并且通过求解L-S拟合问题来从中估计频谱S。然后计算所估计的频谱S的SNR，并且如果SNR小于目标SNR，则从数据流D中选择出附加的样本，并在新构造的L-S矩阵A中将其添加到原始样本。再次确定所估计的频谱S并计算SNR。重复选择和计算的过程，直到实现目标SNR为止。

在某些实施例中，选择样本使得L-S矩阵A非奇异，并且L-S拟合问题的条件数足够大。在这种实施例中，可以采用与L-S矩阵A的条件数成比例的量度。例如，可以采用L-S矩阵A的QR分解来计算或生成一对L-S因子矩阵Q和R。然后，可以采用基于量度的因子矩阵R来确立相对于实现目标SNR的目标的足量的所选样本。具体而言，如果因子矩阵R满足以下条件，则L-S矩阵A和因子短阵R就都是非奇异的：因子矩阵R的主对角线元素的绝对值的最小值元素严格大于零。此外，如果因子矩阵R的主对角线元素的绝对值的最大值元素与因子矩阵R的主对角线元素的绝对值的最小值元素之比小于任意预定限度λ，那么因子矩阵R就是非奇异的。结果，所估计的频谱S将会具有足够大的信噪比(SNR)。换言之，所估计的频谱的SNR将会实现目标SNR。就Matlab型矩阵函数max(·)、min(·)、abs(·)和diag(·)来说，如果因子矩阵R满足由方程(8)给出的条件，则将会满足目标SNR：

min(abs(diag(R)))＞0^max(abs(diag(R)))/min(abs(diag(R)))＜λ    (8)其中min(·)是返回向量自变量的最小元素的按元素最小函数，max(·)是返回向量自变量的最大元素的按元素最大函数、abs(·)是返回矩阵自变量的逐元素绝对值(即幅值)的向量的按元素绝对值函数，diag(·)是返回包含矩阵自变量的主对角线元素的向量的对角线函数，符号^表示有序“与”函数，其中左侧自变量先于右侧自变量被计算。

对于给定情形，可以确定出预定限度λ而无需过度实验。在某些实施例中，约等于“1.1”的预定限度λ(例如，λ＝1.1)本质上保证了给定矩阵R的条件数约为1(即， $\mathrm{cond}\left(R\right)\≅1$ )。当矩阵R的条件数约为1时，所估计的频谱S的SNR近似为输入数据流D的SNR。在其他实施例中，尤其是在要求更快的估计(即更快的估计速度)的情况下，可以采用预定限度λ的更大的值。一般而言，可以采用从约为1到1000量级的预定限度λ(即，1≤λ≤1×10³)，以提供所估计的频谱S的SNR和估计速度之间的折衷。在某些实施例中，任意限度λ的范围从大约1到大约100。

在其他实施例中，实现与方程(7)所给出的几乎相同的结果的关于矩阵R的更简单的条件由方程(9)给出：

min(abs(diag(R)))＞0^1/min(abs(diag(R)))＜λ            (9)

具体而言，如果R是严格非奇异的，那么如果矩阵R的主对角线元素的绝对值的最小值元素的倒数小于任意限度λ(例如小于“1.1”)，就能实现目标SNR。

除了确定要采用的数据流D的一定数目的样本(例如(t_n，y_n))之外，在某些实施例中，还确定样本之间的最小时间步长或时间间距，以确保样本充分地跨越或以其他方式代表所估计的频谱S中的感兴趣的特性。具体而言，在某些实施例中，选择样本之间的最小时间步长Δt_min，以使得从不等于零的最低感兴趣频率f₁的至少两个周期中取得感兴趣的频率集合F中的频率数目的两倍加1(即2|F|+1)个样本。由方程(10)给出的最小时间步长Δt_min一般确保了从最低感兴趣频率中取得足够的样本来执行频谱估计：

             Δt_min＝2/f₁/(2|F|+1)                     (10)

在某些实施例中，选择样本集合或样本群组，而不是个体的样本。在实践中，可以选择毗邻或相邻样本的群组以通过平均来提高所估计的频谱S的SNR，以使其超过和高于通过选择个体样本能提供的。换言之，不是选择第一样本(t_a，y_a)，然后选择第二样本(t_b，y_b)，而是可以选择数据流D中的相邻样本的群组(对、三元组、四元组，等等)。例如，选择第一对{(t_a，y_a)，(t_a+1，y_a+1)}，然后是第二对{(t_b，y_b)，(t_b+1，y_b+1)}，以此类推，其中样本时间t_n和t_n+1(其中n＝a或b)表示数据流D中的相邻样本。在另一个示例中，选择第一三元组{(t_a，y_a)，(t_a+1，y_a+1)，(t_a+2，y_a+2)}，然后是第二三元组{(t_b，y_b)，(t_b+1，y_b+1)，(t_b+2，y_b+2)}，以此类推。在某些实施例中，第一对的“t_a”和第二对的“t_b”之间的时间步长或间距至少是最小时间步长Δt_min。以下，术语“oversample”将被用来表示相邻样本的集合的大小。例如，等于2的oversample值(即oversample＝2)表示2个或一对相邻样本的集合。

图2示出了根据本发明实施例的构造110最小二乘(L-S)矩阵A的流程图。如图2所示，构造110包括从数据流D中选择初始样本集合。选择112初始集合以使得所选样本的数目至少是感兴趣的频率F的数目的两倍加1(即，至少是2|F|+1)。每个所选样本之间的时间间距大于或等于最小时间步长Δt_min。当相邻样本的群组被选择112时(即oversample＞1)，所选群组的数目至少是2|F|+1，群组之间的时间间距(例如每个群组的第一成员)大于或等于最小时间步长Δt_min。

例如，如果被称为“samplesused”的集合代表与所选样本相对应的数据流D的索引列表，并且如果“I”是另一个索引列表，并且“oversample”是每个群组中的样本数目，则在某些实施例中，选择112包括(a)设置列表I等于数据流D的oversample个索引的群组，从而使得群组中的第一样本的时间戳t_i比集合samplesused中的最后一个索引所标识的样本的时间戳t_m至少大最小时间步长Δt_min；(b)将列表I添加到集合samplesused末尾(即，集合samplesused＝samplesusedu∪I，其中∪是并集操作符)，其中(a)和(b)被重复直到(samplesused≥2|F|+1)为止。

构造110还包括利用感兴趣的频率M和所选择112的样本的时间戳t_n来初始化114 L-S拟合的L-S矩阵A。在某些实施例中，初始化114采用以上由方程(5)定义的L-S矩阵A。

构造110还包括对经初始化114的L-S矩阵A执行116 QP，因子分解，以生成因子矩阵Q和R，其中因子矩阵Q是正交矩阵，因子矩阵R是上三角矩阵。例如，执行116 QR因子分解可以通过使用例如但不限于以下函数来完成：基于计算机程序的矩阵计算环境Matlab^的“qr(A)”函数或矩阵函数库LAPACK的函数“DGEQRF”。

构造110还包括利用一个或多个附加的样本来更新118 L-S拟合，直到频谱估计方法100所产生的频谱估计实现目标信噪比为止。在某些实施例中，更新118包括将一个或多个附加的样本添加到L-S矩阵A，然后根据更新后的L-S矩阵A重新计算因子矩阵Q和R。在其他实施例中，更新118包括添加一个或多个附加的样本以更新因子矩阵Q和R，而不先显式地更新L-S矩阵A。在其他实施例中，L-S矩阵A和因子矩阵Q和R都被更新118。

例如，由于任何经更新118的L-S矩阵A包括先前的L-S矩阵

加上所添加的样本，因此经更新118的因子矩阵Q和R可以直接从添加的样本和先前的因子矩阵 和 计算出来。具体而言，更新118包括将与添加的样本相对应的新的或附加的行添加到先前的L-S矩阵

和先前的因子矩阵

和 中的任何一个或两者。这样，可以直接从先前的因子矩阵

和

的行计算出附加的行，而无需显式地更新118 L-S矩阵A并对其执行116 QR因子分解。直接更新118因子矩阵Q和R的方法由上文引用的Golub等所著的Matrix Computations，pp.596-597给出，这里通过引用将其结合进来。此外，更新118因子矩阵Q和R本质上消除了对在执行116 QP分解之后将L-S矩阵A显式地存储在计算机存储器中的任何需求。

在更新118期间，附加的样本被从数据流D中选择出，并被用于更新，直到达到与目标信噪比相关的条件为止。在某些实施例中，该条件由上述方程(8)给出。在某些实施例中，采用由上述方程(9)给出的条件。在其他实施例中，可以采用因子矩阵R和L-S矩阵A中的一个或两个的条件数。

例如，在某些实施例中，选择预定限度λ(例如 $\mathrm{\λ}\≅1.1$ )，并且从开始于之前所使用的最后一个样本的数据流D中选择样本群组(t_n，y_n)。所选择的样本群组在时间上的间距至少是最小时间步长Δt_min，如上文针对选择112所述。此外，这些群组包含的个体样本数目等于oversample，如上文针对选择112所定义和描述。所选样本(t_n，y_n)的时间戳t_n被用来创建Q和R矩阵中的新的行，从而更新118矩阵。选择群组并创建新的行，直到例如符合方程(7)的条件为止。在符合条件的情况下，就实现了所估计的频谱S的目标信噪比。

所使用的样本的索引的集合(例如samplesused)也被更新，以包括添加的所选样本。例如，更新samplesued的索引集合以包括添加的所选样本可以包括(a)设置列表I等于数据流D的oversample个索引的群组，从而使得群组中的第一样本的时间戳t_i比集合samplesused中的最后一个索引所标识的样本的时间戳t_m至少大最小时间步长Δt_min；(b)将I添加到集合samplesused末尾(即，集合samplesused＝samplesused∪I)，其中(a)和(b)被重复直到符合方程(7)的条件为止。

返回参考图1，频谱估计方法100还包括求解120L-S拟合问题，以获得代表傅立叶变换的系数的解答向量x。求解120包括设置向量b等于样本值y_n，其中n是在所使用的样本上索引的，其中包括在更新118期间添加的那些(例如b＝[y_n]，n＝1，…，samplesused)。求解120还包括根据方程(11)利用向量b和因子矩阵Q和R来解出解答向量x：

                     Rx＝Q^Tb                          (11)

其中“Q^T”是因子矩阵Q的转置。本领域的技术人员对于求解由方程(11)例示的矩阵方程的方法是很熟悉的。

频谱估计方法100还包括从解答向量x中提取130傅立叶系数。提取130傅立叶系数产生所估计的频谱S。在某些实施例中，提取130包括(a)设置第一实系数c₀等于解答向量x的第一元素x₁；以及(b)设置第一虚系数s₀等于“0”(即，C₀＝x₁，S₀＝0)。对于从f₁开始的感兴趣的频率的集合F的每个频率，设置(a)包括设置每个连续的实系数ci等于解答向量x的索引为i加上感兴趣的频率的集合F的大小加1的个体元素。此外，对于从f₁开始的感兴趣的频率的集合F的每个频率，设置(b)包括设置每个连续的虚系数s_i等于解答向量x的索引为i加1的个体元素。换言之，对于i∈{1，…，|F|}，设置c_i＝x_i+|F|+1，其中x＝[x₁，x₂，…，x_2|F|+1]。

在某些实施例中，提取130还包括根据实系数和虚系数(c_m，s_m)来计算所估计的频谱S的幅度和相位(A_m，φ_m)。对于从零到感兴趣的频率的集合F的大小的i，第i个幅度A_i是以第i个实系数c_i和第i个虚系数s_i的平方之和的平方根的形式计算的。对于从零到感兴趣的频率的集合F的大小的i，第i个相位φ_i是以第i个虚系数s_i和第i个实系数c_i之比的反正切的形式计算的。换言之，对于i∈{0，…，|F|}，计算幅度和相位包括根据以下方程分别设置第i个幅度A_i和第i个相位φ_i：

$A_i=\sqrt{c_i^2+s_i^2},{\mathrm{\φ}}_i={\tan }^{-1}(s_i/c_i).$

在某些实施例中，在更新118的每次反复之后，因子矩阵Q被存储在存储器中。在其他实施例中，因子矩阵Q不被显式地存储。相反，因子矩阵Q最初被存储，如上文引用的LAPACK的DGEQRF例程的描述中所述。虽然对于本领域的技术人员来说是已知的，但在这里通过引用结合进来的Anderson等人的Lapack Users Guide，3^rd ed.，Society for IndustrialApplied Mathematics Press，2000的题为“QR Factorization”的章节中，提供了对DGEQRF例程以及由此来存储因子矩阵Q的完整描述。然后，在更新118期间，存储描述因子矩阵Q的更新118的一个或多个Givens旋转。从而，取代显式存储和更新118的因子矩阵Q，最初始存储的因子矩阵Q表示和Given旋转集合就代表了经更新118的因子矩阵Q。换言之，不是显式地存储经更新的因子矩阵Q，而是存储矩阵的描述，作为将因子矩阵Q或因子矩阵的转置Q^T乘以向量所需的一组步骤。以这种方式存储因子矩阵Q可以节省存储器和计算时间，因为最终被频谱估计方法100使用的是因子矩阵Q和向量(即b向量)的乘积，而不是因子矩阵Q本身。

图3示出了根据本发明实施例的用于所采样的数据的数据流D的频谱估计方法200。方法200包括从数据流D中选择210第一样本集合。方法200还包括从数据流D中选择220其他样本集合，以相继与第一集合组合，从而创建组合的样本集合。选择其他集合以便组合，直到组合样本集合的所估计频谱的信噪比的量度指示实现了目标信噪比为止。在某些实施例中，选择210与上述方法100中的构造100的选择112初始样本集合本质上类似。此外，在某些实施例中，选择220和方法200采用也是上述方法100的构造110的更新118的步骤的某些方面。

图4示出了根据本发明实施例的在感兴趣的频率集合下估计所采样的数据的频谱的方法300。方法300包括从所采样的数据中选择310带时间戳的样本，以产生样本集合。方法300还包括计算320样本集合的所估计频谱的信噪比的量度。方法300还包括逐步选择330附加的带时间戳样本，以与样本集合相组合，并且重新计算340组合样本集合的所估计频谱的信噪比的量度。重复逐步的选择和重新计算，直到量度小于预定限度λ为止。当信噪比的量度小于预定限度λ时，就实现了目标信噪比。

在某些实施例中，选择310与上述方法100的构造110的选择112初始样本集合本质上类似。在某些实施例中，计算320和方法300采用也是上述方法100的构造110的初始化114和执行116的某些方面。此外，在某些实施例中，逐步选择330和重新计算340与上述方法100的构造110的执行116和更新118的某些方面本质上类似。

图5示出了根据本发明实施例的采用不规则采样的数据的频谱估计的测试系统400的框图。测试系统400对由被测设备(DUT)产生的输出信号进行采样，并且生成输出信号的所估计频谱S。测试系统400以不规则或随机的间隔来对输出信号进行测量和采样。样本是就带时间戳的样本幅度(例如(t_n，y_n))来表达的。在某些实施例中，所估计的频谱S由包括在预定的感兴趣的频率F下频谱S的幅度和相位的向量所代表。在其他实施例中，所估计的频谱S由包括在预定的感兴趣的频率F下的频谱S的实系数和虚系数的向量所代表。具体而言，所估计的频谱S可以由以上方程(3)来描述。

测试系统400包括激励源410。激励源生成被施加到被测设备(DUT)402的激励信号。在某些实施例中，激励源410是生成本质上周期性的激励信号的周期性激励源。在其他实施例中，激励源410是通用激励源(例如任意波形生成器)，其生成非周期性激励信号和本质上周期性的激励信号中的一种或两种。在其他实施例中，激励源生成空信号。例如，当测试结合内部信号生成器(例如时钟生成器，函数生成器等)的DUT 402时，可以采用空信号。这种DUT 402不要求用于测试目的的激励信号。在这种实施例中，可以从测试系统400中省略激励源410，但仍处于本发明的范围内。

测试系统400还包括采样器420。在某些实施例中，采样器420是模数转换器(ADC)。在其他实施例中，基本上任何能够测量物理量并将测得量递送到计算装置的仪器、传感器或其他装置都可以被用作采样器420。示例包括但不限于由具有预定触发时间(例如警报)触发的示波器，可以在网络中工作或者“可联网”的智能传感器(例如参见IEEE1451标准系列)。采样器420对由DUT 402生成的输出信号进行采样，并且将输出信号的样本转换成带时间戳的样本幅度(例如(t_n，y_n))。在某些实施例中，输出信号可以由测量物理过程(例如温度)的传感器402而不是DUT 402生成。

在某些实施例中，采样器420产生包括带时间戳的样本幅度或样本数据点的序列的数据流D。例如，采样器420所产生的数据流D可以由以上针对频谱估计方法100所描述的方程(1)来代表。此外，虽然在这里被描述为显式时间戳，但是与数据流D的样本幅度相关联的时间戳可以是一个或多个显式或隐式时间戳，如以上针对频谱估计方法100所描述。一般来说，采样器420以不规则的间隔采样。但是，采样器420可以以不规则的、随机的和规则的间隔之中的一种或多种间隔来采样。当以规则间隔采样时，不规则样本集合例如可能产生自在从采样器420发送期间的集合中的一个或多个样本的丢失。

测试系统400还包括频谱估计器430。频谱估计器430接收数据流D，并且生成所估计的频谱S。所估计的频谱S是针对任意的预定感兴趣频率的集合F生成的。例如，感兴趣的频率的集合F可以是由以上针对频谱估计方法100描述的方程(2)所定义的有序集合。

图6示出根据本发明实施例的频谱估计器430的框图。在某些实施例中，频谱估计器430确定可用数据到所估计的频谱S的最小二乘拟合。在某些实施例中，频谱估计器430包括计算机处理器432、存储器434以及存储在存储器434中并被计算机处理器432执行的计算机程序436。计算机处理器432可以是通用计算机处理器或者专用处理器，例如但不限于信号处理器。存储器例如可以是随机访问存储器(RAM)、只读存储器(ROM)、闪存或盘存储器中的一种或多种，所述盘存储器包括但不限于硬盘、紧致盘(CD)、软盘和数字视频盘(DVD)。在某些实施例中，计算机程序436可以实现或分类为频谱估计器430的软件或固件。例如，计算机程序436可以用高级计算机编程语言或编程环境来实现，所述高级计算机编程语言或编程环境包括但不限于C/C++、FORTRAN或Matlab^。在其他实施例中，计算机程序436可以实现为计算机逻辑或硬件。例如，计算机程序436本质上可以被硬连线为分立逻辑或专用集成电路(ASIC)。在这种实施例中，计算机逻辑实现方式可以代表计算机处理器432和存储器434以及计算机程序436中的一种或多种。

不论特定实现如何，计算机程序436都包括在被处理器432执行时实现频谱估计的指令。在某些实施例中，频谱估计以最小二乘的方式将来自数据流D的数据拟合到所估计的频谱S。在某些实施例中，由计算机程序436的指令实现的频谱估计包括实现构造最小二乘(L-S)矩阵A的指令和实现求解L-S拟合问题以利用所构造的L-S矩阵A找出所估计的频谱S的指令。所实现的频谱估计还包括计算机程序436的实现计算由L-S拟合所产生的所估计频谱S的信噪比(SNR)的指令以及实现从数据流D中选择附加的样本以添加到所构造的L-S矩阵A直到实现目标或预定SNR的指令为止。

具体而言，在某些实施例中，在L-S问题求解中采用了L-S矩阵A的QR因子分解。采用关于L-S矩阵A的QR因子分解的R因子矩阵的条件来计算是否实现了目标SNR。在某些这样的实施例中，该条件由上述方程(8)给出。在其他实施例中，采用由上述方程(9)给出的条件。在其他实施例中，采用因子矩阵R和L-S矩阵A中的一个或两个的条件数。

在某些实施例中，构造和计算与以上针对频谱估计方法100所描述的构造110本质上类似。具体而言，计算机程序436的指令实现了从数据流D中选择初始样本集合。在某些实施例中，所实现的选择与以上针对频谱估计方法100的构造110所描述的选择112本质上类似。这些指令还包括利用感兴趣的频率F和初始集合中样本的时间戳t_n来实现L-S矩阵A的初始化。在某些实施例中，所实现的初始化与以上针对频谱估计方法100的构造110所描述的初始化114本质上类似。这些指令还包括实现对经初始化的L-S矩阵A执行QR因子分解，以生成因子矩阵Q和R。在某些实施例中，所实现的执行与以上针对频谱估计方法100的构造110所描述的执行116本质上类似。这些指令还包括利用一个或多个附加样本更新L-S拟合直到频谱估计S实现由关于因子矩阵R的条件所给出的目标SNR为止。在某些实施例中，所实现的更新与以上针对频谱估计方法100的构造110所描述的更新118本质上类似。

在某些实施例中，计算机程序436还包括实现求解L-S拟合问题以获得傅立叶变换系数的指令。具体而言，可以采用上述方程(11)形式的方程来解出傅立叶变换系数。具体而言，在某些实施例中，所实现的求解可以与以上针对频谱估计方法100所描述的求解120本质上类似。

在某些实施例中，计算机程序436还包括实现提取傅立叶系数以产生所估计的频谱S的指令。所实现的提取可以与以上针对频谱估计方法100所描述的提取130本质上类似。具体而言，提取可以产生以感兴趣的频率F为索引的实系数和虚系数形式(c_m，s_m)或幅度系数和相位系数形式(A_m，φ_m)的所估计频谱S。

在频谱估计器430的其他实施例中，计算机程序436包括实现从带时间戳的样本的数据流中选择第一样本集合的指令，以及实现随后选择带时间戳样本的其他集合以逐步与第一集合相组合来创建组合样本集合的指令。其他集合被逐步组合，直到组合样本集合的所估计频谱的信噪比的量度指示实现了目标信噪比为止。在某些实施例中，选择第一样本集合与上述方法200中的选择220本质上类似。此外，在某些实施例中，选择带时间戳的样本的其他集合以逐步组合与也是上述方法200的选择220其他样本集合本质上类似。

在频谱估计器430中的另一个实施例中，计算机程序436包括实现从所采样的数据中选择带时间戳的样本以产生样本集合并且计算该样本集合的所估计频谱的信噪比量度的指令。计算机程序436还包括实现逐步选择附加的带时间戳样本以与该样本集合相组合并且重新计算组合样本集合的所估计频谱的信噪比量度直到该量度小于预定限度λ为止的指令。当信噪比量度小于预定限度λ时，就实现了目标信噪比。

在某些实施例中，选择带时间戳的样本与上述方法300中的选择310本质上类似。在某些实施例中，计算信噪比量度与也是上述方法300中的计算320本质上类似。此外，在某些实施例中，逐步选择附加的带时间戳样本以与该样本集合相组合与上述方法300中的逐步选择330基本上类似，而重新计算信噪比量度与也是上述方法300中的重新计算340基本上类似。

返回参考图5，在某些实施例中，测试系统400还包括测试结果评估器440。测试结果评估器440接收来自频谱估计器430的所估计频谱S。接收到的所估计频谱S被与DUT 402的预定频谱标准或规格相比较。如果DUT 402的所估计频谱S符合标准，则结果评估器440将通过结果分配到DUT 402。如果所估计频谱S未能符合一个或多个标准，则失败结果被结果评估器440分配给DUT 402。在某些实施例中，(未示出)结果评估器440被实现为频谱估计器430的计算机程序436的一部分。

从而，已经描述了频谱估计方法、采用频谱估计的测试系统和频谱估计器。应当理解，上述实施例只是对代表本发明原理的许多特定实施例中的某些的说明。很明显，本领域的技术人员可以很容易地设计出许多其他配置，而不会脱离由所附权利要求书限定的本发明的范围。

Claims (10)

1.一种频谱估计器430，包括

计算机处理器432；

存储器434；以及

存储在存储器434中并被计算机处理器432执行的计算机程序436，

其中所述计算机程序436包括这样的指令，这种指令在被所述处理器执行时，实现从带时间戳的样本的数据流中选择110、210、310第一样本集合，还实现随后选择110、220、330带时间戳的样本的其他集合以逐步与所述第一集合相组合以便创建组合样本集合，直到所述组合样本集合的所估计的频谱100、200、300的信噪比的量度指示实现110、220、320、340目标信噪比为止。

2.如权利要求1所述的频谱估计器430，其中，当所述信噪比的量度小于预定限度λ时，实现110、220、320、340所述目标信噪比，所述信噪比的量度与所述组合样本集合到所述所估计的频谱100、200、300的最小二乘拟合的最小二乘(L-S)矩阵A的条件数成比例。

3.如权利要求1-2中任何一个所述的频谱估计器430，其中所述被执行的指令还实现确定110、220、320、340所述组合样本集合到所述所估计的频谱100、200、300的最小二乘拟合的最小二乘(L-S)矩阵A的QR因子分解的右上三角因子矩阵R，当所述因子矩阵R非奇异时，所述信噪比的量度是所述L-S矩阵A的QR因子分解的所选对角线元素的比率。

4.如权利要求3所述的频谱估计器430，其中所述所选对角线元素的比率是以下比率中的一种或两种：从所述因子矩阵R的主对角线的按元素绝对值中选择出来的元素与从所述主对角线的按元素绝对值中选择出来的最小值元素的最大值比率，以及作为从所述因子矩阵R的主对角线的按元素绝对值中选择出来的最小值元素的倒数的比率，并且其中所述预定限度λ等于或大于约1.1，并且小于约1,000。

5.如权利要求1-3中任何一个所述的频谱估计器430，其中，当所述组合样本集合到所述所估计的频谱100、200、300的最小二乘拟合的最小二乘(L-S)矩阵A的QR因子分解的右上三角因子矩阵R符合由下式给出的条件时，实现110、220、320、340所述目标信噪比：

min(abs(diag(R)))＞0∧1/min(abs(diag(R)))＜λ

其中λ是等于或大于约1.1并且小于约1,000的预定限度。

6.如权利要求1-3和5中任何一个所述的频谱估计器430，其中，当所述组合样本集合到所述所估计的频谱100、200、300的最小二乘拟合的最小二乘(L-S)矩阵A的QR因子分解的右上三角因子矩阵R符合由下式给出的条件时，实现110、220、320、340所述目标信噪比：

min(abs(diag(R)))＞0∧max(abs(diag(R)))/min(abs(diag(R)))＜λ

其中λ是等于或大于约1.1并且小于约1,000的预定限度。

7.如权利要求1-6中任何一个所述的频谱估计器430，其中所述数据流包括随机采样的数据和不规则采样的数据中的一种或两种。

8.如权利要求1-7中任何一个所述的频谱估计器430，其中所述被执行的指令还实现通过求解120所述组合集合的最小二乘拟合问题计算所述所估计的频谱100、200、300；并且从所述最小二乘拟合问题的解答中提取130所述所估计的频谱100、200、300的傅立叶系数。

9.如权利要求1-8中任何一个所述的频谱估计器430，被用于采用频谱估计100、200、300的系统400中，所述系统400包括：

激励源410，其连接到被测设备402的输入；以及

采样器420，其连接到所述被测设备402的输出，所述采样器420具有连接到所述频谱估计器430的输出，所述采样器420将来自所述被测设备402的所采样的数据提供给所述频谱估计器430。

10.如权利要求9所述的频谱估计器430，其中所述系统400还包括连接到所述频谱估计器430的输出的测试结果评估器440，所述测试结果评估器440确定所述被测设备402是否符合预定规格。

Images