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>

ADDENDUM SECOND:

Nota Bene the difference between determining the start positions between the first debug session output that is

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

which requires the internal parsing function (that I couldn't see the code for) to have predetermined k in order to set the end value in the colon expansion expression. Then, the end points were derived from those by.

xEnd   = xStart+L-1;

That is the internals of the MATLAB version.

In comparison, given that there is no way to return k from the code; one has to go at it somewhat differently; it is straightforward to find the end positions first directly from

i2=L:noverlap:N;

using the builtin behavior of colon to stop the expansion before or at N rather than the next case after. From having i2 (xEnd) one can then rearrange the above relationship to solve for xStart (i1) .

Why the developer went the other way is indeterminate; the latter above seems more direct, but either produces the same result.

In the above sample function, one could also add some additional logic to be able to pass either the fractional overlap as is done above or a fixed overlap number. Other error checking and enhancements are left for Student <grin>

[Sundeep]

The snippet that extracts the number of segments and their indices was helpful. Thank you for taking the time to explain the nuances of the two approaches.

To the original issue:

1. 'If nfft is greater than the segment length, the data is zero-padded' - Understood.
2. 'If nfft is less than the segment length, the segment is wrapped using datawrap to make the length equal to nfft'. So a circular folding of the data before computing the FFT, which can significantly change the spectral estimate compared to truncation or zero-padding.

So what may be the motivation for applying a data‑wrap instead of simply raising an error when nfft is smaller than the segment length?

Thanks.

[Paul]

Good question.

Also, interesting to me that datawrap is not documented in the sense that it doesn't have its own doc page.

[dpb]

Wrapping ensures continuity of the input signal at the ends -- otherwise the discontinuity will create spectral leakage adding energy across the spectrum. How much will depend on just what the discontinuity matnitude/shape is. Truncating has the effect of a rectangular window, at least partially negating the windowing function.

If you're concerned, you can ensure this case doesn't occur in use with proper inputs.

[Paul]

The undocumented datawrap does not ensure continutity at the ends of a segment, at least now how I see it.

Here's an 8-point segment

s = 1:8

s = 1×8

1 2 3 4 5 6 7 8

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

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

datawrapped to nfft = 5

sdw = datawrap(s,5)

ans = 1×5

7 9 11 4 5

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

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

What would the fft of sdw tell us about s? Maybe there's something more to this when considering multiple segments with overlap?

[Sundeep]

Agree on discontinuity - if  's' is a window function that tapers to zero at both ends, applying the data‑wrap operation forces non‑zero values at the boundaries; consequently, the product nfft .* datawrap(s) shall create discontinuities between segments, leading to leakage.

Perhaps - Maybe there's something more to this when considering multiple segments with overlap?

Thanks.

[dpb]

That's interesting; I confess to having jumped to a conclusion based on the function name rather than actually looking at it in detail. However, consider a "real" signal...

s=randn(8,1);
NFFT=5;
sdw=datawrap(s,NFFT);
plot(1:numel(s),s,'x-',1:numel(sdw),sdw,'+:')
legend('S','SDW')

Internally, the input vector is converted into a set of columns of NFFT length, padding as necessary and then those columns are added -- ...

[buffer(s,NFFT) sum(buffer(s,NFFT),2) sdw]

ans = 5×4

0.0616 0.3689 0.4305 0.4305 1.2941 0.0312 1.3253 1.3253 0.2204 0.4485 0.6688 0.6688 0.9568 0 0.9568 0.9568 0.3120 0 0.3120 0.3120

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

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

A=sum(ans);  % total amplitude
[sum(A(1:2)) A(end)]

ans = 1×2

3.6934 3.6934

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

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

What this does is to put the same amount of energy in the five points that was in the original eight. On reflection, that seems OK and about all one could do. The best course of action is to ensure don't do this.

[Paul]

The energy is not conserved because a^2 + b^2 ~= (a + b)^2

rng(100);
N = 10;
x = rand(1,N);
X = fft(x);
Ndw = 7;
Xdw = fft(datawrap(x,Ndw));
[sum(x.*x), sum(conj(X).*X)/N, sum(conj(Xdw).*Xdw)/Ndw]

ans = 1×3

2.7628 2.7628 4.2247

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

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

The datawrap also changes the apparent frequency content (see this answer for an explanation of what's happening here):

figure
plot((0:N-1)/N,abs(X),'-o',(0:Ndw-1)/Ndw,abs(Xdw),'-o'),grid

Maybe these effects aren't as pronounced with a tapered window?

I did the tiniest bit of searching to try to understand if this datawrap approach is justified. I didn't find anything (again, not expending much effort). But I did find two other packages (R and Octave) that simply ignore nfft if it is smaller than the window length.

[dpb]

Note I said the signal amplitude was conserved, not energy. <g>

One would have to ask Mathworks speciically why they chose to compute some approximation for the case instead of bailing; that's indeterminate otherwise. It can be avoided to not be an issue so it is somewhat of a moot issue other than simply National Enquirer "inquiring minds" kind of curiosity.

[Sundeep]

Maybe I am misinterpreting the term 'segment' here.

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

And the below text in italics is from pwelch doc.

nfft — Number of DFT points

max(256,2^nextpow2(length(window))) (default) | integer | []

Number of DFT points, specified as a positive integer. For a real-valued input signal, x, the PSD estimate, pxx has length (nfft/2 + 1) if nfft is even, and (nfft + 1)/2 if nfft is odd. For a complex-valued input signal,x, the PSD estimate always has length nfft. If nfft is specified as empty, the default nfft is used.

1. If nfft is greater than the segment length, the data is zero-padded.
2. If nfft is less than the segment length, the segment is wrapped using datawrap to make the length equal to nfft.

To me the word segment is the number of samples contained in the window vector. Am I misinterpreting the segment defintion here?

While both of you have both presented examples where nfft > window length and are applying the datawrap to the segment of the signal of original length nfft to get to window length, right?

Thanks.

[dpb]

The segment length is the window length, yes.

You asked specifically about the case in which NFFT is less than the segment length, not greater than; that's the only time in which datawrap is used to shorten the requested length down to NFFT; that is what the examples above illustrate; @Paul used the integers thinking of them as indices(?) whereas I just created a dummy time series since what is done is to combine the whole segment into one NNFT-long short segment that is the sum of the segment entities modulo NFFT and then the DFT is applied to that composite vector.

But, you can avoid this ever happening if you always ensure NFFT >= numel(window) which goes back to your earlier commen/queryt of why didn't Mathworks decide to just barf/error in the short NFFT case? -- which is what Paul noted that at least a couple of other implementations do. I commented that is indeterminate unless submit an official support request (or a Mathworks staffer who was around when the decision was made happens by and comments). I think while perhaps somewhat interesting, it is of no real consequence since you can prevent it from ever being an issue.

[Sundeep]

ok, yes.

Please help clarify one more time (for case 2) -

a) The datawrap is to be applied on the larger window vector to shorten to nfft (to match signal length) and not to the signal itself.

b) Your example called out a "signal", which is why the datawrap step seemed puzzling at first. As noted applying a datawrap to a typical window, the endpoints of the window may no longer remain zero.

