对一组信号使用 pwelch：一些问题（Matlab）答案

【问题标题】：Using pwelch to a set of signals: some questions (Matlab)对一组信号使用 pwelch：一些问题（Matlab）
【发布时间】：2016-02-21 05:27:44
【问题描述】：

我想对一组信号使用 pwelch，但我有一些问题。

首先，假设我们有 32 个持续时间为 30 秒的 (EEG) 信号。采样频率为fs=256samples/sec，因此每个信号的长度为 7680。我想使用pwelch 来估计这些信号的功率谱密度 (PSD)。

问题 1： 基于pwelch的documentation，

pxx = pwelch(x) 返回使用 Welch 的重叠段平均估计器找到的输入信号 x 的功率谱密度 (PSD) 估计值 pxx。当 x 是向量时，它被视为单个通道。当x为矩阵时，对每一列独立计算PSD并存储在pxx的对应列中。

但是，如果调用pwelch如下

% ch_signals: 7680x32; one channel signal per each column
[pxx,f] = pwelch(ch_signals);

生成的pxx 的大小为1025x1，而不是我所期望的1025x32，因为文档指出，如果 x 是矩阵，则 PSD 会为每一列独立计算并存储在 pxx 的相应列中.

问题 2： 假设我克服了这个问题，并且我独立地计算了每个信号的 PSD（通过将pwelch 应用于ch_signals 的每一列），我想知道这样做的最佳方法是什么。假设信号是一个 30 秒的信号，采样频率为 fs=256，我应该如何调用 pwelch（带有什么参数？）这样 PSD 才有意义？

问题 3：如果我需要将我的 32 个信号中的每一个拆分为多个窗口并将 pwech 应用于每个窗口，那么最好的方法是什么？假设我想将每个 30 秒的信号分成 3 秒的窗口，重叠 2 秒。我应该如何为每个窗口调用pwelch？

【问题讨论】：

该文档适用于 R2015 测试版，您的可能不是 R2015。例如我的是 R2012，当x 是一个矩阵时它什么也没说！ pwelch(x,window) 中的 window 不是你的意思吗？
@Rashid 文档指的是最新版本，即 R2015b，其中 b 代表今年的第二个版本。该方案自 2006 年开始使用。
@Rashid，正确的，我在 Matlab R2014a 上，并且 - 确实 - 它没有说明这个案例。谢谢！所以，在这种情况下，我只是遍历通道，这完全没问题。我真正的问题是关于问题 2 和 3：在这些情况下我应该如何应用 pwelch（参数？）。特别是在第二种情况（Q3）中，信号需要事先在多个窗口中进行拆分。这里的窗口数量越大越好。非常感谢您的回复！
在 R2014b 中，他们引入了对 pwelch 和其他功能 (Reference) 的多通道支持。 Q2 和 Q3 并不是真正关于编程，而是更多关于参数化。或许可以在SPSE 中提问以获得具体答案。
感谢 cmets @Matt

标签： matlab spectral-density

【解决方案1】：

这是一个例子，就像你的情况一样，

结果表明，该算法指示的信号频率恰到好处。

矩阵的每一列，y 都是一个正弦曲线来检查它是如何工作的。

窗口为 3 秒，重叠 2 秒，

Fs = 256;                       
T = 1/Fs;                    
t = (0:30*Fs-1)*T;           
y = sin(2 * pi * repmat(linspace(1,100,32)',1,length(t)).*repmat(t,32,1))';
for i = 1 : 32
[pxx(:,i), freq] = pwelch(y(:,i),3*Fs,2*Fs,[],Fs); %#ok
end
plot(freq,pxx);
xlabel('Frequency (Hz)'); 
ylabel('Spectral Density (Hz^{-1})');

【讨论】：

非常感谢您的回答！还有一件事。您的回答涉及我的第二个问题。第三个呢？如果我想在 3 秒的窗口中分割原始信号并有 2 秒的重叠，那么 pwelch 有什么用处？我的意思是每个窗口应该调用pwelch 什么？从某种意义上说，将pwelch 应用于 3 秒信号可能会导致较差的结果......可能是这种情况吗？
@nullgeppetto，你为什么不像我的帖子那样使用 3 秒的汉明窗口和 2 秒的重叠？
我不确定我是否理解你。让我们忘记 32 个信号。假设我们只有一个 30 秒的信号，fs=256。假设我不想将 pwelch 应用于整个信号，但我需要先将其拆分为多个窗口，然后将 pwelch 应用于每个窗口。那会有意义吗？我们这里有两种窗口，我知道这可能会令人困惑......
@nullgeppetto，您可以使用reshape(y,3*Fs,32,[]); 分割您的信号，然后两个for 循环可以解决问题。但我不确定这是否是你想要的。
@nullgeppetto，不客气，感谢您接受答案。