Matlab – How to obtain Energy spectrum of a signal after FFT in Matlab

EDIT:

I stumbled on this explanation to obtain the energy spectrum from an IEEE paper(Open Circuit Fault Diagnosis in 3 phase uncontrolled rectifiers, Rahiminejad, Diduch, Stevenson, Chang).
"A recorded sample of the signal containing a number of samples equivalent to 4T is captured and its FFT is determined using an FFT size equal to the record length (where T is the fundamental period).
Assuming the FFT size is matched to 4 periods of a periodic waveform, every 4th FFT bin will coincide with a harmonic frequency, in particular the center of FFT bin 4k+1 will coincide with the kth harmonic frequency.

The energy of the kth harmonic is calculated as the sum of the squared magnitudes of the 5 consecutive FFT values centred at bin 4k+1. The additional FFT values are included in the harmonic energy calculation so as to reduce the sensitivity of the calculated energy to an error in the frequency estimate which oculd result in the kth harmonic peak shifting away from bin 4k+1."

I do not fully understand the passage above. In my limited understanding, the bold line
refers to the sum of the squared magnitudes of the output of function fft(), i.e. the complex fourier series coefficients.

Can someone please show some light into obtaining the energy spectrum ?
@fpe : I am not sure if ESD performs the same as energy spectrum. BTW, thanks alot for your answer:)

I am trying to plot the energy spectrum of a signal to look at the for example Normalised energy contained first three harmonics, energy ratio of fundamental to 2nd harmonics etc….

Here I have managed to get the Hanning window FFT amplitude-Hz and power-Hz. But, I have no idea how to get energy-Hz for each frequency components.

Any help is much appreciated !

function [f,Xall_Wnd]=fftplotExxx(time,X_input)
Fs = 20000; % Sampling frequency
x = X_input; 

% Fast Fourier Transform 
    L = length (X_input); % Length of FFT
    nfft = 2^nextpow2(L); % Next power of 2 from length of signal

%wave = wave.*hamming(length(wave));
x_HammingWnd = x.*hamming(L);

% Take fft, padding with zeros so that length(X) 
%is equal to nfft 

Xall_Wnd = fft(x_HammingWnd, nfft)/L;  %hamming window fft

% FFT is symmetric, throw away second half

% Take the magnitude of fft of x 
mx_Wnd = 2*abs(Xall_Wnd(1:nfft/2+1)); 

% To get Power of x(t) by sqr of magnitude
m2x_Wnd = mx_Wnd.^2; 

% I am Not sure how to get energy spectrum
for i=1:L:nfft-L
E(i) = sum(Xall_Wnd(1:nfft/2+1).^2);
end

% Frequency vector
    f = Fs/2*linspace(0,1,nfft/2+1);

% Generate the plot, title and labels. 
subplot(2,2,1)
plot(time,X_input);
title('Time Domain')
xlabel('Time(s)')

subplot(2,2,2)
plot(f,m2x_Wnd); 
title('Power Spectrum of x(t)'); 
xlabel('Frequency (Hz)'); 
ylabel('Normalised Power of fft');

subplot(2,2,3)
plot(f,mx_Wnd); 
title('Hamming Window_ Spectrum of x(t)'); 
xlabel('Frequency (Hz)'); 
ylabel('Normalised Magunitude of fft');

subplot(2,2,4)
plot(f,E); 
title('Energy Spectrum of x(t)'); 
xlabel('Frequency (Hz)'); 
ylabel('Energy');
end