[Matlab FFT]: FFT and ESD from Literature Definition

Hi all 

I am complete beginner in DSP and FFT in matlab. Thanks for helping guys !

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) in the following definition:
_____________________________________________________________________________

I stumbled on this explanation to obtain the energy spectrum, which spells out my intention, from an IEEE paper(Open Circuit Fault Diagnosis in 3 phase uncontrolled rectifiers, 2012, 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+1will 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 could result in the kth harmonic peak shifting away from bin 4k+1."
________________________________________________________________________________

I do not fully understand the passage above. In my limited understanding, do they refer  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 ?

>> I am complete beginner in DSP and FFT in matlab

stop reading IEEE papers. 
They give you loads of equations about denaturation properties of proteines
on lipophile surfaces in high-temperature environments. 
What you need is a recipe for scrambled eggs.

There are many articles, tutorials, free DSP textbooks on the web.	 

_____________________________		
Posted through www.DSPRelated.com

On Mon, 16 Sep 2013 01:29:19 -0700 (PDT), kit <cktang90@gmail.com>
wrote:

>Hi all 
>
>I am complete beginner in DSP and FFT in matlab. Thanks for helping guys !
>
>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) in the following definition:
>_____________________________________________________________________________
>
>I stumbled on this explanation to obtain the energy spectrum, which spells out my intention, from an IEEE paper(Open Circuit Fault Diagnosis in 3 phase uncontrolled rectifiers, 2012, 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+1will 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 could result in the kth harmonic peak shifting away from bin 4k+1."
>________________________________________________________________________________
>
>I do not fully understand the passage above. In my limited understanding, do they refer  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 ?

Hello cktang,

Let's review the quote one paragraph at a time.

   "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)."

That paragraph, by itself is not clear at all.  After 
reading the other pragraphs I think I understand what
they're saying.

They capture samples over a time duration of 4T seconds.
So if there's a periodic wave whose period is T in
the signal, then there will be four repetitions of 
the T-period waveform in the 4T-second duration 
sampled signal.

   "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+1will coincide with the 
   kth harmonic frequency."

If the 4T-second duration digital signal contained 
exactly four cycles of a single sine wave (whose 
period is T seconds), then the spectral energy would 
show up at the 4+1 = 5th bin.  (They are considering 
the very first FFT bin to be the zeroth, 0th, bin!)

If the 4T-second duration digital signal contained 
exactly four cycles of a single sine wave (whose 
period is T seconds) and exactly eight cycles of another 
sine wave (whose period is T/2 seconds), then the spectral 
energy would show up at the 4+1 = 5th bin *AND* 
the 8+1 = 9th bin.


If the 4T-second duration digital signal contained 
exactly four cycles of a single sine wave (whose 
period is T seconds) and exactly eight cycles of another 
sine wave (whose period is T/2 seconds) and exactly 
twelve cycles of a single sine wave (whose period is 
T/3 seconds), then the spectral energy would show up at 
the 4+1 = 5th bin *AND* the 8+1 = 9th bin*AND* the 
12+1 = 13th bin.

That last scenerio is when we have one cycle of one sine 
wave, two cycles of another sine wave, and three cycles of 
a 3rd sine wave in atime duration of T seconds (one fourth 
the duration of the total digital signal).

So in that last scenario there's spectral energy at 
the 5th bin (k=1, so 4k+1=5), the 9th bin (k=2, so 4k+1=9), 
and the 13th bin (k=3, so 4k+1=13).  So by their 
definition:

   5th bin is the 1st harmonic (k=1)
   9th bin is the 2nd harmonic (k=2)
   13th bin is the 3rd harmonic (k=3)

The quote's 3rd & 4th paragraphs (which should 
have been combined into a single paragraph) are:

   "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 could result in 
   the kth harmonic peak shifting away from bin 4k+1."

What they're saying is that if the three sine waves in 
our 3rd scenario had slightly more, or slightly less, 
than four, eight, and twelve cycles in the T-second 
duration didigtal signal, then  spectral leakage will 
occur and the spectral energy from a single sine wave 
will "smear" above and below the 5th, 9th, and 13th 
bins. 

To "collect" most of the spectral energy of the 
lowest-freq sine wave, they add the bin values of 
five bins centered at the 5th bin.  To "collect" most 
of the spectral energy of the highest-freq sine wave, 
they add the bin values of five bins centered at 
the 13th bin.

I'm too lazy to try to descibe spectral leakage in 
my reply here.  It's covered in every DSP textbook.

Good Luck,
[-Rick-]

Dear Rick 

 Thank you very much for your *detailed* explanation ! It clears my doubt and affirms my understanding ! Feels so grateful to have a textbook author replying me ! I think this is a very good supplementary platform for learning, besides reading !

 Have been staring at your explanation for a few days and got some help for the following matlab code (thx nispio from dspstackexchange). However, I do understand the idea of applying a window function at time domain before fft to reduce spectral leakage !
----------------------
%% 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 |x|^2
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
----------------------

My objective is that I am working on categorizing operating conditions of a given current signal,x(t), that the basic idea is: 
" Given an operating condition, the current signal will show a unique periodic time domain signal.Hence, this will translate into a unique and more informative frequency domain signal."

Using its resulting spectral energy, and the features described above, the operating condition of the machine can be determined. 

My plan is:
1. I have a periodic time domain signal, f0 = 50Hz sampled (Fs = 20000Hz) from a machine experiment, x(t).
2. I am using a temporal sliding hamming window to split x(t) into frames, calculate each frame's spectral energy, continue for the rest of the signal, then average all frame's energy by the total sum of window.
3.  Then, I wish to calculate the spectral's energy ratio, E1(1st harmonic to total energy of the frame) and E123(first 3 harmonics to total energy of the frame). These are the features that are used to categorized the operating condition of a signal.

----------------------------
%continue from the snippet above
Pyk_tot = sum(Pyk(1:end)); % Sum of energies in signal segment 

E1 = Pyk(1)/Pyk_tot  % normalised energy ratio 1st harmonic to total energy in segment
E3 = (Pyk(1)+Pyk(2)+Pyk(3))/Pyk_tot % normalised energy ratio first 3 harmonics to total energy in segment
E1st = Pyk(1) % Energy of the 1st harmonic 
E2nd = Pyk(2) % Energy of the 2nd harmonic
E3rd = Pyk(3) %  Energy of the 3rd harmonic 
E1_2= E1st/E2nd  %  energy ratio 1st harmonic to 2nd harmonic 
E1_3= E1st/E3rd  %  energy ratio 1st harmonic to 3rd harmonic 
---------------

Another part is 
--------------
%% Buffer it and window
win = hamming(nfft);%//chosen window type based on application
y = buffer(y, nfft, nfft/2); %// 50% overlap between frames in this instance
y = y(:, 2:end-1); %//optional step to remove zero padded frames
y = (  y' * diag(win)  )'; %//efficiently window each frame using matrix algebra

% Calculate mean FFT
F = abs(  fft(y, nfft)  /  sum(win)  );
F = mean(F,2);
--------------

 The codes are not working perfectly as described in the plan above (pt 1,2,3).
 
 Or, better, is there any Matlab library that can be used to work on this   categorizing operating conditions as stated in objectives ?
 Things like: 
 - FFT(Autocorrelation of x(t)) = PSD, Is x(t) from DAQ experimentally is WSS random process ?
 - Hpsd=dspdata.psd(Data)
   Hs=spectrum.estmethod(input1,...)
   periodogram
   pwelch etc etc ?


 Sorry for the long reply with lots of question !
 Thank you, may you be well and happy !

~ck