Question about Matlab's implementation of Welch PSD when window size is larger than NFFT

Hi, when I looked into Matlab's PWELCH function, I noticed something that I don't understand.

% Suppose we have some data x which has 1000 samples.
x = randn(1000,1);

% And we use a window of size 200.
nWin = 200;
window = hamming(nWin);

% Now we also specify nFFT to 128
nFFT = 128;

% what Matlab's PWELCH does, is like the following:

% grab one frame
oneFrame = x(1:nWin);
% apply the window;
oneFrameWindowed = window.*oneFrame;

% datawrap it. WHY SUM but not MEAN?
oneFrameWindowedWrapped = sum(buffer(x,nFFT),2);

% compute periodogram for one frame.... etc.
psdOneFrame = fft(oneFrameWindowedWrapped,nFFT);



What confuses me is the datawrap uses SUM but not MEAN. Furthermore, why after all do we need a datawrap stage?  Shouldn't we use fft(oneFrameWindowed,nFFT) to get a length nFFT periodogram directly?

I appreciate your kindly help.

The FFT operator is periodic with period nFFT anyway, so using the buffer command, the summing the 2nd dimension, saves operation and time. 

Consider:
x = randn(1024,1);
xFFT = fft(x,512);

Is equivalent to:
xFFT = fft(sum(buffer(x,512),2),512);


On Monday, November 25, 2013 7:23:01 PM UTC-5, mr.siy...@gmail.com wrote:
> Hi, when I looked into Matlab's PWELCH function, I noticed something that I don't understand.
> 
> 
> 
> % Suppose we have some data x which has 1000 samples.
> 
> x = randn(1000,1);
> 
> 
> 
> % And we use a window of size 200.
> 
> nWin = 200;
> 
> window = hamming(nWin);
> 
> 
> 
> % Now we also specify nFFT to 128
> 
> nFFT = 128;
> 
> 
> 
> % what Matlab's PWELCH does, is like the following:
> 
> 
> 
> % grab one frame
> 
> oneFrame = x(1:nWin);
> 
> % apply the window;
> 
> oneFrameWindowed = window.*oneFrame;
> 
> 
> 
> % datawrap it. WHY SUM but not MEAN?
> 
> oneFrameWindowedWrapped = sum(buffer(x,nFFT),2);
> 
> 
> 
> % compute periodogram for one frame.... etc.
> 
> psdOneFrame = fft(oneFrameWindowedWrapped,nFFT);
> 
> 
> 
> 
> 
> 
> 
> What confuses me is the datawrap uses SUM but not MEAN. Furthermore, why after all do we need a datawrap stage?  Shouldn't we use fft(oneFrameWindowed,nFFT) to get a length nFFT periodogram directly?
> 
> 
> 
> I appreciate your kindly help.

Thank you for the reply, but apparantly the 2 ffts listed in your example are not equivalent?

x = [2;4;3;5];

f1 = fft(x,2)
f1 = 6, -2


f2 = fft(sum(buffer(x,2),2),2)
f2 = 14, -4

Maybe of interest: 

pwelch() is useful when the data stream is too long to be processed "in one
shot". With the computing power of a modern PC, these cases tend to get
quite rare. 

Instead, it's typically used to reduce the number of bins in the output. 
An alternative is to filter and decimate the power spectrum estimate.

