137 views (last 30 days)

Show older comments

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

Moses
on 20 Sep 2016

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

Walter Roberson
on 30 Apr 2018

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

Peter
on 7 Nov 2018

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

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?

Celso Coslop Barbante
on 26 Jul 2012

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.

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!