Using pwelch to a set of signals: some questions (

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?

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