A piece of code that does this (which I've been using for a very long time)
can be found here: 
http://www.dsprelated.com/showcode/239.php

If the signal is not cyclic, use any window of your choice on the input
signal, as with conventional FFT. 

Used correctly, the results are very close to readings from a spectrum
analyzer instrument (which has some vendor dependency, i.e. the exact shape
of the RBW filter).
	 

_____________________________		
Posted through www.DSPRelated.com

Thanks for the info, but this does not answer my question.



On Wednesday, November 27, 2013 11:21:13 AM UTC-8, mnentwig wrote:
> Maybe of interest: 
> 
> 
> 
> pwelch() is useful when the data stream is too long to be processed "in one
> 
> shot". With the computing power of a modern PC, these cases tend to get
> 
> quite rare. 
> 
> 
> 
> Instead, it's typically used to reduce the number of bins in the output. 
> 
> An alternative is to filter and decimate the power spectrum estimate.
> 
> 
> 
> A piece of code that does this (which I've been using for a very long time)
> 
> can be found here: 
> 
> http://www.dsprelated.com/showcode/239.php
> 
> 
> 
> If the signal is not cyclic, use any window of your choice on the input
> 
> signal, as with conventional FFT. 
> 
> 
> 
> Used correctly, the results are very close to readings from a spectrum
> 
> analyzer instrument (which has some vendor dependency, i.e. the exact shape
> 
> of the RBW filter).
> 
> 	 
> 
> 
> 
> _____________________________		
> 
> Posted through www.DSPRelated.com

no, it doesn't. That's why the first line reads: "maybe of interest" :-)	 

_____________________________		
Posted through www.DSPRelated.com

On Monday, November 25, 2013 4:23:01 PM UTC-8, mr.siy...@gmail.com wrote:
> Hi, when I looked into Matlab's PWELCH function, I noticed something that I don't understand.
> ...

I don't have access to Matlab R2013b to test but, the syntax is here:
http://www.mathworks.com/help/signal/ref/pwelch.html

Three comments on the following parts of your description:

> % datawrap it. WHY SUM but not MEAN?
> oneFrameWindowedWrapped = sum(buffer(x,nFFT),2);

1) The second argument to buffer should be nWIN not nFFT

> % compute periodogram for one frame.... etc.
> psdOneFrame = fft(oneFrameWindowedWrapped,nFFT);

2) The fft should be before the sum or mean for Welch's method

3) you need to calculate the psd after the fft and before the sum or mean for this to calculate Welch's method

> What confuses me is the datawrap uses SUM but not MEAN. Furthermore, why after all do we need a datawrap stage?  Shouldn't we use fft(oneFrameWindowed,nFFT) to get a length nFFT periodogram directly?

Given the errors in your description so far, are you sure that pwelch doesn't use mean as the documentation says?

> ...

Dale B. Dalrymple

Thanks for the reply, but I did no mistake in my description, if you have a chance please open Matlab and you will see.  

Anyway:
1) the second argument IS nFFT, this stage occurs after windowing is done: you have nWin and nFFT and you need to apply both, right?

2) That is what confuses me, the FFT occurs AFTER the SUM in time domain. 

I know this for SURE because I implemented it in c# and have to use Matlab's peculiarity in order to get identical results.

On Monday, December 2, 2013 3:13:11 PM UTC-8, dbd wrote:
> On Monday, November 25, 2013 4:23:01 PM UTC-8, mr.siy...@gmail.com wrote:
> 
> > Hi, when I looked into Matlab's PWELCH function, I noticed something that I don't understand.
> 
> > ...
> 
> 
> 
> I don't have access to Matlab R2013b to test but, the syntax is here:
> 
> http://www.mathworks.com/help/signal/ref/pwelch.html
> 
> 
> 
> Three comments on the following parts of your description:
> 
> 
> 
> > % datawrap it. WHY SUM but not MEAN?
> 
> > oneFrameWindowedWrapped = sum(buffer(x,nFFT),2);
> 
> 
> 
> 1) The second argument to buffer should be nWIN not nFFT
> 
> 
> 
> > % compute periodogram for one frame.... etc.
> 
> > psdOneFrame = fft(oneFrameWindowedWrapped,nFFT);
> 
> 
> 
> 2) The fft should be before the sum or mean for Welch's method
> 
> 
> 
> 3) you need to calculate the psd after the fft and before the sum or mean for this to calculate Welch's method
> 
> 
> 
> > What confuses me is the datawrap uses SUM but not MEAN. Furthermore, why after all do we need a datawrap stage?  Shouldn't we use fft(oneFrameWindowed,nFFT) to get a length nFFT periodogram directly?
> 
> 
> 
> Given the errors in your description so far, are you sure that pwelch doesn't use mean as the documentation says?
> 
> 
> 
> > ...
> 
> 
> 
> Dale B. Dalrymple

On Tuesday, December 3, 2013 9:59:35 AM UTC-8, Siyi wrote:
> Thanks for the reply, but I did no mistake in my description, if you have a chance please open Matlab and you will see.  
> 

