quetsions about matlab pWelch implementation

Hi All,

I have two questions about matlab's pWelch implementation. I apologize if my questions have been asked elsewhere, but I have looked long and hard for an answer to my questions and I haven't stumbled upon them yet. So I decided to start a new question.

Some background: I was using matlab to for prototyping. Now I want to deploy my code in C# so I have to translate everything from matlab. Before I take that step, I thought I would first break down matlab's pWelch into its "sub-steps", and make sure I get the same result as the original function. I'm very close, but not getting identical results. According to the documentation for pWelch, the folks at mathworks follow these steps:

1.The input signal vector x is divided into k overlapping segments according to window and noverlap (or their default values). If the window size is larger than the number of FFT points (NFFT) , the signal is divided into NFFT–length segments and then, the last segment is padded with zeros.

2. The specified (or default) window is applied to each segment of x. (No preprocessing is done before applying the window to each segment.)

3. An nfft-point FFT is applied to the windowed data.

4. The (modified) periodogram of each windowed segment is computed.

5. The set of modified periodograms is averaged to form the spectrum estimate S(ejω).

6. The resulting spectrum estimate is scaled to compute the power spectral density as S(ejω)/F, where F is -2π when you do not supply the sampling frequency -fs when you supply the sampling frequency

So I just followed each of their steps, one by one, to reproduce the results. However, I'm off by some "normalization factor". My results are all multiplied by this factor. And another difference: the 0Hz, 1Hz, and 256Hz components are wrong. Everything else is correct. This is really baffling.

I digged a little deeper and found this article: <http://www.mathworks.com/support/solutions/en/data/1-17M0Q/index.html?product=SG&solution=1-17M0Q>

I don't really grasp the the whole "biased vs unbiased" concept, but I tried multiplying by norm(w)^2/sum(w)^2 and it still didn't fix my problem. I noticed they referenced a book, page 419: http://i.imgur.com/7mvak.png

So I took a look at that book and I see that they did have a different normalization factor. I tried that too and still no luck.

Below is the code I wrote to compare. My data is sampled at 512Hz. I am doing a pWelch of window length 512 points, and FFT length 512 points, with 50% overlap. Using the standard hamming window. I used Handel as the fake data, just so it's easier for you guys to run the code.

    load handel;
    data = y(1:2048)';

    %make a hamming window
    w = hamming(512);

    %preallocate the space for the psd results
    mypsd = zeros(512,7);

    %step 1: loop through the data, 512 points at a time, with 256 points overlap
    for i = 0:6
    %step 2: apply a hamming window
    temp = data(1+256*i:512+i*256)'.*w; 

    %step 3: calculate FFT and take just the first half
    temp = fft(temp);

    %step 4: calculate the "periodogram" by taking the absolute value squared
    temp = abs(temp).^2;   

    %save the results in the storage variable
    mypsd(:,i+1) = temp;
    end

    %step 5: take the average of all the periodograms
    mypsd_v1 = mean(mypsd,2);

    %step 6: divide by sampling rate (512Hz for my data)
    mypsd_v1 = mypsd_v1/512;

    %throw away the 2nd half of mypsd
    mypsd_v1 = mypsd_v1(1:257);

    %the matlab equivalent, using the pWelch function
    [p,f]=pwelch(data,512,.5,512,512);

    %now calculate the factor by which to multiply the data, so my approach and
    %the matlab approach match up
    norm_factor = mypsd_v1(115)/p(115);
    mypsd_v1 = mypsd_v1/norm_factor;

    %plot the results
    plot(0:length(p)-1,log(p),'b'), hold on, plot(0:length(p)-1,log(mypsd_v1),'r');
    title(sprintf('normalization factor: %f',norm_factor));
    xlabel('Frequency (Hz)');
    ylabel('Log Power');

    %version 2: apply the U factor. L = window length and K = # of sections
    U = (1/512)*sum(w.^2);
    K = 7;
    L = 512;
    mypsd_v2 = (1/(K*L*U))*sum(mypsd,2);
    mypsd_v2 = mypsd_v2(1:257);

    norm_factor = mypsd_v2(115)/p(115);
    mypsd_v2 = mypsd_v2/norm_factor;

    %plot the results
    figure, plot(0:length(p)-1,log(p),'b'), hold on, plot(0:length(p)-1,log(mypsd_v2),'r');
    title(sprintf('normalization factor: %f',norm_factor));
    xlabel('Frequency (Hz)');
    ylabel('Log Power');

Any help/pointers would be greatly appreciated!

Thanks, Neraj

Answer 1

[Moses]

Hi Neraj,

Did you ever get your code to work? You sub-steps are correct but it seems like some parts of your code are off, in particular the normalization factor. I made some changes to your code so that it now displays the same result.

For reference on the why this normalization factor was used, please refer equation #24 in this document. It is a very helpful and detailed paper on spectral density estimation and windows. http://holometer.fnal.gov/GH_FFT.pdf

