MATLAB的答案

pwelch计算置信区间如何

16个视图(30天)
斯蒂芬。
斯蒂芬。 2020年4月30日
回答: 威廉•罗斯 2021年2月26日
Mathworks能指出哪些公式被用于计算置信区间在pwelch吗?也许国家代码中所使用的公式?这将包括exacty自由度的数量是如何被计算。
3评论
显示特殊照顾年长的评论
斯蒂芬。
斯蒂芬。 2021年2月15日
有一些好书在这本书——最初的詹金斯和瓦特很好(虽然绝版,而不是evailble pdf表单)是一本针对海洋学家汤森和金刚砂(数据分析方法在物理海洋学)。我最初问这个问题的原因是我怀疑计算置信区间由pwelch并不完全正确。可以检查这个通过光谱的随机过程有一个已知的频谱-詹金斯和瓦特给一个很好的例子的这样一个过程的书。我写了一个脚本——下面生成时间序列使用自回归过程J&W使用的例子。道歉为穷人编码技术(66年我第一次学会了用Fortran编程),您将看到,如果一个颠倒的CIs Matlab计算和交换上下CIs之间的结果,一个人看起来很不错,而如果一个棍棒与CIs pwelch,上下边界都高,这意味着下界不够低和上限太高了。这就是让我认为pwelch出事了。
%
%检查谱置信区间使用随机过程与一个已知的频谱
%通过AR1过程所使用的时间序列生成详细的詹金斯和瓦269页
% s Monismith
%斯坦福大学
% 2/15/21
%
清晰的
M = 16;
seriesize = 2 ^ M;%的时间长度
t = 0: seriesize-1;% *
X = 1 (seriesize, 1);%
Z = randn (seriesize, 1);
系数= 0.5;
为j = 3: seriesize
(j) = X (j - 1)因子* X (j2) + Z (j);
结束
(1)= X (seriesize)因子* X (seriesize-1) + Z (1);
(2)= X (1) X (seriesize) + Z因子* (2);
idothese = seriesize;
%
%关键调整——有多少段
%
nsegs = 2;
%
%现在设置和做pwelch估计
%
windowlength = idothese / nsegs;
nfft = windowlength;
[PXX f, PXXc] = pwelch (X, windowlength 0.5 * windowlength nfft, 1,“ConfidenceLevel”,0.95);
%
%计算边界——注意,他们现在乘法因素对已知的频谱
%,我们将使用pwelch给什么
%
因子1 =值(PXXc (: 1)。/ PXX);
factor2 =值(PXXc (:, 2)。/ PXX);
次要情节(1、2、1)
PXX semilogy (f,“柯”,“markeredgecolor”,0.7 * (1,1,1),…
“markerfacecolor”,0.7 * (1,1,1),“markersize”4)
持有在
Pxxtheory = var (X) * (0.834) / (2.25 - 3 * cos(2 *π* f) + cos(4 *π* f));
持有在
Pxxtheory semilogy (f,“r”,“线宽”,2)
semilogy (f, Pxxtheory *因子1,“b——”,“线宽”,2)
semilogy (f, Pxxtheory * factor2,“r——”,“线宽”,2)
持有从
ylim ([1 e 1 1 e2))
网格
集(gca),“字形大小”,20)
ylabel (“\ Phi_ {XX}’,“字形大小”,20)
包含(“f”,“字形大小”,20)
标题(“由pwelch独联体”,“字形大小”,20)
%
%现在重新置信区间的建议修正
%
因子1 =值(PXX. / PXXc (: 1));
factor2 =值(PXX. / PXXc (:, 2));
%
次要情节(1、2、2)
PXX semilogy (f,“柯”,“markeredgecolor”,0.7 * (1,1,1),…
“markerfacecolor”,0.7 * (1,1,1),“markersize”4)
持有在
Pxxtheory = var (X) * (0.834) / (2.25 - 3 * cos(2 *π* f) + cos(4 *π* f));
持有在
Pxxtheory semilogy (f,“r”,“线宽”,2)
semilogy (f, Pxxtheory *因子1,“b——”,“线宽”,2)
semilogy (f, Pxxtheory * factor2,“r——”,“线宽”,2)
持有从
ylim ([1 e 1 1 e2))
网格
集(gca),“字形大小”,20)
ylabel (“\ Phi_ {XX}’,“字形大小”,20)
包含(“f”,“字形大小”,20)
标题(独联体倒和交换,“字形大小”,20)