I would be useful to know what test signal and argument values you have used to verify your description as equivalent to pwelch.

> Anyway:
> 1) the second argument IS nFFT, this stage occurs after windowing is done: you have nWin and nFFT and you need to apply both, right?
> 
> 2) That is what confuses me, the FFT occurs AFTER the SUM in time domain. 
> 
> I know this for SURE because I implemented it in c# and have to use Matlab's peculiarity in order to get identical results.
> 

My comments were an attempt to make your description consistent with Welch's method. The snippets you have provided do not describe Welch's method. You can find Welch's paper here:
http://www.utdallas.edu/~cpb021000/EE%204361/Great%20DSP%20Papers/Welchs%20Periodogram.pdf

Matlab's R2013b documentation for pwelch:
http://www.mathworks.com/help/signal/ref/pwelch.html
describes processing consistent with Welch's paper.

If you can provide a more complete description of your interpretation of pwelch we could discuss what the process you are describing actually does. That might allow us to determine what is going on.

FFT after the sum does not appear in Welch's method, but is used in operations sometimes described in the literature as "time-aliasing".

Dale B. Dalrymple

% Suppose we have some data x which has 1000 samples. 
x = randn(1000,1); 

% And we use a window of size 200. 
nWin = 200; 
window = hamming(nWin); 

% Now we also specify nFFT to 128 
nFFT = 128; 

nSample = numel(x);

%% Lets re-create what Matlab's PWELCH does:

nOverlap = 0;
nFrame = nSample/nWin; 

psd0 = 0;
for k = 1:nFrame

    % grab one frame 
    oneFrame = x((1:nWin)+(k-1)*nWin);

    % apply the window; 
    oneFrameWindowed = window.*oneFrame; 

    % datawrap it.
    oneFrameWindowedWrapped = sum(buffer(oneFrameWindowed,nFFT),2); 

    % compute periodogram for one frame.... etc. 
    psdOneFrame = fft(oneFrameWindowedWrapped,nFFT);
    
    psdOneFrame = abs(psdOneFrame).^2;
    
    psd0 = psd0+psdOneFrame;
end

psd0 = psd0(1:(nFFT/2+1));

% normalization, etc.
psd0 = psd0/sum(window.^2)/pi/nFrame;

%% Matlab's PWELCH;

psd1 = pwelch(x,nWin,nOverlap,nFFT);

figure; plot([psd0,psd1]); % identical except for dc components;





On Tuesday, December 3, 2013 12:26:21 PM UTC-8, dbd wrote:
> On Tuesday, December 3, 2013 9:59:35 AM UTC-8, Siyi wrote:
> 
> > Thanks for the reply, but I did no mistake in my description, if you have a chance please open Matlab and you will see.  
> 
> > 
> 
> 
> 
> I would be useful to know what test signal and argument values you have used to verify your description as equivalent to pwelch.
> 
> 
> 
> > Anyway:
> 
> > 1) the second argument IS nFFT, this stage occurs after windowing is done: you have nWin and nFFT and you need to apply both, right?
> 
> > 
> 
> > 2) That is what confuses me, the FFT occurs AFTER the SUM in time domain. 
> 
> > 
> 
> > I know this for SURE because I implemented it in c# and have to use Matlab's peculiarity in order to get identical results.
> 
> > 
> 
> 
> 
> My comments were an attempt to make your description consistent with Welch's method. The snippets you have provided do not describe Welch's method. You can find Welch's paper here:
> 
> http://www.utdallas.edu/~cpb021000/EE%204361/Great%20DSP%20Papers/Welchs%20Periodogram.pdf
> 
> 
> 
> Matlab's R2013b documentation for pwelch:
> 
> http://www.mathworks.com/help/signal/ref/pwelch.html
> 
> describes processing consistent with Welch's paper.
> 
> 
> 
> If you can provide a more complete description of your interpretation of pwelch we could discuss what the process you are describing actually does. That might allow us to determine what is going on.
> 
> 
> 
> FFT after the sum does not appear in Welch's method, but is used in operations sometimes described in the literature as "time-aliasing".
> 
> 
> 
> Dale B. Dalrymple