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