pwelch method of power spectral density (psd) calculation

I have a question on the pwelch method for power spectral density computation in Matlab that is not clarified in the documentation for pwelch.

My usage : pxx = pwelch(x, window, noverlap, nfft)

1. If nfft is larger than the length of window, I believe MATLAB pads each windowed segment with zeros so that the FFT is taken on nfft points. Please confirm.
2. What happens when nfft is smaller than the window length? Is the window size truncated to match nfft?
3. pwelch appears to scale the DFT here appropriately with number of samples of the nfft (N) and frequency spacing as opposed to a regular DFT in Matlab that is not scaled, pl confirm.
4. Is there a way to output the number of averages performed when pwelch is called as above?

Thanks.

[Paul]

Hi Sundeep,

Have you checked the doc page pwelch? I'm pretty sure that the answers can be found there, certainly at least for items 1-3.

[Sundeep]

Hi Paul, in Plain sight, thanks for the response.

[dpb]

Item 4 can be inferred/calculated from the statement in the More About section "Welch’s method computes a modified periodogram for each segment and then averages these estimates to produce the estimate of the power spectral density."

The number is not returned as an available output from the function, no, it is dependent upon the values chosen for the window length, overlap and nfft that combine to determine the segment length which then determines how many segments there are going to be.

[Paul]

I don't think nfft would affect the number of segments.

Unfortunately, the doc page for pwelch doesn't provide much information about a) how the segment length is computed, and b) what it does if the segments don't fit "nicely" with original input signal.

As far as b) is concerned, I did a quick code check and it seems to me that pwelch ignores data at the end of the input that can't be used to form a full segment. In contrast, pspectrum states explicitly "pspectrum zero-pads the signal if the last segment extends beyond the signal endpoint. The function returns t, a vector of time instants corresponding to the centers of the segments." i.e., numel(t) is the number of segments.

[dpb]

From further reading, I agree that NFFT won't matter.

As noted in the help doc, PWELCH divides each column of X into overlapping sections of the same length as the window. If the length of the X vector is such that it cannot be divided exactly into an integer number of sections with the defined overlap, X is truncated to that length.

I didn't go digging deep enough to check whether if the window is 128 and the overlap is 50% whether the second section starts at index 64 (the overlap number) or 65 (the next index after the overlap) although I presume it would probably be the latter. Not accounting for that might end up with an "off by one" number of segments causing a truncation earlier in the latter case.

Answer 1

[dpb]

OK, I dove into the bowels...it's straightforward although I haven't worked out a simple algebraic relationship, for the example in the doc of a 512-pt seriies with a window of 132 points and a 50% overlap, the start and end points over which the segments are defined are

LminusOverlap = L-noverlap;
xStart = 1:LminusOverlap:k*LminusOverlap;
xEnd   = xStart+L-1;

K>> L 
L =
   132
K>> LminusOverlap
LminusOverlap =
    66
K>
K>> [xStart;xEnd]
ans =
     1    67   133   199   265   331
   132   198   264   330   396   462
K>>

As can be seen, there are k=6 segments another segment of 132 points would go past the 512 signal length.

An internal inaccessible routine parses the inputs and returns L and k; to compute the number of elements in k for any given signal length, window length and overlap, compute the start and end locations in a loop until the end matches or exceeds the signal length. If it exceeds, k is one less than the loop counter, if it turns out even, then they're the same.

ADDENDUM:

It is simpler to compute the value of k and the start, end points of the segments using colon addressing --

function [k,i1,i2]=computeSegments(N,L,f)
% return number of segments and start, end of segments for PWELCH
% Inputs: N - length of time series
%         L - window length
%         f - fraction of window to overlap
% Outputs:k - number of segments
%         i1- segment start (optional)
%         i2- segment end (optional)

  noverlap=fix(L*(1-f));              % number points of overlap
  i2=L:noverlap:N;                    % end of segments
  i1=i2-L+1;                          % start of segments
  k=numel(i2);                        % how many segments
end

[k,~,iend]=computeSegments(512,132,0.5) % call ask for only number and end points

k = 6

iend = 1×6

132 198 264 330 396 462

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>