And yes, I agree this is of not great consequence since one can avoid the case where nfft < length of the window segment (... can be dangerous).

Thanks again.

[Paul]

@Sundeep

"While both of you have both presented examples where nfft > window length and are applying the datawrap to the segment of the signal of original length nfft to get to window length, right? "

In my example in this comment, the segment s = 1:8 was intended to represent an 8-point segment, i.e., window length of 8, of a signal that has length greater than or equal to 8. I then applied datawrap with nfft = 5, so nfft was less than the window length and therefore less than the original length of the signal.

@dpb

This comment states: "What this does is to put the same amount of energy in the five points that was in the original eight."

I don't see anywere where this thread is talking about signal amplitude being conserved (I'm not even sure what that means).

[Paul]

"a) The datawrap is to be applied on the larger window vector to shorten to nfft (to match signal length) and not to the signal itself."

As I understand it, datawrap is applied to the product of the window function and the segment of x to shorten the product to have length nfft. I verified this by stepping through the code.

However, the documentation states: " If nfft is less than the segment length, the segment is wrapped using datawrap to make the length equal to nfft." I think it would be more accurate to say "If ..., the windowed segment is wrapped ...." because at the top of the doc page the term "segment" means an extracted portion of x before applying the window.

I'm not sure what you mean by "to match signal length." The "signal" in context of pwelch is the entirety of x and nothing is being shortened to match the length of x.

[Sundeep]

I did not oginally think that the datawrap was perfomed on the "windowed segment" - i.e. product of the window function and an extracted segment of x of length of the window vector" and subsequently shorten the product to a length of nfft.

My thoughts were a nfft portion of the signal is extracted - that is a signal 'segment' taken, datawrap performed on the window vector (only) to reduce to nfft size, product taken, and a DFT performed. This description felt a bit convoluted, and the wording in the documentation for the nfft section of the pwelch help may have contributed to this.