登录置评。

在回答这个问题。

答案(2)

Kiran Felix罗伯特
Kiran Felix罗伯特 2021年2月5日
你好斯蒂芬,
参考的 pwelch 函数文档细节。
2的评论
显示 1年长的评论
威廉•罗斯
威廉•罗斯 2021年2月15日
Kiran,
我读pwelch的源代码。米(C: \ Program Files \ MATLAB工具箱\ R2018b \ \ \信号\ pwelch.m)。它没有透露如何计算置信区间。pwelch。韦尔奇m调用。米(C: \ Program Files \ MATLAB工具箱\ R2018b \ \ \信号\ \ welch.m +光谱),我检查。它调用另一个版本的pwelch()我本来没有找到,可能因为它是在一个编译的库。因此我想欣赏斯蒂芬的问题实际的答案。谢谢你！

登录置评。

威廉•罗斯
威廉•罗斯 2021年2月26日
双向的置信区间(C.I.)与功率谱denisty概率p (p.s.d),由pwelch返回(),是由
Pxxhat (f) * k / chi2 ((1 + p) / 2, k) < Pxx (f) < Pxxhat (f) * k / chi2 ((1 - p) / 2, k)
Pxxhat (f)的实验估计p.s.d.在频率f和Pxx (f)是真的,但未知,p.s.d价值。,k是自由度。这是analagous方差的置信区间。自由度是由
k = 2 * k
大写字母K是段的数量或windows esitmating p.s.d时使用。
如果选择窗口长度和偏移量,以便windows信号没有任何剩余的部分,然后
K = (N-L) / D + 1
其中N =窦的长度L =窗口长度,D =偏移量。(不要混淆重叠和抵消。= L-D重叠。pwelch()以重叠作为参数,但韦尔奇(1967)使用抵消= D在他的公式。
例子:信号长度= 1024 =汉明窗,长度为256,half-overlapped。然后N = 1024, L = 256, D =抵消= 128。因此K = 7,由上面的公式,我们可以确认7一半完全重叠窗口覆盖信号。然后k = 2 * k = 14。
概率{Pxxhat * 14 / chi2 (.975 14) < Pxx < Pxxhat * 14 / chi2 (.025, 14)} = 0.95
= >概率{Pxxhat * 0.536 < Pxx < Pxxhat} = 0.95 * 2.487
我有检查上述公式的各种组合N, L、D和概率p。公式正确再现置信区间pwelch报道()。
上面的公式不正确从统计的角度来看,当重叠大于零,但它们的公式pwelch()使用。错误是,k,自由度,应该有点小于2 * k时窗口重叠。韦尔奇的精确公式, IEEE反式音频电声学 AU-15 1967年:70 - 73, https://ieeexplore.ieee.org/document/1161901 。它很复杂所以我不会复制在这里。它有一个变量的数量条款,根据窗口形状和数量的重叠。从pwelch C.I.()不考虑窗口形状或重叠。pwelch()的置信区间是正确当窗户不重叠。置信区间的错误报道pwelch()使用汉明时相对较小或损害或类似windows half-overlap。一般来说,C.I. pwelch报道()比它应该更窄。见下表:
置信区间的错误报道pwelch()变得更加严重,如果重叠更大。更严重的,我的意思是,从pwelch C.I.()比真正的C.I.更窄了,

类别

找到更多的在谱估计在帮助中心和文件交换

标签

s manbetx 845


释放

R2019a

社区寻宝

找到宝藏在MATLAB中央,发现社区如何帮助你!

开始狩猎!