Scaling of the PWELCH function in MATLAB

Question

I am trying to compute the power spectral density of a random signal using the PWELCH function in MATLAB.

Since I think have not understood properly how Pwelch scales the PSD, I wrote a sample program, in which I generate a sum of two sinusoids with given amplitudes - A1 and A2 - and given frequencies - f1 and f2 - superimposed to random noise.

I computed the PSD using the FFT and then using Pwelch. It seems to me that the PSD computed via FFT is correct, because it returns the power corresponding to the sinusoids (peaks at $A1^2$ and $A2^2$).

However Pwelch seems to have a different scaling because the peak values depend on the number of samples I put in the PWELCH window. This should make sense up to a certain extent, but the variation seems too large to me. I think I haven't properly understoon either how to use PWELCH or how to scale the PSD returned by PWLECH.

Can you help me understand this? Thanks. The sample code is reported below.

Thank you E.

---

clear all
close all
clc

rng default


Fs = 10000;            % Sampling frequency                    
T = 1/Fs;             % Sampling period       
L = 50000;             % Length of signal

t = (0:L-1)*T; 

A1=2;
A2=3;

f1=5;
f2=20;

X=A1*sin(2*pi*f1*t) + A2*sin(2*pi*f2*t)+ 2*randn(size(t));

figure(1)
hold on
grid on
plot(t,X)
xlabel('time (s)')
ylabel('signal (Unit)')

%% Power Spectrum via FFT
% Compute the two-sided spectrum P2. Then compute the single-sided spectrum P1 based on P2 and the even-valued signal length L

Y = fft(X);

P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);

f = Fs*(0:(L/2))/L;

%% via FFT

figure(2)
hold on
grid on
plot(f,P1.^2) 
title('Single-Sided Power Spectrum of X(t)')
xlim([0 50])
xlabel('f (Hz)')
ylabel('PSD (Unit^2/Hz)')
set(gca,'xscale','log')
% set(gca,'yscale','log')

%% via PWELCH

nfft= L;
window= rectwin(nfft);

[Pxx_f,f]=pwelch(X,window,[],nfft,Fs);

figure(2)
plot(f,Pxx_f,'r','linewidth',1.5);

legend('via FFT','via Pwelch','location','best')

Answer 1

I believe, there are two aspects of your question.

For the FFT calculation, the power density spectrum is the square of the FFT modulus multiplied by the sampling period (in your notation: T) and divided by the number of samples (in your notation:L): T/L * abs(FFT(X))^2 or T/L *abs(Y(f))^2

Here is now where I believe you have an error, you have to add the power of Y(f) and (Y(-f)) and not add first Y(f) and (Y(-f) and then perform the power calculation. By adding the amplitude first, in your notation, P1(2:end-1) = 2*P1(2:end-1) then after calculating the power, you have the double of the correct value.

A far as the Welch method for the calculation of the Power spectrum, it takes into account all the necessary normalization. You can chose between the power spectrum (option "power") or the power density spectrum (default option "psd")

Trying to use the option "power" for the welch calculation, it will gives as answer the amplitude of the sin squared divided by 2 (a sin as a power of 1/2): (A1^2)/2.

To calculate the power density spectrum, you have to divide the previous results by the applied frequency resolution, that is the sampling frequency by the number of samples, in your notation Fs/L. With the welch function, you have to use the option 'psd'. Sorry for the long answer and hoping, this helps.