Thanks for stepping thru the code and clarifying. Thanks dpb for the clarifications as well.

[Paul]

If you think the doc page can be improved, click on a rating star at the bottom of the page and then a dialog box will pop up where you can type in your feedback.

We can illustrate how it works.

Start with short signal

rng(100);
N = 10;
x = rand(N,1);

Define the window length for the simplest case where the entire signal forms one segment and the FFT length is shorter than the segment length.

nfft = 7;win = N;

Call pwelch

% this usage of pwelch is technically undocumented because noverlap is supposed to be
% a positive integer
[P,f] = pwelch(x,hamming(N),0,nfft,1,'centered');

By-hand calculations:

Define the window

h = hamming(N,'symmetric');

Multiply the segment, which is the full signal in this example, by the window and datawrap the product

sw = datawrap(x.*h,nfft);

Compute the nfft-point DFT of the datawrapped, windowed segment

Sw = fftshift(fft(sw));

Compute the PSD

PSw = Sw.*conj(Sw)/norm(h)^2;

Compare to the result from pwelch

[f,P,PSw]

ans = 7×3

-0.4286 0.0332 0.0332 -0.2857 0.3107 0.3107 -0.1429 0.0052 0.0052 0 1.1529 1.1529 0.1429 0.0052 0.0052 0.2857 0.3107 0.3107 0.4286 0.0332 0.0332

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

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

[Sundeep]

@Paul - Thank you for sharing the code and show datawrap is performed on the 'windowed segment' to get it to a nfft size.

I am leaving for an event this morning and want to do the same for the case 1 -

1. If nfft is greater than the segment length, the data is zero-padded.

It appears trivial but want to confirm pwelch is doing what one is expecting here.

[dpb]

I don't see anywere where this thread is talking about signal amplitude being conserved (I'm not even sure what that means)"

It falls out of the datawrap function algorithm and is what was illustrated in the example -- the total amplitude of the wrapped signal is the same as that of the original segment -- that's the result of (and one presumes the reason for having doine) the summation across columns. If you want to call it somnething else, so be it...

[Paul]

@dpb

I must not have seen your edit to this comment. It sounds like you're using the term "total amplitude" to mean the sum of the elements, correct?

rng(100)
s = rand(8,1);            % original segment
nfft = 5;
sdw = datawrap(s,nfft);   % wrapped segement
[sum(s),sum(sdw)]

ans = 1×2

3.7140 3.7140

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

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

I'm not familiar with that definition of "total amplitude." The only other thing I would call it is the "sum of the elements."

I'm still not sure what you meant by "What this does is to put the same amount of energy in the five points that was in the original eight." (emphasis added). Were you using "energy" synonymously with "total amplitude" synonymously with "sum of the elements"?

[Paul]

@Sundeep

My understanding of case 1 is as follows.

rng(100);
N = 10;
x = rand(N,1);
nfft = 13;win = N;
% this usage of pwelch is technically undocumented because noverlap is supposed to be
% a positive integer
[P,f] = pwelch(x,hamming(N),0,nfft,1,'centered');
h = hamming(N,'symmetric');
sw = x.*h;
Sw = fftshift(fft(sw,nfft));  % zero pads sw and then takes the FFT of the padded result
PSw = Sw.*conj(Sw)/norm(h)^2;
[f,P,PSw]

ans = 13×3

-0.4615 0.0474 0.0474 -0.3846 0.0688 0.0688 -0.3077 0.2795 0.2795 -0.2308 0.2629 0.2629 -0.1538 0.0238 0.0238 -0.0769 0.3626 0.3626 0 1.1529 1.1529 0.0769 0.3626 0.3626 0.1538 0.0238 0.0238 0.2308 0.2629 0.2629 0.3077 0.2795 0.2795 0.3846 0.0688 0.0688 0.4615 0.0474 0.0474

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

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

You may want to try recreating pwelch results for more than one segement, with or without overlap.

[dpb]

Well, if one adds up the amplitudes of the signals, I'd think "total amplitude": is a pretty good name for it... You see any other explanation for having done so? I suppose they could have computed the sum of squares and renormalized to conserve that quantity instead.

If I wrote "energy" somewhere, sorry, it was an inadvertent use of the word.

[Paul]

