# 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

• Can you give an explanation for this code? It does not seem directly related to the main part of your question. – nispio Sep 8 '13 at 2:40
• Thank you @nispio for your reply, I appreciate your time :) The codes above is just what I have done so far. I know it is very fundamental, but I need to get started somehow. Here, from this code, I did not manage to get the energy spectrum as desired, like the reference given. But, I got the amplitude and power spectrum (a few versions of methodology, which gives slightly different results). For eg, amplitudes 1. 'mY = 2*abs(Ydft(1:nfft/2+1))' 2. 'mY_HamWnd = 2*abs(Ydft_HamWnd(1:nfft/2+1))' 3. 'mY2 = abs(Ydft); mY2= mY2(1:nfft/2+1); mY2 (2:end-1) = 2* mY2 (2:end-1);' – kitkit Sep 8 '13 at 3:23
• Can you provide any more details on what you are trying to do? Your question is about computing the energy spectrum, yet the passage you cited is about computing the energy in the harmonics of a given frequency. Since these are not the same thing, maybe a little extra context would be helpful. What goal are you trying to achieve, and what are the inputs to your system? – nispio Sep 8 '13 at 4:45
• Dear nispio, I am actually working on a operating condition (OC) categorization of a given current signal from a machine experiment. The basic idea is that given an operating condition, the current signal will show a unique time domain signal.Hence, this will translate into a unique and more informative frequency domain signal (in this case, I choose Fourier Transform as my tool). I think choosing an Energy spectrum will provide a more robust perception as compared to amplitude spectrum. That is the reason why I choose energy spectrum to look for signatures for the "OC" categorization. – kitkit Sep 8 '13 at 9:23
• I am unclear the difference on the difference between "energy spectrum" and "energy in harmonics of a given frequency". As said in the previous comment, I would like to obtain signatures that are able to differentiate the operating conditions of a given current signal. For example, to extract a contributing 6th harmonics signature, I would take its normalized ratio such as "energy ratio of the 2nd harmonic to the 6th harmonic" or something like "Normalized energy within the first four harmonics" etc. – kitkit Sep 8 '13 at 9:39

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>0)-(x<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)')

• The total (normalized) energy of the segment can be computed as simply Eyy = y(1:N)*y(1:N)'/N. Then the ratios can be computed as 2*Pyk/Eyy where the factor of 2 accounts for the fact that you are only looking at the positive frequency support. – nispio Sep 11 '13 at 20:32