I need to calculate time-evolving power spectral density using Matlab periodogram function

I need to calculate time-evolving power spectral density using Matlab periodogram function based on Welch theory to estimate the PSD of a moving 400-kyr boxcar filter with an overlap of 85%. In all the figures, values have to be plotted at the center of each 400-kyr window over which they are calculated.

Answer 1

[Wayne King]

You cannot have a PSD estimate on the same time scale as the original signal. That would mean that you have a PSD estimate based on 1 sample. You must sum over some interval to produce a PSD estimate.

I tried to help you by telling you that you need to create a time vector that essentially takes the midpoints of each time interval. Remember, you time jump between intervals is really the length of your window minus the number of overlapped samples, so it's 30 samples with a sampling frequency of 1000 Hz -- 0.030 sec.

I also told you that buffer() will prepend zeros and append zeros so you need to take those into account.

Fs = 1000;
t = 0:0.001:4-0.001;
x = cos(2*pi*10*t)+randn(size(t));
winsize = 200;
numoverlap = round(0.85*winsize);
win = hamming(200);
X = buffer(x,200,numoverlap);
for nn = 1:size(X,2)
[Pxx(:,nn),F] = pwelch(X(:,nn),win,length(win)/2,length(win),Fs);
end

 % create a time vector
 idxbegin = find(X(:,1) == 0);
 numpresteps = length(idxbegin);
 idxend = find(X(:,end) == 0);
 numpoststeps = length(idxend);
 tbegin = -(numpresteps*dt)/2;
 tend = t(end)+((numpoststeps*dt))/2;
 tvec = linspace(tbegin,tend,size(Pxx,2));

 surf(tvec,F,10*log10(abs(Pxx)),'EdgeColor','none');   
 axis xy; axis tight; colormap(jet); view(0,90);
 xlabel('Time (sec)');
 ylabel('Frequency (Hz)');

[Biljana basarin]

thank you... I am new at all this Thank you for you help and time spent. Sorry for inconvenience.

Answer 2

[Wayne King]

If you want to use Welch's method in a time-evolving manner, use buffer() to segment the signal with overlap and obtain Welch estimates on those overlapped segments.

I would caution you against using a boxcar filter, here I'll give you an example with a Hamming window. You can substitute your boxcar filter as needed.

Further, you have not been clear about whether the overlap of 85% applies to both the time-evolving PSD or the overlap in the Welch's estimate, I'll use 50% for the latter.

Fs = 1000;
t = 0:0.001:4-0.001;
x = cos(2*pi*10*t)+randn(size(t));
winsize = 200;
numoverlap = round(0.85*winsize);
win = hamming(200);
X = buffer(x,200,numoverlap);
for nn = 1:size(X,2)
  [Pxx(:,nn),F] = pwelch(X(:,nn),win,length(win)/2,length(win),Fs);
end

The columns of Pxx give you the time-varying Welch PSD estimates. You may want to avoid using the last column of Pxx because that is computed on the last column of X, which may contain a lot of zeros.

[Biljana basarin]

I tried to apply what you suggested, I am still not sure if I got the results I need as I dont know how to plot the results. The values have to be plotted at the center of each window over which they are calculated. Could you please explain how to do this? Thank you very much.

Answer 3

[Wayne King]

Just use surf(), you can easily work out a "meaningful" time vector but you have to look at the why buffer() has prepended and appended zeros to get the matrix right.

 Fs = 1000;
 t = 0:0.001:4-0.001;
 x = cos(2*pi*10*t)+randn(size(t));
 winsize = 200;
 numoverlap = round(0.85*winsize);
 win = hamming(200);
 X = buffer(x,200,numoverlap);
 for nn = 1:size(X,2)
  [Pxx(:,nn),F] = pwelch(X(:,nn),win,length(win)/2,length(win),Fs);
 end
 surf(1:size(Pxx,2),F,10*log10(abs(Pxx)),'EdgeColor','none');   
 axis xy; axis tight; colormap(jet); view(0,90);
 xlabel('Time');
 ylabel('Frequency (Hz)');

[Biljana basarin]

It still doesn't give me the result needed. Actually I need to show the variations of one signal component on time scale. I first, filtered the data at certain frequency, then I apply moving boxcar filter with the overlap of 85% and then calculate time evolving PSD using Matlab periodogram based on Welch theory. I need to plot PSD on the same time scale as original signal. After all of this I am now completely confused. Please help.