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

16 views (last 30 days)
Biljana basarin on 16 Mar 2013
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.

Wayne King on 19 Mar 2013
Edited: Wayne King on 19 Mar 2013
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 on 20 Mar 2013
thank you... I am new at all this Thank you for you help and time spent. Sorry for inconvenience.

Wayne King on 16 Mar 2013
Edited: Wayne King on 16 Mar 2013
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 on 17 Mar 2013
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.

Wayne King on 18 Mar 2013
Edited: Wayne King on 18 Mar 2013
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 on 18 Mar 2013
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.

### Categories

Find more on Spectral Estimation in Help Center and File Exchange

### Tags

No tags entered yet.

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!