2

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?

nullgeppetto
  • 539
  • 1
  • 8
  • 25
  • 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 Feb 20 '16 at 21:03
  • @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 Feb 20 '16 at 21:08
  • @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 Feb 20 '16 at 21:16
  • In R2014b they introduced multichannel support for `pwelch` and other functions ([Reference](http://www.mathworks.com/help/signal/release-notes.html)). *Q2* and *Q3* is not really regarding programming but more about the parametrization. Perhaps ask this in [SPSE](http://dsp.stackexchange.com) to get a specific answer. – Matt Feb 20 '16 at 21:23
  • 1
    Thanks for comments @Matt – nullgeppetto Feb 20 '16 at 21:25

1 Answers1

1

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

enter image description here

Rashid
  • 4,326
  • 2
  • 29
  • 54
  • Thank you very much for your answer! Just one more thing. Your answer concerns my 2nd question. What about the third? What would be a good use of `pwelch` in case I want to split the original signal in windows of 3 seconds with 2 seconds overlap? I mean what should be the call of `pwelch` for each of those windows? In the sense that applying `pwelch` to a 3-second signal may lean to poor results... Could that be the case? – nullgeppetto Feb 20 '16 at 21:29
  • @nullgeppetto, Why don't you use 3sec windows of Hamming with 2sec overlapping, like in my post? – Rashid Feb 20 '16 at 21:32
  • I'm not sure I understand you. Let's forget about the 32 signals. Let's say we have just a single signal of 30 seconds with fs=256. Let's assume that I don't want to apply pwelch to the whole signal, but I need to split it first into windows and apply pwelch to each window. Would that be meaningful? We have two kinds of windows here, I know it may be confusing... – nullgeppetto Feb 20 '16 at 21:37
  • @nullgeppetto, you can use `reshape(y,3*Fs,32,[]);` to split your signal then two `for` loops can solve the problem. But I'm not sure if that is what you want. – Rashid Feb 20 '16 at 21:37
  • @nullgeppetto, it seems you want to see which frequencies happen in each 3sec window. You can consider time-frequency distributions and wavelet maybe. – Rashid Feb 20 '16 at 21:41
  • @nullgeppetto, see (https://en.wikipedia.org/wiki/Short-time_Fourier_transform) for short time Fourier transform, check out the graphs (time vs frequency). – Rashid Feb 20 '16 at 21:44
  • Thank you very much for all your help! – nullgeppetto Feb 20 '16 at 21:45
  • 1
    @nullgeppetto, you're welcome, thanks for accepting the answer. – Rashid Feb 20 '16 at 21:48