calculate stop band energy from PSD graph

I want to find the value of the stop band energy i.e. the total leaked energy to the neighbouring sub-carriers. The value is calculated with the normalized total energy but I don't know how to find it. Below is the code of a simple rectangular pulse/filter.

%% Rectangular Filter

Ts = 0.01;                           % sampling time

time = -20:0.01:20;                  % time axis       

NFFT = 2048;                         % number of samples

freq = linspace(-0.5,0.5,NFFT);      % frequency axis

ht_Rectangular = zeros(size(time));

%% In Time Domain

for n = 1:length(time)
if time(n)>=-0.5 && time(n)<=0.5
ht_Rectangular(n) = 1;
end
end

figure;
plot(time,ht_Rectangular);
xlabel('t/T_s');
ylabel('h(t)');
title('Impulse Response of Rectangular Filter');
grid on,

%% In Frequency Domain

Hf_Rectangular = fftshift(fft(ht_Rectangular,NFFT));  % FFT with FFTshift for both negative & positive frequencies

figure;
plot(freq, (Hf_Rectangular(1:length(freq)))/max(Hf_Rectangular));
xlabel('f/F');
ylabel('H(F)');
title('Frequency Response of Rectangular Filter');
grid on,

%% Power Spectral Density

[PSD_Rectangular, F_Rectangular] = periodogram (ht_Rectangular/max(ht_Rectangular));

figure;
plot (F_Rectangular, 10 * log10 (PSD_Rectangular));
xlabel('f/F');
ylabel('PSD (dB)');
title('Power Spectral Density of Rectangular Filter');
grid on,