Using pwelch to a set of signals: some questions (Matlab)

Question

I would like to use pwelch on a set of signals and I have some questions.

First, let's say that we have 32 (EEG) signals of 30 seconds duration. The sampling frequency is fs=256 samples/sec, and thus each signal has length 7680. I would like to use pwelch in order to estimate the power spectral density (PSD) of those signals.

Question 1: Based on the pwelch's documentation,

pxx = pwelch(x) returns the power spectral density (PSD) estimate, pxx, of the input signal, x, found using Welch's overlapped segment averaging estimator. When x is a vector, it is treated as a single channel. When x is a matrix, the PSD is computed independently for each column and stored in the corresponding column of pxx.

However, if call pwelch as follows

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

the resulting pxx is of size 1025x1, not 1025x32 as I would expect, since the documentation states that if x is a matrix the PSD is computed independently for each column and stored in the corresponding column of pxx.

Question 2: Let's say that I overcome this problem, and I compute the PSD of each signal independently (by applying pwelch to each column of ch_signals), I would like to know what is the best way of doing so. Granted that the signal is a 30-second signal in time with sampling frequency fs=256, how should I call pwelch (with what arguments?) such that the PSD is meaningful?

Question 3: If I need to split each of my 32 signals into windows and apply pwech to each one of those windows, what would be the best approach? Let's say that I would like to split each of my 30-second signals into windows of 3 seconds with an overlap of 2 seconds. How should I call pwelch for each one of those windows?

The document is for R2015 beta, probably yours isn't R2015. For example mine is R2012 and it says nothing when x is a matrix! The window in pwelch(x,window) isn't what you mean? — Rashid
@Rashid The documentation refers to the latest release which is R2015b where b stands for the second release of the year. This scheme is being used since 2006. — Matt
@Rashid, correct, I'm on Matlab R2014a and -indeed- it doesn't say anything about this case. Thanks! So, in this case I just iterate over channels, this is totally ok. What my issue really is concerns the questions 2 and 3: how should I apply pwelch in these cases (arguments?). Especially in the second case (Q3), where the signal needs to be split in a number of windows beforehand. Here the larger the number of windows, the better. Many thanks for your response! — nullgeppetto
In R2014b they introduced multichannel support for pwelch and other functions (Reference). Q2 and Q3 is not really regarding programming but more about the parametrization. Perhaps ask this in SPSE to get a specific answer. — Matt

Rashid Rashid · Accepted Answer · 2016-02-20T21:20:47

Here is an example, just like your case,

The results show that the algorithm indicates the signal frequencies just right.

Each column of matrix, y is a sinusoidal to check how it works.

The windows are 3 seconds with 2 seconds of overlapping,

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})');

Using pwelch to a set of signals: some questions (Matlab)

1 Answers