非参数方法
下面的章节将讨论<一个href="//www.tatmou.com/it/it/help/signal/ug/nonparametric-methods.html" class="intrnllnk">周期图一个>,<一个href="//www.tatmou.com/it/it/help/signal/ug/nonparametric-methods.html" class="intrnllnk">改进的周期图一个>,<一个href="//www.tatmou.com/it/it/help/signal/ug/nonparametric-methods.html" class="intrnllnk">韦尔奇一个>,<一个href="//www.tatmou.com/it/it/help/signal/ug/nonparametric-methods.html" class="intrnllnk">多窗口一个>非参数估计方法,以及相关的<一个href="//www.tatmou.com/it/it/help/signal/ug/nonparametric-methods.html" class="intrnllnk">运行CPSD函数一个>,<一个href="//www.tatmou.com/it/it/help/signal/ug/nonparametric-methods.html" class="intrnllnk">传递函数估计一个>,<一个href="//www.tatmou.com/it/it/help/signal/ug/nonparametric-methods.html" class="intrnllnk">相干函数一个>。
周期图
概括地说,一种估算PSD的过程就是找到样本的离散时间傅里叶变换的过程(通常在一个网格FFT)和适当的规模大小的平方的结果。这估计是所谓的<年代pan class="emphasis">周期图年代pan>。
PSD的周期图估计的一个信号<年代pan class="inlineequation"> 的长度<年代pan class="emphasis">l年代pan>是
Fs的采样频率。
在实践中,实际的计算<年代pan class="inlineequation"> 可以只在有限数量的执行频率点,通常采用FFT。大多数计算周期图的实现方法<年代pan class="inlineequation"> 分PSD估计的频率
在某些情况下,周期图的计算通过FFT算法更高效的频率是2的幂。因此并不少见垫与零输入信号,延长其长度为2的幂。
周期图作为一个例子,考虑以下1001 -元素信号xn
,它由两个正弦信号加噪声:
fs = 1000;<年代pan style="color:#228B22">%采样频率年代pan>t = (0: fs) / fs;<年代pan style="color:#228B22">% 1秒的样本年代pan>一个= (1 - 2);<年代pan style="color:#228B22">%正弦信号振幅(行向量)年代pan>f = (150; 140);<年代pan style="color:#228B22">%正弦信号频率(列向量)年代pan>xn = A * sin(2 *π* f * t) + 0.1 * randn(大小(t));<年代pan style="color:#228B22">%的最后三行是等价的年代pan>% xn =罪(2 *π* 150 * t) + 2 *罪(2 *π* 140 * t) + 0.1 * randn(大小(t));年代pan>
PSD的周期图估计可以计算使用周期图
。在这种情况下,数据向量乘以一个汉明窗生产周期图修改。
[Pxx F] =周期图(xn汉明(长度(xn)),长度(xn), fs);情节(F, 10 * log10 (Pxx))包含(<年代pan style="color:#A020F0">“赫兹”年代pan>)ylabel (<年代pan style="color:#A020F0">“数据库”年代pan>)标题(<年代pan style="color:#A020F0">改进的周期图估计功率谱密度的年代pan>)
算法年代trong>
周期图计算和鳞片FFT产生的输出功率与频率图如下。
如果输入信号是实值,由此产生的FFT的大小是对称的零频率(DC)。的就是FFT,只有第一个(1 +
nfft
/ 2)点是独一无二的。确定惟一的值的数量,只保留那些独特的点。独特的FFT的平方大小的值。规模的平方大小(DC)除外<年代pan class="inlineequation"> ,在那里<年代pan class="emphasis">N年代pan>信号的长度是之前任何补零。规模的直流值<年代pan class="inlineequation"> 。
创建一个频率向量的数量的独特点,nfft和采样频率。
情节的大小的平方FFT频率。
周期图的性能
下面的章节将讨论周期图的性能对泄漏的问题,决议,偏见,和方差。
频谱泄漏年代trong>
考虑一个有限长度的PSD(长度<年代pan class="inlineequation"> )信号<年代pan class="inlineequation"> 。它通常是有用的解释<年代pan class="inlineequation"> 无限信号相乘的结果,<年代pan class="inlineequation"> 一个有限长矩形窗口,<年代pan class="inlineequation"> :
因为乘法在时域卷积对应于在频域中,周期图在频域的期望值
显示的期望值的周期图的卷积是真的PSD的平方狄利克雷内核。
卷积的效果是最好的理解为正弦数据。假设<年代pan class="inlineequation"> 由一笔吗<年代pan class="inlineequation"> 复杂的正弦曲线:
它的光谱
对于一个有限长序列变成了
等于前面的方程
所以在有限长信号的频谱,取代了狄拉克δ函数的形式<年代pan class="inlineequation"> ,对应于一个矩形窗口的频率响应为中心频率<年代pan class="inlineequation"> 。
一个矩形窗口的频率响应周期sinc的形状:
L = 32;[h, w] = freqz (rectwin(左)/ L, 1);y = diric (w、L);情节(w /π,20 * log10 (abs (h)))<年代pan style="color:#A020F0">在年代pan>情节(w /π,20 * log10 (abs (y)),<年代pan style="color:#A020F0">“——”年代pan>)举行<年代pan style="color:#A020F0">从年代pan>ylim([-40,0])传说(<年代pan style="color:#A020F0">的频率响应年代pan>,<年代pan style="color:#A020F0">“周期性Sinc”年代pan>)包含(<年代pan style="color:#A020F0">“ω\ / \π”年代pan>)
情节显示mainlobe和几个旁瓣,其中最大的是大约13.5分贝以下mainlobe高峰。这些叶占称为频谱泄漏的影响。而无限长的信号有其权力集中在离散频率点<年代pan class="inlineequation"> 视窗化(或截断)信号连续的权力“泄露”在离散频率点<年代pan class="inlineequation"> 。
因为短矩形窗口的频率响应是一个穷得多的近似狄拉克δ函数比一个时间窗口,谱泄漏数据记录短时尤其明显。考虑以下的100个样本序列:
fs = 1000;<年代pan style="color:#228B22">%采样频率年代pan>t = (0: fs / 10) / fs;<年代pan style="color:#228B22">%十分之一秒的样本年代pan>一个= (1 - 2);<年代pan style="color:#228B22">%正弦信号的振幅年代pan>f = (150; 140);<年代pan style="color:#228B22">%正弦信号频率年代pan>xn = A * sin(2 *π* f * t) + 0.1 * randn(大小(t));周期图(xn rectwin(长度(xn)), 1024年,fs)
重要的是要注意,频谱泄漏的影响完全取决于数据记录的长度。它不是一个事实的结果计算周期图频率有限数量的样本。
决议年代trong>
决议年代pan>是指区分光谱特征的能力,是一个关键的概念分析的谱估计的性能。
为了解决两个正弦信号频率相对接近,两者的区别是必要的频率要大于mainlobe泄露的光谱的宽度为其中一个正弦曲线。mainlobe宽度定义是mainlobe的宽度的权力是mainlobe权力(即峰值的一半。3 dB宽度)。这个宽度约等于<年代pan class="inlineequation"> 。
换句话说,两个正弦信号的频率<年代pan class="inlineequation"> 和<年代pan class="inlineequation"> 可分解性条件要求
在上面的示例中,两个正弦曲线仅仅相距10 Hz,必须大于100个样本数据记录决议允许两个截然不同的正弦曲线的周期图。
考虑这样一种情况:不满足这一标准,至于下面的67个样本序列:
fs = 1000;<年代pan style="color:#228B22">%采样频率年代pan>t = (0: fs / 15) / fs;<年代pan style="color:#228B22">% 67个样本年代pan>一个= (1 - 2);<年代pan style="color:#228B22">%正弦信号的振幅年代pan>f = (150; 140);<年代pan style="color:#228B22">%正弦信号频率年代pan>xn = A * sin(2 *π* f * t) + 0.1 * randn(大小(t));周期图(xn rectwin(长度(xn)), 1024年,fs)
上面的讨论决议没有考虑噪声的影响,因为信噪比(信噪比)迄今为止相对较高。当信噪比低,更难区分真正的光谱特性,噪声工件出现在光谱估计基于周期图。下面的例子说明了这一点:
fs = 1000;<年代pan style="color:#228B22">%采样频率年代pan>t = (0: fs / 10) / fs;<年代pan style="color:#228B22">%十分之一秒的样本年代pan>一个= (1 - 2);<年代pan style="color:#228B22">%正弦信号的振幅年代pan>f = (150; 140);<年代pan style="color:#228B22">%正弦信号频率年代pan>xn = A * sin(2 *π* f * t) + 2 * randn(大小(t));周期图(xn rectwin(长度(xn)), 1024年,fs)
偏见的周期图年代trong>
PSD的周期图是一个有偏估计量。它以前是期望值
渐近无偏周期图,这是明显的从前面观察,随着数据记录长度趋于无穷时,矩形窗口的频率响应更紧密地接近狄拉克δ函数。然而,在某些情况下,周期图是一个贫穷的估计量的PSD即使数据记录。这是由于周期图的方差,如下解释。
方差的周期图年代trong>
的方差可以证明是周期图
这表明方差不倾向于零的数据长度<年代pan class="inlineequation"> 趋向于无穷。在统计方面,PSD的周期图不一致的估计量。然而,周期图的谱估计可以是一个有用的工具情况下,信噪比高,特别是如果数据记录。
修改后的周期图
的<年代pan class="emphasis">改进的周期图年代pan>窗户前的时域信号计算DFT以平滑信号的边缘。这样做的效果是减少旁瓣的高度或频谱泄漏。这一现象产生的解释引入虚假频率信号的旁瓣突然截断,发生在一个矩形窗口。对于非矩形窗口,截断信号衰减顺利的结束点,因此引入的杂散频率不太严重。另一方面,也非矩形窗口扩大mainlobe,这导致减少决议。
的周期图
允许你计算一个改进的周期图通过指定窗口使用的数据。例如,比较一个默认的矩形窗和汉明窗。DFT指定相同数量的点在这两种情况下。
fs = 1000;<年代pan style="color:#228B22">%采样频率年代pan>t = (0: fs / 10) / fs;<年代pan style="color:#228B22">%十分之一秒的样本年代pan>一个= (1 - 2);<年代pan style="color:#228B22">%正弦信号的振幅年代pan>f = (150; 140);<年代pan style="color:#228B22">%正弦信号频率年代pan>nfft = 1024;xn = A * sin(2 *π* f * t) + 0.1 * randn(大小(t));周期图(xn rectwin(长度(xn)), nfft, fs)
周期图(xn汉明(长度(xn)), nfft, fs)
您可以验证,虽然不太明显的旁瓣Hamming-windowed周期图,两个主要的山峰是广泛的。事实上,3 dB mainlobe对应一个汉明窗的宽度大约是一个矩形窗口的两倍。因此,对于一个固定数据长度、PSD分辨率可以达到与汉明窗是大约一半,实现一个矩形窗口。mainlobe宽度和旁瓣高度的相互竞争的利益在一定程度上可以解决通过使用变量窗口如凯撒窗口。
不规则窗口影响信号的平均功率,因为有些时候样品衰减时乘以窗口。为了弥补这个缺陷,周期图
和pwelch
规范化的窗口有一个团结的平均功率。这确保了测量平均功率通常是独立于窗口的选择。如果频率成分并不好解决PSD估计,窗口选择并影响平均功率。
修改后的PSD的周期图估计
在哪里<年代pan class="emphasis">U年代pan>窗口归一化常数:
对于大的值l
,U
往往成为独立的窗口长度。添加U
归一化常数确保修改后的周期图是渐近无偏。
韦尔奇的方法
一种改进的估计量的PSD韦尔奇提出的。的方法包括将时间序列数据划分为(可能是重叠的)段,计算每个部分的修改周期图,然后平均PSD的估计。韦尔奇的PSD估计结果。工具箱函数pwelch
韦尔奇的实现方法。
修改后的平均周期图往往减少估计的方差相对于一个单一的周期图估计整个数据记录。虽然重叠部分引入了冗余信息,这种效应被非矩形窗口的使用减少,从而减少了或重要性<年代pan class="emphasis">重量年代pan>给的最终样品段(重叠)的样本。
然而,正如上面提到的,短的结合使用数据记录和非矩形窗口的结果在降低分辨率的估计量。总之,减少方差和分辨率之间的权衡。可以操纵参数在韦尔奇的方法获得改善相对于周期图估计,特别是在信噪比很低。这是下面的示例中所示。
考虑一个信号包含301个样本:
fs = 1000;<年代pan style="color:#228B22">%采样频率年代pan>t = (0:0.3 * fs) / fs;<年代pan style="color:#228B22">% 301个样本年代pan>一个= 8 [2];<年代pan style="color:#228B22">%正弦信号振幅(行向量)年代pan>f = (150; 140);<年代pan style="color:#228B22">%正弦信号频率(列向量)年代pan>xn = A * sin(2 *π* f * t) + 5 * randn(大小(t));周期图(xn rectwin(长度(xn)), 1024年,fs)
我们可以获得3段韦尔奇的谱估计有50%的重叠使用矩形窗。
pwelch (xn rectwin(150), 50512年,fs)
在上面的周期图中,噪声和正弦信号的泄漏使人本质上的人工山峰。相比之下,虽然产生的PSD韦尔奇的方法具有更广泛的山峰,你仍然可以区分两个正弦曲线,脱颖而出“噪声地板上。”
然而,如果我们试图进一步降低方差,分辨率的损失原因的一个正弦曲线完全丢失。
pwelch (xn rectwin(100), 75512年,fs)
在韦尔奇的方法偏差和规范化
韦尔奇的PSD的方法产生一个有偏估计量。PSD的期望值估计是:
在哪里
韦尔奇的估计量的方差是很难计算,因为它取决于所使用的窗口和重叠部分的数量。基本上,方差是段的数量成反比的平均周期图被修改。
多窗口方法
周期图可以解释为滤波长度<年代pan class="inlineequation"> 信号,<年代pan class="inlineequation"> 通过滤波器组(一组并行过滤器)<年代pan class="inlineequation"> 冷杉带通滤波器。3 dB的带宽的带通滤波器可以约等于<年代pan class="inlineequation"> 。每一个带通滤波器的大小反应类似于一个矩形窗口。的周期图可以被看作是一个计算的力量(即每个过滤信号。,the output of each bandpass filter) that uses just one sample of each filtered signal and assumes that the PSD of<年代pan class="inlineequation"> 在每一个带通滤波器的带宽是恒定的。
随着信号的长度增加,每一个带通滤波器的带宽减少,使其更有选择性滤波器,改善常数的近似PSD滤波器的带宽。这提供了另一种解释为什么PSD的周期图估计提高随着信号的长度增加。然而,从这个角度来看,有两个因素明显妥协周期图估计的准确性。首先,矩形窗口收益率可怜的带通滤波器。其次,计算每个带通滤波器的输出功率的依赖于单个样本的输出信号,产生一个非常粗糙的近似。
韦尔奇的方法可以得到一个类似的诠释一个过滤器银行。在韦尔奇的实现中,几个样品是用来计算输出功率,从而降低方差的估计。另一方面,每一个带通滤波器的带宽大于相应的周期图方法,导致损失的决议。滤波器组模型从而提供了一种新的解释的方差和分辨率之间的妥协。
汤普森的<年代pan class="emphasis">多窗口方法年代pan>(MTM)基于这些结果提供一种改进的PSD估计。而不是使用带通滤波器,本质上是矩形窗口(如周期图方法),最理想的带通滤波器的MTM方法使用银行计算估计。这些优化FIR滤波器来自一组序列称为离散长球状序列(也称为离散长<年代pan class="emphasis">斯莱皮恩序列年代pan>)。
此外,MTM的方法提供了一个时间带宽参数来平衡方差和决议。这个参数是由时间带宽积,<年代pan class="inlineequation"> 它直接关系到使用的蜡烛数量计算频谱。总有<年代pan class="inlineequation"> 蜡烛用来估计形式。这意味着,<年代pan class="inlineequation"> 增加,有更多的功率谱的估计,估计的方差减少。然而,每个锥的带宽也是成正比的<年代pan class="inlineequation"> ,<年代pan class="inlineequation"> 增加,更多的谱泄漏(即每个估计展品。,wider peaks) and the overall spectral estimate is more biased. For each data set, there is usually a value for<年代pan class="inlineequation"> 允许偏差和方差之间的最佳平衡。
信号处理工具箱™函数,实现了MTM方法pmtm
。使用pmtm
计算的PSD信号。
fs = 1000;<年代pan style="color:#228B22">%采样频率年代pan>t = (0: fs) / fs;<年代pan style="color:#228B22">% 1秒的样本年代pan>一个= (1 - 2);<年代pan style="color:#228B22">%正弦信号的振幅年代pan>f = (150; 140);<年代pan style="color:#228B22">%正弦信号频率年代pan>xn = A * sin(2 *π* f * t) + 0.1 * randn(大小(t));pmtm (xn 4 [], fs)
通过降低时间带宽积,可以增加分辨率为代价较大的方差。
pmtm (xn 1.5, [], fs)
这种方法比韦尔奇的方法由于计算昂贵的成本计算离散长球状序列。长时间数据系列(10000分以上),它是有用的计算离散长一次,MAT-file拯救他们。dpsssave
,dpssload
,dpssdir
,dpssclear
继续提供的数据库保存在MAT-file离散长dpss.mat
。
交叉谱密度函数
PSD的一个特例<年代pan class="emphasis">交叉谱密度年代pan>(运行CPSD)函数,定义两个信号之间的关系
一样的相关性和协方差序列,工具箱<年代pan class="emphasis">估计年代pan>PSD和运行CPSD因为信号长度是有限的。
估计两个等长的信号的互谱密度x
和y
使用韦尔奇的方法,<一个href="//www.tatmou.com/it/it/help/signal/ref/cpsd.html">运行cpsd
函数形式的周期图的FFT的产物x
和FFT的共轭y
。不同于实值PSD,运行CPSD是一个复杂的功能。运行cpsd
处理的切片和窗口x
和y
在相同的方式pwelch
功能:
Sxy =运行cpsd (x, y, nwin、noverlap nfft, fs)
传递函数估计
韦尔奇的一个应用非参数系统辨识的方法。假设<年代pan class="emphasis">H年代pan>是一个线性、时不变系统和<年代pan class="emphasis">x年代pan>(<年代pan class="emphasis">n年代pan>),<年代pan class="emphasis">y年代pan>(<年代pan class="emphasis">n年代pan>)的输入和输出<年代pan class="emphasis">H年代pan>,分别。的功率谱<年代pan class="emphasis">x年代pan>(<年代pan class="emphasis">n年代pan>)的运行CPSD有关<年代pan class="emphasis">x年代pan>(<年代pan class="emphasis">n年代pan>),<年代pan class="emphasis">y年代pan>(<年代pan class="emphasis">n年代pan>)
估计之间的传递函数<年代pan class="emphasis">x年代pan>(<年代pan class="emphasis">n年代pan>),<年代pan class="emphasis">y年代pan>(<年代pan class="emphasis">n年代pan>)是
这种方法估计大小和相位信息。的tfestimate
函数使用韦尔奇的方法来计算运行CPSD和功率谱,然后形式传递函数的系数估计。使用tfestimate
以同样的方式使用运行cpsd
函数。
生成一个信号组成的两个正弦信号嵌入在高斯白噪声。
rng (<年代pan style="color:#A020F0">“默认”年代pan>fs = 1000;<年代pan style="color:#228B22">%采样频率年代pan>t = (0: fs) / fs;<年代pan style="color:#228B22">% 1秒的样本年代pan>一个= (1 - 2);<年代pan style="color:#228B22">%正弦信号的振幅年代pan>f = (150; 140);<年代pan style="color:#228B22">%正弦信号频率年代pan>xn = A * sin(2 *π* f * t) + 0.1 * randn(大小(t));
滤波器的信号xn
冷杉移动平均滤波器。计算实际的大小反应和估计响应。
h =(10) / 10的;<年代pan style="color:#228B22">%移动平均滤波器年代pan>yn =过滤器(h 1 xn);[命令f] = tfestimate (xn yn 256128256 fs);H = freqz (H 1 f, f);
策划的结果。
次要情节(2,1,1)情节(f, abs (H))标题(<年代pan style="color:#A020F0">实际传递函数级的年代pan>)yl = ylim;网格次要情节(2,1,2)情节(f, abs(命令))标题(<年代pan style="color:#A020F0">“传递函数估计级”年代pan>)包含(<年代pan style="color:#A020F0">的频率(赫兹)年代pan>)ylim (yl)网格
相干函数
两个信号之间的平方连贯性<年代pan class="emphasis">x年代pan>(<年代pan class="emphasis">n年代pan>),<年代pan class="emphasis">y年代pan>(<年代pan class="emphasis">n年代pan>)是
这个商是一个实数在0和1之间,措施之间的关系<年代pan class="emphasis">x年代pan>(<年代pan class="emphasis">n年代pan>),<年代pan class="emphasis">y年代pan>(<年代pan class="emphasis">n年代pan>)的频率<年代pan class="inlineequation"> 。
的mscohere
函数序列xn
和yn
计算功率谱和运行CPSD,并返回系数的平方级运行CPSD和产品的功率谱。它的选项和操作类似运行cpsd
和tfestimate
功能。
生成一个信号组成的两个正弦信号嵌入在高斯白噪声。在1 kHz信号采样1秒。
rng (<年代pan style="color:#A020F0">“默认”年代pan>fs = 1000;t = (0: fs) / fs;一个= (1 - 2);<年代pan style="color:#228B22">%正弦信号的振幅年代pan>f = (150; 140);<年代pan style="color:#228B22">%正弦信号频率年代pan>xn = A * sin(2 *π* f * t) + 0.1 * randn(大小(t));
滤波器的信号xn
冷杉移动平均滤波器。计算的相干函数和情节xn
和滤波器输出yn
作为频率的函数。
h =(10) / 10的;yn =过滤器(h 1 xn);mscohere (xn yn 256128256 fs)
如果输入序列长度,窗口长度,和重叠的数据点的数量在一个窗口是这样mscohere
只作用于一个记录,函数返回所有的。这是由于相干函数为线性相关的数据是1。
另请参阅
应用程序
功能
运行cpsd
|<年代pan itemscope itemtype="//www.tatmou.com/help/schema/MathWorksDocPage/SeeAlso" itemprop="seealso">mscohere
|<年代pan itemscope itemtype="//www.tatmou.com/help/schema/MathWorksDocPage/SeeAlso" itemprop="seealso">周期图
|<年代pan itemscope itemtype="//www.tatmou.com/help/schema/MathWorksDocPage/SeeAlso" itemprop="seealso">pmtm
|<年代pan itemscope itemtype="//www.tatmou.com/help/schema/MathWorksDocPage/SeeAlso" itemprop="seealso">pwelch
|<年代pan itemscope itemtype="//www.tatmou.com/help/schema/MathWorksDocPage/SeeAlso" itemprop="seealso">tfestimate