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

16 views (last 30 days)

Show older comments

##### 0 Comments

### Accepted Answer

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

### More Answers (2)

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.

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

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!