Amplitude is typically considered to be a positive quantity. In your example, all of the elements of the segment were positive so the sum of the amplitudes of the elements and the sum of the elements themselves are the same. But what if the segment has both positive and negative elements? datwrap will preserve the sum of the elements but not the sum of the amplitudes. So "total sum" might be a better term in the context of datawrap than "total amplitude", which would more naturally be defined as sum(abs(x)), which some might call the absolute sum of x (as related to absolute summability, which is a very important concept in control theory and signal processing).

[dpb]

Whatever you say....sorry I let myself get dragged into this...

Answer 2

[Paul]

As to why pwelch applies datawrap for the case nfft < numel(win) ...

I did a bit more searching to try to find the motivation for the pwelch implementation and found nothing other than DATAWRAP - MATLAB Answers - MATLAB Central which is essentially asking the same question as @Sundeep.

After reading that thread, I have a hypothesis as to what's going on.

Let sw be an eight-point windowed segment.

sw = rand(1,8);

One period of its DTFT is

[dtftsw,w] = freqz(sw,1,4096,'whole');

Plotting just the amplitude for simplicity

hf1 = figure;
plot(w,abs(dtftsw),'DisplayName','DTFT'),grid

As must be the case, the 8-point DFT of sw are samples of its DTFT

dftsw = fft(sw);
hold on
stem((0:7)/8*2*pi,abs(dftsw),'DisplayName','DFT-8')

If we want 14 samples of the DTFT (i.e., nfft > numel(win)), we zero pad and then take the DFT, which again are samples of he DTFT.

dftsw = fft(sw,14);
stem((0:13)/14*2*pi,abs(dftsw),'DisplayName','DFT-14')

But suppose we want fewer than eight samples of the DTFT of sw, perhaps five. How would we get those samples of the DTFT? We could just use freqz, specifying the five frequencies of interest. Or, we could use fft and take advanage of the fact that undersampling in the frequency domain leads to aliasing in the time time domain. It is that aliasing in the time domain that is implemented by datawrap.

dftsw = fft(datawrap(sw,5)); % five point DFT after datawarapping
stem((0:4)/5*2*pi,abs(dftsw),'DisplayName','wDFT-5');
legend('Location','North')

In summary, datawrap is applied to the windowed segment when nfft < numel(win) so that the DFT of the output of datawrap are frequency domain samples of the DTFT of the windowed segment that is the input to datawrap.

Keep in mind that the aliasing in the time domain could be problematic. Consider a corner case where the windowed segment has a certain symmetry

rng(100);
z = rand(1,4);
sw = [z,-z]  % windowed segment of x

sw = 1×8

0.5434 0.2784 0.4245 0.8448 -0.5434 -0.2784 -0.4245 -0.8448

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

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

The 4-point output of datawrap is all zero

wrappedsw = datawrap(sw,4)

wrappedsw = 1×4

0 0 0 0

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

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

which is consistent with the fact that every other frequency domain sample of the DTFT of sw is zero

abs(fft(sw))

ans = 1×8

0 2.4541 0 2.0274 0 2.0274 0 2.4541

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

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

[Sundeep]

It’s great to learn that the aliasing we observe in the time domain is what a data‑wrap operation may be implementing here. Thank you for digging into this detail with an example and formulating a hypothesis.

May be a rigorous proof can be constructed by examining the inverse Discrete Fourier Transform of a frequency‑domain signal that has been undersampled. I am not asking for one - just thinking. Thanks.

[Paul]

There is no question, in my mind, that datawrap implements the aliasing in the time domain that results from N-periodic summation of a finite duration signal of length L with N < L.

Undersampling in the frequency domain leading to aliasing in the time domain would be covered in any useful textbook on discrete-time signal processing. No need to construct a proof, except as a learning exercise, which would be useful.

My hypothesis wasn't about the relationship between undersampling and aliasing, which is a fact. Rather, I was hypothesizing about the motivation of the developers of pwelch for taking this approach. I'm quite certain their approach is motivated by what I've shown, i.e., use datawrap if nfft < numel(win) to get the correct frequency domain samples needed for the spectral estimation, but I didn't want to say that w/ absolute certainty because the doc is silent on their reasoning.

[Sundeep]

ok, understand your point, is good to know.

I will let MATLAB know —via feedback—about the segmented‑window issue (the fact that the product of the window function and an extracted segment of x whose length matches the window vector is performed first, followed by the data‑wrap). It is ambiguous and not clear from their current documentation.

I will go ahead and add a question about the motivation for using the datawrap: ask why such an internal implementation for case 2 ( nfft < window‑segment length ) is not described - i.e. arrive at the correct spectal estimates albeit under sampled in the frequency domain (or what additional motivation).

Thanks so much for showing the work.