如何在 MATLAB 中平均使用 Welch 方法估计的相干性

信息处理 matlab 功率谱密度 连贯性 平均
2022-02-13 11:52:32

我正在尝试估计 MATLAB 中信号不连续片段的平均相干性。对于 Welch 的频谱估计，以下工作：

w = blackmanharris(50);  % window
n = numel(w);            % fft length
o = 0;                   % overlap

x = rand(200,1); x1 = x(1:100); x2 = x(101:200);

s  = pwelch(x,  w, o, n);
s1 = pwelch(x1, w, o, n);
s2 = pwelch(x2, w, o, n);
sa = (s1+s2)/2;

max(abs(s-sa))  % almost zero

但是，当使用 Welch 方法估计相干性时，算术平均值会产生不正确的结果：

y = rand(200,1); y1 = y(1:100); y2 = y(101:200);

c  = mscohere(x,  y,  w, o, n);
c1 = mscohere(x1, y1, w, o, n);
c2 = mscohere(x2, y2, w, o, n);
ca = (c1+c2)/2;

max(abs(c-ca))  % very different from zero

如何正确平均两个相干估计？

1个回答

如果您查看mscohere包含幅度平方相干性定义的文档，您会看到它是交叉谱密度的绝对平方除以两个单独的谱密度的乘积。由于功率密度输入了表达式的分子和分母，因此来自完整时间序列的 MSC 估计不仅仅是部分估计的平均值。

要获得良好的“平均”MSC 估计，请执行以下操作：分别估计每个段中的交叉和单光谱密度，在段上平均这三个，然后将它们最终组合成一个平均 MSC。继续你的代码：

sx1 = pwelch(x1, w, o, n);
sy1 = pwelch(y1, w, o, n);
sxy1 = cpsd(x1, y1, w, o, n);

sx2 = pwelch(x2, w, o, n);
sy2 = pwelch(y2, w, o, n);
sxy2 = cpsd(x2, y2, w, o, n);

sxa = (sx1 + sx2) / 2;
sya = (sy1 + sy2) / 2;
sxya = (sxy1 + sxy2) / 2;
ca = abs(sxya) .^ 2 ./ sxa ./ sya;

现在max(abs(c-ca))几乎为零。

其它你可能感兴趣的问题
上一篇选择性地提高 FFT (DFT) 的频率分辨率 下一篇如何模拟磁带饱和度（音频 DSP）？