非参数方法-MATLAB和Simulink万博1manbetx - 万博1manbetx,s manbetx 845,万博尤文图斯

试用软件

产品更新

非参数方法

下面几节将讨论周期图,修正周期图,韦尔奇,多用户非参数估计方法，以及相关的CPSD函数,传递函数估计,相干函数．

周期图

打开生活的脚本

一般来说，估计过程PSD的一种方法是简单地找到过程样本的离散时间傅里叶变换（通常在网格上使用FFT进行），并适当地缩放结果的幅度平方。这个估计称为周期图．

信号PSD的周期图估计 $x_{L} (N)$ 的长度L是

$P_{x x} (F) = \frac{1.}{L F_{s}} {| \sum_{N = 0}^{L - 1.} x_{L} (N) E^{- J 2. π F N / F_{s}} |}^{2.},$

其中Fs为采样频率。

在实践中，实际计算 $P_{x x} (F)$ 只能在有限的频率点上执行，通常采用FFT。周期图方法的大多数实现计算 $N$ -频率点PSD估计

$F_{K} = \frac{K F_{s}}{N}, K = 0, 1., \dots, N - 1. ．$

在某些情况下，如果频率数为二次方，则通过FFT算法计算周期图的效率更高。因此，用零填充输入信号以将其长度扩展到二次方并不少见。

作为周期图的一个例子，考虑以下1001个元素信号xn，由两个正弦波和噪声组成：

fs = 1000;%采样频率t=（0:fs）/fs；%一秒钟的样品A = [1 2];%正弦振幅（行向量）f=[150；140]；%正弦频率(列向量)xn=A*sin（2*pi*f*t）+0.1*randn（尺寸（t））；%最后三行相当于%xn=sin（2*pi*150*t）+2*sin（2*pi*140*t）+0.1*randn（尺寸（t））；

PSD的周期图估计值可以使用周期图．在这种情况下，数据向量乘以一个汉明窗口来产生一个修改的周期图。

[Pxx，F]=周期图（xn，hamming（长度（xn）），长度（xn），fs）；绘图（F，10*log10（Pxx））xlabel(“赫兹”) ylabel (“dB”)头衔(“修正周期图功率谱密度估计”)

算法

周期图计算并缩放FFT的输出，产生如下功率与频率图。

如果输入信号是实值的，则得到的FFT的幅度相对于零频率（DC）是对称的。对于偶数长度FFT，仅第一个（1+非规则采样快速傅里叶变换/2)点是唯一的。确定唯一值的数量，并只保留那些唯一的点。
取唯一FFT值的平方大小。将平方大小（DC除外）按 $2. / (F_{s} N)$ 哪里N是任何零填充之前的信号长度。按比例缩放DC值 $1. / (F_{s} N)$ ．
根据唯一点的数量、nfft和采样频率创建频率向量。
根据频率绘制得到的幅度平方FFT。

周期图的性能

打开生活的脚本

下面几节讨论关于泄漏、分辨率、偏差和方差问题的周期图的性能。

频谱泄漏

考虑有限长度的PSD（长度） $L$ )信号 $x_{L} (N)$ .口译通常很有用 $x_{L} (N)$ 作为无限信号相乘的结果， $x (N)$ ，通过一个有限长的矩形窗口， $W_{R} (N)$ :

$x_{L} (N) = x (N) W_{R} (N) ．$

因为时域中的乘法对应于频域中的卷积，所以频域中周期图的期望值为