The 2 comes from ignoring the redundant negative frequencies. Dividing by fs comes from wanting the power spectral density expressed in V^2/Hz. The sum of of the square of the window comes from the fact that the power is the two spectrums multiplied to each other (Pxx = Px ⋅ Px*, where ⋅ represents pointwise product and ∗ represents the complex conjugate). Hope this helps future viewers of this answer as well.

 load handel;
 data = y(1:2048)';

 % make a hamming window
 fs = 512;
 w = hamming(512);

 % preallocate the space for the psd results
 mypsd = zeros(512, 7);

 % step 1: loop through the data, 512 points at a time, with 256 points overlap
 for i = 0:6
    % step 2: apply a hamming window
    temp = data(1+256*i:512+i*256)'.*w;

    % step 3: calculate FFT and take just the first half
    temp = fft(temp);

    % step 4: calculate the "periodogram" by taking the absolute value squared
    temp = abs(temp).^2;

    % save the results in the storage variable
    mypsd(:, i+1) = temp;
 end

 % step 5: take the average of all the periodograms
 mypsd_v1 = mean(mypsd,2);

 % throw away the 2nd half of mypsd
 mypsd_v1 = mypsd_v1(1:257);

 % normalizing factor
 mypsd_v1 = mypsd_v1/(fs*sum(w.^2));

 % ignore the DC and Nyquist value
 mypsd_v1(2:end-1) = mypsd_v1(2:end-1) * 2;

 % the matlab equivalent, using the pWelch function
 [pxx, f] = pwelch(data, w, [], 512, fs);

 figure;
 subplot(2, 1, 1);
 plot(f, mypsd_v1);

 subplot(2, 1, 2);
 plot(f, pxx);

[Simon Wiede]

Hi! Is 6 in "for i = 0:6" the number of fft's?

[Walter Roberson]

Four full windows of 512 makes the length 2048. Add three windows at the half-window offsets, total 7. The 0 based index of the last of those is 6.

Answer 2

[Invizible Soul]

Hello,

Why in the following line of the code

xper = xper./(2*pi);

there is a division by 2*pi ?

cordially.

Answer 3

[Wayne King]

Hi Neraj, Your question and code are quite long, so I'll just give you a simple example using two windows with no overlap on a random vector. I'll set the random number generator first to its default so I can give you the final result, but that won't matter much. You'll see that the below agrees with the MATLAB pwelch().

rng default;
x = randn(1024,1); 
Pxx = pwelch(x,512,0,512);

xper = zeros(512,2);
for k = 0:1
x1 = x(1+k*512:(k+1)*512);
xw = x1.*hamming(512);
xper(:,k+1) = abs(fft(xw)).^2./norm(hamming(512),2)^2;
end
xper = mean(xper,2);
xper = xper./(2*pi);
xper = xper(1:length(xper)/2+1);
xper(2:end-1) = 2*xper(2:end-1);
max(abs(Pxx-xper))

You see that the maximum difference in absolute value is only on the order of 10^(-16).

Hope this helps

[Neraj]

Hi Wayne,

I copied and pasted your code (except for the first line, which gives an error). When I plotted the log of Pxx and xper on the same graph, I see the same offset issue. For me,

max(abs(Pxx-xper)) = 7.1751

Also, 0Hz and 1Hz do not match.

Am I missing something?

Thanks,

Neraj

[Peter]

hi Wayne, I tried to generalize the above code for cases where segment length and nfft are different successfully for the case where nfft is larger than segment length. See below, case 2. The case where nfft is smaller than the segment length (case 1) puzzles me though. Matlab 2017b help says: "If nfft is less than the segment length, the segment is wrapped using datawrap to make the length equal to nfft." but I can't understand what this means. Could you shine some light? (perhaps by means of a piece of code where the question marks are)

clear all
casenr=2
nt=1024;
switch casenr
    case 1 %segments longer than nfft, 2 segments                
        segment_length=nt/2;
        nfft=400;       
    case 2 %segments shorter than nfft, 2 segments                
        segment_length=400;
        nfft=512;        
end

nzeropad=nfft-segment_length;

x = randn(nt,1); 
[Pxx,w] = pwelch(x,segment_length,0,nfft);
xper = zeros(nfft,2);

for k = 0:1
    %take segment
    i1=1+k*segment_length;
    i2=(k+1)*segment_length;
    x1 = x(i1:i2);    
    %windowing + zero padding:
    if nfft>=segment_length
        windowlength=segment_length;
        xw = [x1.*hamming(windowlength);zeros(nzeropad,1)];
    elseif nfft<segment_length  
        %?????????    
    end
xper(:,k+1) = abs(fft(xw)).^2./norm(hamming(windowlength),2)^2;
end
xper = mean(xper,2);
xper = xper./(2*pi);
xper = xper(1:length(xper)/2+1);
xper(2:end-1) = 2*xper(2:end-1);
max(abs(Pxx-xper))

Answer 4

[Wayne King]

You copied and pasted this code?

x = randn(1024,1); 
Pxx = pwelch(x,512,0,512);
xper = zeros(512,2);
for k = 0:1
x1 = x(1+k*512:(k+1)*512);
xw = x1.*hamming(512);
xper(:,k+1) = abs(fft(xw)).^2./norm(hamming(512),2)^2;
end
xper = mean(xper,2);
xper = xper./(2*pi);
xper = xper(1:length(xper)/2+1);
xper(2:end-1) = 2*xper(2:end-1);
max(abs(Pxx-xper))

First clear everything in your workspace, then try it.

I get perfect agreement. What version of MATLAB are you using?

[Abusom]

What is the purpose of this line:

xper(2:end-1) = 2*xper(2:end-1);

I understand the code but want to know why you did so?

[Wayne King]

I am doubling everything in the one-sided spectral estimate except DC and the Nyquist.

[Abusom]

Yes, I understand that. But, why did you do so? I don't see this in the algorithm. Thanks

[Celso Coslop Barbante]

Also got agreement. This code was very nice to understand how to mimic pwelch PSD estimator with FFT. I am trying to implement a pwelch-like function in hardware using FPGA and FFT cores. Thanks.