Why do results differ between pwelch and dsp.SpectrumEstimator (random noise example) ?

I am refering to an already existing question : https://fr.mathworks.com/matlabcentral/answers/338948-why-do-results-differ-between-pwelch-and-dsp-spectrumestimator

for which an interesting example was provided by Mathworks Support Team using a sine wave as the input test signal. This example clearly shows that there is a right way to implement the dsp.SpectrumEstimator to match the pwelch function.

However, in this example I tried to change the test signal to a random noise using :

x = randn(nPoints,1);

I was expecting a very good match between the two methods (like a perfect match), but there are some differences, and I would like to understand where they are coming from.

Any ideas ?

Thanks !

Answer 1

[Shivam Gothi]

Hello @Olivier Schevin,

As per my understanding, you're attempting to calculate the Power Spectral Density (PSD) of a signal using both the “pwelch” function and the “dsp.SpectrumEstimator” object method. However, it seems that the output signals from these two approaches are not identical, even though both utilize the same method to determine the PSD.

The explanation for this is as follows:

• In “dsp.SpectrumEstimator” object, we need to manually segment the signal and step this segment through the object.
• “pwelch” method automatically divides the signal into segments, equal in length. The length of each segment is specified using the “window” property of “pwelch” function.

But, there is an additional input argument for “pwelch” function called “noverlap”. Refer to the below documentation:

https://www.mathworks.com/help/releases/R2024a/signal/ref/pwelch.html?s_tid=doc_ta

The “pwelch” function uses the value of “noverlap” and decides the number of samples to be overlapped between the segments. The script provided in the MATLAB answer link attached by you, specifies “noverlap” as “empty array”. If you do not specify “noverlap”, or specify “noverlap” as empty, the default number of overlapped samples is 50% of the window length. But in the case of “dsp.SpectrumEstimator” object, no samples are overlapped.

Therefore, in the present script, 50% of the samples from segments are getting overlapped in case of “pwelch” function, but no samples are overlapped in case of “dsp.SpectrumEstimator” object.

Due to above stated reason, you are not getting the same output for both the methods.

To resolve the issue, I specified “noverlap” input argument for “pwelch” function as “0”. I am attaching the code below to demonstrate this. Here, I am taking input signal “X” as suggested by you in the question.

fs = 8192; % sampling rate (Hz)
dt = 1/fs;  % sample time (s)

nfft = 8192;
fWindow = hann(nfft, 'periodic');
freqRange = 'onesided';
spectrumType = 'power'; % 
spectrumTypeDsp = 'Power'; %'Power density';

nrec = 100;  %  Obtain the spectral average of nrec segments 

% total number of samples
nPoints = nrec*nfft;

%% Create test signal 
% sine wave at 1000 Hz with amplitude 5.3 with white noise
t = dt*(0:nPoints-1)';
A = 5.3;
fSine = 1000;

rng(42)
x = randn(nPoints,1);

%%
% Calculate pwelch estimate. set "noverlap" to 0
[PxxTrue,fTrue] = pwelch(x,fWindow, 0,nfft,fs,freqRange,spectrumType); 


% demonstration of correct implementation 


seObject2 = dsp.SpectrumEstimator('SampleRate',fs,...
    'SpectrumType',spectrumTypeDsp,'Window', 'hann',...
    'FFTLengthSource','Property','FFTLength',nfft,...
    'SpectralAverages',nrec,'FrequencyRange',freqRange,Method='Welch');

%%

% Correct method to obtain the sepctral estimate 
for idx = 1:nrec
    PxxEst = seObject2( x((idx-1)*nfft+1 : idx*nfft) ) ; 
end

fEst = getFrequencyVector(seObject2);

%%
% Plot the results
figure
plot(fTrue,db(PxxTrue),fEst,db(PxxEst))
legend('pwelch','dsp.SpectrumEstimator')
xlabel('Frequency (Hz)')
ylabel('Power (dB)')
title('Correct Welch Power Spectrums ')

Conclusion:

You can see that both the outputs “match perfectly” as desired by you.

I hope this answers your question!

[Olivier]

Many thanks @Shivam Gothi, it matches perfectly indeed !