$E {{P_{}^{ˆ}}_{x x} (F)} = \frac{1.}{F_{s}} \int_{- F_{s} / 2.}^{F_{s} / 2.} \frac{罪^{2.} (L π (F - F^{'}) / F_{s})}{L 罪^{2.} (π (F - F^{'}) / F_{s})} P_{x x} (F^{'}) D F^{'},$

表明周期图的期望值是真实PSD与Dirichlet核平方的卷积。

对于正弦数据，卷积的影响是最容易理解的 $x (N)$ 由以下各项组成： $M$ 复杂的正弦曲线:

$x (N) = \sum_{K = 1.}^{N} {A.}_{K} E^{J ω_{K} N} ．$

它的光谱是

$X (ω) = \sum_{K = 1.}^{N} {A.}_{K} δ (ω - ω_{K}),$

对于一个有限长的序列

$X (ω) = \int_{- π}^{π} (\sum_{K = 1.}^{N} {A.}_{K} δ (ε - ω_{K})) W_{R} (ω - ε) D ε ．$

前面的方程式等于

$X (ω) = \sum_{K = 1.}^{N} {A.}_{K} W_{R} (ω - ω_{K}) ．$

所以在有限长度信号的频谱中，狄拉克函数被形式项所取代 $W_{R} (ω - ω_{K})$ ，对应于以频率为中心的矩形窗口的频率响应 $ω_{K}$ ．

矩形窗口的频率响应具有周期sinc的形状：

L = 32;[h, w] = freqz (rectwin(左)/ L, 1);y = diric (w、L);情节(w /π,20 * log10 (abs (h)))在图（w/pi，20*log10（绝对值（y）），“——”)举行从ylim（[-40,0]）图例(的频率响应,“周期性Sinc”)xlabel(“\omega/\pi”)

该图显示了一个主瓣和几个副瓣，其中最大的一个副瓣比主瓣峰值低约13.5 dB。这些瓣解释了所谓的频谱泄漏效应。而无限长信号的功率恰好集中在离散频率点 $F_{K}$ ，加窗(或截断)信号在离散频率点周围有连续的功率“泄漏” $F_{K}$ ．

由于短矩形窗口的频率响应比长窗口的频率响应更接近Dirac delta函数，因此当数据记录较短时，频谱泄漏尤其明显。考虑以下100个样本序列：

fs = 1000;%采样频率t = (0: fs / 10) / fs;%十分之一秒值的样品A = [1 2];%正弦振幅f=[150；140]；%正弦频率xn=A*sin（2*pi*f*t）+0.1*randn（大小（t））；周期图（xn，rectwin（长度（xn）），1024，fs）

值得注意的是，频谱泄漏的影响仅取决于数据记录的长度。这不是周期图在有限数量的频率样本下计算的结果。

决议

决议是指识别光谱特征的能力，是分析光谱估计器性能的关键概念。

为了分辨频率相对接近的两个正弦波，两个频率之间的差值必须大于其中任何一个正弦波泄漏频谱的主瓣宽度。主瓣宽度定义为功率为1/2时的主瓣宽度峰值主瓣功率（即3 dB宽度）。该宽度约等于 $F_{s} / L$ ．

换句话说，对于两个频率的正弦波 $F_{1.}$ 和 $F_{2.}$ ，可分辨性条件要求

$F_{2.} - F_{1.} > \frac{F_{s}}{L} ．$

在上面的例子中，两个正弦波仅以10赫兹分开，数据记录必须大于100个样本，以允许两个不同的正弦波通过周期图的分辨率。

考虑不满足此标准的情况，如以下67个样本序列:

fs = 1000;%采样频率t = (0: fs / 15) / fs;%67个样本A = [1 2];%正弦振幅f=[150；140]；%正弦频率xn=A*sin（2*pi*f*t）+0.1*randn（大小（t））；周期图（xn，rectwin（长度（xn）），1024，fs）

由于噪声信噪比（SNR）一直比较高，所以关于分辨率的讨论没有考虑噪声的影响，当信噪比较低时，真实的谱特征难以区分，并且基于周期图的谱估计中出现了噪声伪影。

fs = 1000;%采样频率t = (0: fs / 10) / fs;%十分之一秒值的样品A = [1 2];%正弦振幅f=[150；140]；%正弦频率xn=A*sin（2*pi*f*t）+2*randn（大小（t））；周期图（xn，rectwin（长度（xn）），1024，fs）

周期图的偏差

周期图是PSD的有偏估计量。它的期望值是

