MATLAB Answers

Neraj
0

quetsions about matlab pWelch implementation

Asked by Neraj
on 20 Feb 2012
Latest activity Edited by Moses
on 8 May 2019
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 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

  0 Comments

Sign in to comment.

Tags

4 Answers

Answer by Moses
on 20 Sep 2016
Edited by Moses
on 8 May 2019

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

  2 Comments

Hi! Is 6 in "for i = 0:6" the number of fft's?
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.

Sign in to comment.


Answer by Wayne King
on 20 Feb 2012

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

  2 Comments

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

Sign in to comment.


Answer by Wayne King
on 22 Feb 2012

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?

  4 Comments

Show 1 older comment
I am doubling everything in the one-sided spectral estimate except DC and the Nyquist.
Yes, I understand that. But, why did you do so? I don't see this in the algorithm. Thanks
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.

Sign in to comment.


Answer by Invizible Soul on 12 Nov 2012

Hello,
Why in the following line of the code
xper = xper./(2*pi);
there is a division by 2*pi ?
cordially.

  0 Comments

Sign in to comment.