Using pwelch to set signals: some questions (Matlab)
I would like to use pwelch for a set of signals and I have some questions.
First, let's say we have 32 (EEG) signals of 30 seconds duration. The sampling rate is fs=256
samples / sec, and therefore each signal has a length of 7680. I would use pwelch
to estimate the power spectral density (PSD) of these signals.
Question 1:
Based on the pwelch
documentation ,
pxx = pwelch (x) returns an estimate of the power spectral density (PSD), pxx, of the input signal x, found using an average estimate over a Welch segmented segment. When x is a vector, it is treated as one channel. When x is a matrix, the PSD is calculated independently for each column and stored in the corresponding pxx column.
However, if called pwelch
like this
% ch_signals: 7680x32; one channel signal per each column
[pxx,f] = pwelch(ch_signals);
the resultant pxx
is sized 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 pxx column.
Question 2:
Let's say that I am overcoming this problem and I am calculating the PSD of each signal myself (applying pwelch
to each column ch_signals
), I would like to know what is the best way to do this. Granted that the signal is a 30 second signal in time with a sampling rate of fs=256
what should I call it pwelch
(with what arguments?) So does PSD make sense?
Question 3: If I need to split each of my 32 signals into windows and apply a batch to each 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 by 3 seconds with an overlap of 2 seconds. How can I call pwelch
for each of these windows?
source to share
Here is an example, just like your case,
The results show that the algorithm indicates the frequencies of the signals to the right.
Each column of the matrix y
is sinusoidal to test how it works.
Windows are 3 seconds with 2 seconds overlap,
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})');
source to share