$E {{P_{}^{ˆ}}_{x x} (F)} = \frac{1.}{F_{s}} \int_{- F_{s} / 2.}^{F_{s} / 2.} \frac{罪^{2.} (L π (F - F^{'}) / F_{s})}{L 罪^{2.} (π (F - F^{'}) / F_{s})} P_{x x} (F^{'}) D F^{'} ．$

周期图是渐近无偏的，这从早期观察中可以明显看出，随着数据记录长度趋于无穷大，矩形窗口的频率响应更接近狄拉克δ函数。然而，在某些情况下，即使数据记录很长，周期图也不能很好地估计PSD。这是由于周期图的变化，如下所述。

周期图的方差

周期图的方差可以表示为

$v A. R ({P_{}^{ˆ}}_{x x} (F)) = {\begin{array}{c} P_{x x}^{2.} (F), & 0 < F < F_{s} / 2., \\ 2. P_{x x}^{2.} (F), & F = 0, F_{s} / 2., \end{array}$

这表明方差不会随着数据长度而趋于零 $L$ 趋于无穷大。在统计术语中，周期图不是PSD的一致估计值。然而，在信噪比较高的情况下，尤其是在数据记录较长的情况下，周期图可以作为谱估计的有用工具。

修正周期图

打开生活的脚本

这个修正周期图在计算DFT之前打开时域信号的窗口，以便平滑信号的边缘。这具有降低旁瓣高度或频谱泄漏的效果。这种现象导致将旁瓣解释为使用矩形窗口时发生的突然截断引入信号的杂散频率。对于非矩形窗口，截断信号的端点平滑衰减，因此引入的杂散频率不太严重。另一方面，非矩形窗口也会使主瓣变宽，从而导致分辨率降低。

这个周期图允许您通过指定要在数据上使用的窗口来计算修改后的周期图。例如，比较默认矩形窗口和Hamming窗口。在两种情况下指定相同数量的DFT点。

fs = 1000;%采样频率t = (0: fs / 10) / fs;%十分之一秒值的样品A = [1 2];%正弦振幅f=[150；140]；%正弦频率nfft=1024；xn=A*sin（2*pi*f*t）+0.1*randn（大小（t））；周期图（xn，rectwin（长度（xn）），nfft，fs）

周期图（xn，hamming（长度（xn）），nfft，fs）

你可以证实，虽然旁瓣在汉明窗周期图中不那么明显，但两个主峰更宽。事实上，与汉明窗相对应的主瓣宽度约为矩形窗的两倍。因此，对于固定的数据长度，使用汉明窗口可获得的PSD分辨率大约是矩形窗口的一半。利用凯撒窗等可变窗口，可以在一定程度上解决主瓣宽度和副瓣高度之间的相互冲突问题。

非矩形加窗影响信号的平均功率，因为当与窗相乘时，一些时间样本会衰减。为了弥补这一点，周期图和普韦奇使窗口具有统一的平均能力。这确保了测量的平均功率一般独立于窗口的选择。如果PSD估计器不能很好地解析频率分量，窗口的选择确实会影响平均功率。

PSD的修正周期图估计为

${P_{}^{ˆ}}_{x x} (F) = \frac{{| X (F) |}^{2.}}{F_{s} L U},$

哪里U为窗口归一化常数:

$U = \frac{1.}{L} \sum_{N = 0}^{N - 1.} {| W (N) |}^{2.} ．$

对于较大的L,U倾向于独立于窗口长度。添加U作为一个归一化常数，保证了修正的周期图是渐近无偏的。

韦尔奇方法

打开生活的脚本

Welch提出了一种改进的PSD估计方法。该方法包括将时间序列数据划分为（可能重叠的）段，计算每个段的修正周期图，然后对PSD估计进行平均。结果是Welch的PSD估计。工具箱函数普韦奇韦尔奇的实现方法。

相对于整个数据记录的单个周期图估计，修改周期图的平均值倾向于减小估计值的方差。虽然段之间的重叠引入了冗余信息，但使用非矩形窗口会减少这种影响，从而降低重要性或可用性重量给定段的结束样本（重叠的样本）。

然而，如上所述，短数据记录和非矩形窗口的组合使用会导致估计器的分辨率降低。总之，方差减少和分辨率之间存在折衷。可以操纵Welch方法中的参数，以获得相对于周期图的改进估计，尤其是当NR较低。下面的示例说明了这一点。

考虑一个由301个样本组成的信号：

fs = 1000;%采样频率t=（0:0.3*fs）/fs；%301个样本A = [2 8];%正弦振幅（行向量）f=[150；140]；%正弦频率(列向量)xn=A*sin（2*pi*f*t）+5*randn（大小（t））；周期图（xn，rectwin（长度（xn）），1024，fs）

我们可以用一个矩形窗口得到3段50%重叠的Welch谱估计值。

普韦尔奇（xn，rectwin（150），50512，fs）

在上面的周期图中，噪声和泄漏使得其中一个正弦曲线基本上无法与人工峰值区分开来。相比之下，虽然韦尔奇方法产生的PSD具有更宽的峰值，但您仍然可以区分两个正弦，这两个正弦与“噪声地板”不同

然而，如果我们试图进一步减小方差，分辨率的损失会导致其中一个正弦信号完全丢失。

pwelch (xn rectwin(100), 75512年,fs)

韦尔奇方法中的偏差与归一化

韦尔奇的方法产生了PSD的有偏估计量。PSD估计的期望值为:

$E {P_{韦尔奇} (F)} = \frac{1.}{F_{s} L U} \int_{- F_{s} / 2.}^{F_{s} / 2.} {| W (F - F^{'}) |}^{2.} P_{x x} (F^{'}) D F^{'},$

哪里L为数据段的长度，U修改后的周期图的定义中是否存在相同的归一化常数W（f）是窗函数的傅里叶变换。与所有周期图一样，Welch估计量是渐近无偏的。对于固定长度的数据记录，Welch估计的偏差大于周期图的偏差，因为段的长度小于整个数据样本的长度。

韦尔奇估计量的方差很难计算，因为它既取决于使用的窗口，也取决于段之间的重叠量。基本上，方差与被平均的修正周期图的段数成反比。

多窗口方法

打开生活的脚本

周期图可以解释为过滤长度 $L$ 信号 $x_{L} (N)$ ，通过 $L$ FIR带通滤波器。每个带通滤波器的3 dB带宽可以显示为大约等于 $F_{s} / L$ . 每个带通滤波器的幅度响应类似于矩形窗口的幅度响应。因此，可以将周期图视为每个滤波信号（即，每个带通滤波器的输出）的功率的计算，其仅使用每个滤波信号的一个样本，并且假设PSD为 $x_{L} (N)$ 在每个带通滤波器的带宽上是恒定的。

随着信号长度的增加，每个带通滤波器的带宽减小，使其成为一个更有选择性的滤波器，并改进了常数PSD在滤波器带宽上的近似。这提供了另一个解释，为什么周期图的PSD估计随着信号长度的增加而改善。然而，从这个观点来看，有两个明显的因素会损害周期图估计的准确性。首先，矩形窗口产生一个差的带通滤波器。其次，每个带通滤波器输出功率的计算依赖于输出信号的单个样本，产生一个非常粗略的近似。

Welch的方法可以从滤波器组的角度给出类似的解释。在Welch的实现中，使用多个样本来计算输出功率，从而减少了估计的方差。另一方面，每个带通滤波器的带宽大于对应于周期图方法的带宽，这导致了损耗因此，滤波器组模型为方差和分辨率之间的折衷提供了新的解释。

汤普森的多用户方法(MTM)建立在这些结果的基础上，提供改进的PSD估计。MTM方法不使用本质上是矩形窗口的带通滤波器(如周期图方法)，而是使用一组最优带通滤波器来计算估计值。这些最优FIR滤波器是由一组被称为离散长椭球序列(DPSSs，也称为斯莱皮恩序列)．

此外，MTM方法还提供了一个时间带宽参数来平衡方差和分辨率。这个参数由时间-带宽乘积给出， $N W$ 它与用于计算光谱的圆锥的数量直接相关 $2. N W - 1.$ 用于形成估算的锥度。这意味着 $N W$ 增加时，对功率谱有更多的估计，估计的方差减小。然而，每个锥度的带宽也成正比 $N W$ ，以便 $N W$ 增加时，每个估计显示出更多的频谱泄漏（即，更宽的峰值），并且总体频谱估计更具偏差。对于每个数据集，通常都有一个 $N W$ 这样就可以在偏差和方差之间进行最佳权衡。

信号处理工具箱™ 实现MTM方法的函数是pmtm．使用pmtm计算信号的PSD。

fs = 1000;%采样频率t=（0:fs）/fs；%一秒钟的样品A = [1 2];%正弦振幅f=[150；140]；%正弦频率xn=A*sin（2*pi*f*t）+0.1*randn（尺寸（t））；pmtm（xn，4，[]，fs）

通过降低时间带宽乘积，您可以以更大的方差为代价提高分辨率。

pmtm（xn，1.5，[]，fs）

由于计算离散长椭球序列的成本，该方法比Welch方法的计算成本更高。对于长数据序列（10000点或更多），一次计算DPSS并将其保存在MAT文件中非常有用。dpsssave,dpssload,dpssdir,dpssclear用于在MAT文件中保留已保存DPSS的数据库dpss.mat．

互谱密度函数

屏蔽门是屏蔽门的一个特例交叉谱密度（CPSD）功能，在两个信号之间定义x(N)及Y(N)作为

$P_{x Y} (ω) = \frac{1.}{2. π} \sum_{M = - \infty}^{\infty} R_{x Y} (M) E^{- J ω M} ．$

与相关和协方差序列的情况一样，工具箱估计PSD和CPSD，因为信号长度是有限的。

估计两个等长信号的互谱密度x和Y用韦尔奇的方法，我们得到了cpsd函数将周期图作为x以及FFT的共轭Y与实值PSD不同，CPSD是一个复杂函数。cpsd处理的剖面和窗口x和Y以同样的方式普韦奇功能:

Sxy=cpsd（x，y，nwin，noverlap，nfft，fs）

传递函数估计

打开生活的脚本

韦尔奇方法的一个应用是非参数系统辨识。假设H是一个线性、时不变系统x(N)及Y(N)是的输入和输出H分别地然后分析了系统的功率谱x(N)与的CPSD有关x(N)及Y(N)

$P_{Y x} (ω) = H (ω) P_{x x} (ω) ．$

二者之间传递函数的估计x(N)及Y(N)是

$H_{}^{ˆ} (ω) = \frac{{P_{}^{ˆ}}_{Y x} (ω)}{{P_{}^{ˆ}}_{x x} (ω)} ．$

这种方法估计幅度和相位信息。这个tfestimate函数使用Welch方法计算CPSD和功率谱，然后形成传递函数估计的商。使用tfestimate与您使用cpsd作用

产生一个信号，由嵌入在白色高斯噪声中的两个正弦信号组成。

rng (“默认”) fs = 1000;%采样频率t=（0:fs）/fs；%一秒钟的样品A = [1 2];%正弦振幅f=[150；140]；%正弦频率xn=A*sin（2*pi*f*t）+0.1*randn（尺寸（t））；

滤波器的信号xn使用FIR移动平均滤波器。计算实际幅值响应和估计响应。

h =(10) / 10的;%移动平均滤波器yn=滤波器（h，1，xn）；[HEST，f]=tfestimate（xn，yn，256128256，fs）；H=频率z（H，1，f，fs）；

绘制结果。

次要情节(2,1,1)情节(f, abs (H))标题(“实际传递函数大小”)yl=ylim；网格子地块（2,1,2）地块（f、abs（HEST））名称(“传递函数幅值估计”)xlabel(‘频率（Hz）’)ylim（yl）网格

相干函数

打开生活的脚本

两个信号之间的振幅平方相干性x(N)及Y(N)是

$C_{x Y} (ω) = \frac{{| P_{x Y} (ω) |}^{2.}}{P_{x x} (ω) P_{Y Y} (ω)} ．$

该商是一个介于0和1之间的实数，用于测量x(N)及Y(N)频率 $ω$ ．

这个mscohere函数接受序列xn和伊恩，计算它们的功率谱和CPSD，并返回CPSD的幅度平方与功率谱的乘积的商。其选项和操作类似于cpsd和tfestimate功能。

产生一个信号，由嵌入在白色高斯噪声中的两个正弦信号组成。信号以1khz采样1秒。

rng (“默认”) fs = 1000;t = (0: fs) / fs;A = [1 2];%正弦振幅f=[150；140]；%正弦频率xn=A*sin（2*pi*f*t）+0.1*randn（尺寸（t））；

滤波器的信号xn使用FIR移动平均滤波器。计算并绘制图像的相干函数xn和滤波器输出伊恩作为频率的函数。

h =(10) / 10的;yn =过滤器(h 1 xn);mscohere (xn yn 256128256 fs)

如果输入序列长度、窗口长度和窗口中重叠数据点的个数为mscohere只对一条记录进行操作，函数返回所有记录。这是因为线性相关数据的相干函数为1。

另见