You need to use Welch’s method for power spectral density estimation. Read the details here, and see the demo below.
clc, clear
x = randn(1024*16, 1);
% length of each segment
seg = 1024;
% FFT length
nfft = nextpow2(seg);
% Overlap 25% of window
ov = 0.25*1024;
% Analysis window
wind = hanning(seg);
% Welch’s power spectral density estimate
[pxx,f] = pwelch(x, wind, ov, [], 1, 'centered', 'psd');
% Plot
plot(f, 10*log10(pxx))
xlabel('Normalized frequency - Hz')
ylabel('PSD - dB')
title('Welch’s power spectral density')