Energy Spectrum of a signal after FFT in Matlab I have a vector of signal x(t) with its time vector. I want to obtain a frequency representation of the signal, in particular the energy spectrum of x(t).



Can someone please show some light into obtaining the energy spectrum ?
Is energy spectrum and energy spectral density ,ESD produce the same result ?



Any help is much appreciated !

 t = 0:.001:.25;
x = 5*sin(2*pi*50*t) + 1.2*sin(2*pi*120*t);
y = x + 2*randn(size(t));
y=y';
Fs = 1000; % Sampling frequency

% Fast Fourier Transform 
    L = length (y); % Window Length of FFT
    nfft = 2^nextpow2(L); % Transform length , Next power of 2 from length of signal
                          %  Exactly same as nfft = pow2(nextpow2(L)) 

% Take fft, padding with zeros so that length(input) is equal to nfft     
    Ydft = fft(y,nfft)/L;

    y_HamWnd = y.*hamming(L);           %hamming window fft
    Ydft_HamWnd = fft(y_HamWnd,nfft)/L;

%FFT is symmetric, throw away second half

    % Take the magnitude of fft of x
%Amplitude Spectrum 01 (Single Sided)
mY = 2*abs(Ydft(1:nfft/2+1)) ;

%Amplitude Spectrum 02 (Single Sided HamWmd)
mY_HamWnd = 2*abs(Ydft_HamWnd(1:nfft/2+1));

%Amplitude Spectrum 03 (Single Sided; Exclude DC-- FoldingFreq)
mY2 = abs(Ydft);
mY2= mY2(1:nfft/2+1);
mY2 (2:end-1) = 2* mY2 (2:end-1); 


    %To get the Power
% P ; single sided; then coefficient square**********Multiply by 2 ?
YdftCoeffi_Sqr = Ydft(1:nfft/2+1).^2;   

% P = y*(conj y) == y^2  ;     %********DO WE NEED TO MUL 2 ?
Pyy_NormNfft = Ydft.*conj(Ydft);   

% P ; single sided; First square, then Multiply by 2
m1YL = abs(Ydft(1:nfft/2+1));  
mSqrYL = m1YL.^2;  
m2YL = mSqrYL*2; 

% Frequency vector
f1 = Fs /2 * linspace(0,1, nfft/2 + 1);
f2 = (0:nfft-1)*(Fs/nfft);    

% Generate the plot, title and labels. 

figure;
plot(t,y); title('Time Domain'); xlabel('Time(s)'); 

figure;
subplot(3,1,1)
plot(f1,mY)
title('Amp Spec 01 --mY = 2*abs(Ydft(1:nfft/2+1)) /L; '); xlabel('Frequency (Hz)'); ylabel('|mY|')
        %axis ([0 500 0 5]); %Zoom in 
subplot(3,1,2)
plot(f1,mY_HamWnd)
title('Amp Spec 02 -- mYHamWnd =2*abs(YdftHamWnd(1:nfft/2+1))/L;'); xlabel('Frequency (Hz)'); ylabel('|mY_HamWnd|')
        %axis ([0 500 0 5]); %Zoom in 
subplot(3,1,3)
plot(f1,mY2)
title('Amp Spec 03 DC-- FoldingFreq -- mY2 ;First, single sided; then square, lastly Multiply by 2'); xlabel('Frequency (Hz)'); ylabel('|mY2|')
        %axis ([0 500 0 5]); %Zoom in 

figure;
subplot(3,1,1)
plot(f1,YdftCoeffi_Sqr); 
title('Ydft Complex Coefficient Square; Single Sided-- YdftCoeffi_Sqr = Ydft(1:nfft/2+1).^2');  xlabel('Frequency (Hz)'); ylabel('YdftCoeffi_Sqr');
subplot(3,1,2)
plot(f1,m2YL); 
title('Ydft Complex Coefficient Magnitude Square; Single sided; then Multiply by 2'); 
xlabel('Frequency (Hz)'); 
ylabel('Power m2YL');
        %axis ([0 500 0 5]); %Zoom in
The cited passage gives a method for obtaining the energy in the harmonics of a specific frequency.  If your frequency of interest is $f_0$, then take a section of samples of length $N$ where $N = \frac{4 F_s}{f_0}$ and compute the $N$-point FFT.  Now you can find the energy of the second harmonic by squaring the magnitude of the $4k+1 = 9$th FFT bin. To account for errors in estimating $f0$ you can sum the squared magnitudes of bins 7 through 11.

Below is some example code that might give an idea of what the algorithm in question does:

%% Generate an example signal
Fs = 5e3;                     % Sampling rate
f0 = 151;                     % Frequency of interest
t = (0:Fs-1)/Fs;              % Time samples
x = cos(2*pi*f0*t);           % Pure sinusoid @ f0
noise = 0.1*randn(size(x));   % Gaussian noise
y = (x&gt;0)-(x&lt;0) + noise;      % Introduce harmonics and noise

%% Compute energy in the harmonics of f0
N = round(4*Fs/f0);           % Determine record length
Y = fft(y(1:N))/N;            % Compute N-FFT (normalized)
Y2 = Y.*conj(Y);              % Magnitude squared
K = floor(N/2/4);             % Number of harmonic "bins"
for k = 1:K
    idx =  4*k+1+(-2:2);      % Five indices centered at 4k+1
    Pyk(k) = sum(Y2(idx));    % Sum of energies around 4k+1
end

%% Plot the results
[Pyy,Fyy] = pwelch(y,[],[],[],Fs);
plot(Fyy,10*log10(Pyy));
hold on;
F = linspace(0,Fs-Fs/N,N);
Fk = F(1+4*(1:K));
plot(Fk,10*log10(Pyk),'ro');
title('Power Spectral Density');
xlabel('Frequency (Hz)')
ylabel('Power/Frequency (dB/rad/sample)')
data2 = sprintf('Energy in Harmonics of %d Hz',f0);
legend('Power Spectral Density Estimate', data2);
hold off;




As far as finding a "signature" in the data, you can see that the algorithm could be very useful for that.  However, it is very important to note that the algorithm presumes some fundamental frequency of the signal in question.  If the signal in question has a fixed period, but the shape of the signal over that period changes with time, then this algorithm would be useful.  However, if the frequency varies with time, then this algorithm is not a good fit for your application.

If you are simply looking to estimate the spectrum of a signal, then maybe this page will be of help: FFT for Spectral Analysis (mathworks.com)