Spectral Analysis of Electricity Demand Using Hilbert–Huang Transform

The large amount of sensors in modern electrical networks poses a serious challenge in the data processing side. For many years, spectral analysis has been one of the most used approaches to extract physically meaningful information from a sea of data. Fourier Transform (FT) and Wavelet Transform (WT) are by far the most employed tools in this analysis. In this paper we explore the alternative use of Hilbert–Huang Transform (HHT) for electricity demand spectral representation. A sequence of hourly consumptions, spanning 40 months of electrical demand in Spain, has been used as dataset. First, by Empirical Mode Decomposition (EMD), the sequence has been time-represented as an ensemble of 13 Intrinsic Mode Functions (IMFs). Later on, by applying Hilbert Transform (HT) to every IMF, an HHT spectrum has been obtained. Results show smoother spectra with more defined shapes and an excellent frequency resolution. EMD also fosters a deeper analysis of abnormal electricity demand at different timescales. Additionally, EMD permits information compression, which becomes very significant for lossless sequence representation. A 35% reduction has been obtained for the electricity demand sequence. On the negative side, HHT demands more computer resources than conventional spectral analysis techniques.


Introduction
Modern electrical networks are one of the main scenarios in which a widespread deployment of sensors is being achieved [1], thus acquiring a large amount of information of different types. Processing that information fosters making more efficient and intelligent decisions, in such a way that resulting power networks are usually known as smart grids [2].
The large number of sensors used in today s smart grids has become economically and technically feasible due to a significant reduction in their cost, the availability of higher capacity communication systems, and the greater accessibility to storage and processing equipment [3].
The straightforward availability of greater datasets poses a challenge on the data analytic side [4] making possible new and more efficient approaches to several applications such as fault detection [5], predictive maintenance [6], transient stability analysis [7], electric device state estimation [8], power quality monitoring [9], topology identification [10], renewable energy forecasting [11], and non-technical loss detection [12].
Especially relevant are the smart grid applications in which a high ratio of renewable and distributed energies has to be considered. Autoregressive integrated moving average (ARIMA) has In Figure 1b, a ten week evolution (March and April 2016) is depicted, showing a typical weekly behavior with high electricity consumption on labor days and lower values on weekends. Finally, Figure 1c shows an example for a 15-day period, where the two daily peaks (typically at noon and mid-afternoon) are shown.

Conventional Techniques for Spectral Analysis
Frequency-domain analysis of the dataset is primarily tackled by employing Fourier Transform (FT). Calling ( ) the value of the electricity demand at time , its frequency representation ( ) is defined as: When ( ) is digitized at intervals, that is, at a sampling rate = 1/ , a sequence of values is obtained, where ( ) denotes the electricity demand at the -th sampling time. For discrete sequences, frequency-domain representation is obtained using the Discrete Fourier Transform (DFT): where ( ) is a complex value representing the amplitude and phase of the harmonic component at frequency / . The graph representing the absolute value | ( )| is usually known as the spectrum of . Instead of DFT, a faster version of the algorithm, called Fast Fourier Transform (FFT), is commonly used. As the dataset described in the previous subsection contains hourly load values, sampling period is = 1 h and, then, sampling frequency is = 1 hour = 8760 year . In Figure 1b, a ten week evolution (March and April 2016) is depicted, showing a typical weekly behavior with high electricity consumption on labor days and lower values on weekends. Finally, Figure 1c shows an example for a 15-day period, where the two daily peaks (typically at noon and mid-afternoon) are shown.

Conventional Techniques for Spectral Analysis
Frequency-domain analysis of the dataset is primarily tackled by employing Fourier Transform (FT). Calling x(t) the value of the electricity demand at time t, its frequency representation X( f ) is defined as: ( When x(t) is digitized at T s intervals, that is, at a sampling rate f s = 1/T s , a sequence of N values is obtained, where x(n) denotes the electricity demand at the n-th sampling time. For discrete sequences, frequency-domain representation is obtained using the Discrete Fourier Transform (DFT): where X(k) is a complex value representing the amplitude and phase of the harmonic component at frequency k f s /N. The graph representing the absolute value X(k) is usually known as the spectrum of x. Instead of DFT, a faster version of the algorithm, called Fast Fourier Transform (FFT), is commonly used. As the dataset described in the previous subsection contains hourly load values, sampling period is T s = 1 h and, then, sampling frequency is f s = 1 hour −1 = 8760 year −1 . If x(n) cannot be considered a stationary sequence, it is more convenient to consider the DFT of a full sequence's short slice: where w(n − m) is a certain window function centered at the m-th sampling time. In this research, the Hamming window function has been used [39]. X(m, k) is a discrete complex sequence that it is commonly represented by the square of its absolute value. The resulting graph is called the sequence's spectrogram. Although Short-Time Fourier Transform (STFT) is widely used in many technical fields, it suffers a serious limitation on time-frequency resolution. Calling n w the number of values (width) of the window function, then time interval ∆t (time resolution) and frequency interval ∆ f (frequency resolution), are: Joint time-frequency resolution is limited because: that is, the greater the resolution in time, the smaller in frequency. Making window function wider (greater value of n w ), time interval ∆t is increased (time resolution decreased) and frequency interval ∆ f is decreased (frequency resolution increased).
To overcome this problem many researchers use the Wavelet Transform (WT), defined as: where ψ[·] is a time-limited function called a wavelet, which is shifted at time mT s (m is an integer) and time-scaled by a factor s (a real positive value). X(m, s) is a function, discrete in m and continuous in s, that it is commonly represented by the square of its value. The resulting graph is called the sequence's scalogram where the m-axis represents the time and the s-axis is inversely proportional to the frequency. In this research, the Morlet wavelet has been used, which is defined as: ψ[ξ] ≡ e − ξ 2 2 cos(5ξ).
Wavelet transform is widely used for analyzing non-stationary sequences. In the case of stationary sequences, we are no longer interested in the evolution over time of its spectral behavior as they do not change. Thus, it is better to obtain a time-invariant spectrum using the Marginal Wavelet Transform (MWT) which is defined as:

The Hilbert Transform for Spectral Analysis
In some harmonic analysis applications, it is sometimes preferred to rely on an alternative approach called Hilbert Transform (HT), defined as: Sensors 2020, 20, 2912 In most applications, the resulting Hilbert Transformx(n) is combined with the original function, obtaining the so-called analytic function, defined as: which is a time-dependent complex sequence. Its module is called instantaneous amplitude, that is, Sensors 2020, 20, x FOR PEER REVIEW 5 of 25 In most applications, the resulting Hilbert Transform ̂( ) is combined with the original function, obtaining the so-called analytic function, defined as: which is a time-dependent complex sequence. Its module is called instantaneous amplitude, that is, ( ) ≡ |̃( )| . An instantaneous frequency ( ) can also be obtained using the phase of the analytic function: Combining the instantaneous amplitude and frequency in a single vector of two functions, the Hilbert Spectrum (HS) is obtained, that is: For the case of stationary sequences, it is better to obtain the Marginal Hilbert Spectrum (MHS), which is defined as: (14)

Empirical Decomposition and the Hilbert-Huang Transform for Spectral Analysis
In many practical sequences it is not uncommon that some Hilbert instantaneous frequencies (and their derived HS and MHS) have negative values, which are meaningless [40]. For the instantaneous frequency to always be positive, the sequence must satisfy two conditions [41]: 1. It is symmetric with respect to a null local average value; and 2. It has the same number of zero-crossing and extrema.
A sequence fulfilling these requirements is called an Intrinsic Mode Function (IMF). Therefore, to obtain meaningful Hilbert Transform, the original sequence ( ) should be first decomposed in several IMFs ( ). A method to achieve this goal is named Empirical Mode Decomposition (EMD) which is based on two nested loops: An outer iteration to obtain the -th IMF ( ); and an inner iteration to obtain ℎ , ( ), the -th trial for the -th IMF. The EMD algorithm (also called a sifting process) can be summarized in the following steps [42]: In most applications, the resulting Hilbert Transform ̂( ) is combined with the original function, obtaining the so-called analytic function, defined as: which is a time-dependent complex sequence. Its module is called instantaneous amplitude, that is, ( ) ≡ |̃( )| . An instantaneous frequency ( ) can also be obtained using the phase of the analytic function: Combining the instantaneous amplitude and frequency in a single vector of two functions, the Hilbert Spectrum (HS) is obtained, that is: For the case of stationary sequences, it is better to obtain the Marginal Hilbert Spectrum (MHS), which is defined as: (14)

Empirical Decomposition and the Hilbert-Huang Transform for Spectral Analysis
In many practical sequences it is not uncommon that some Hilbert instantaneous frequencies (and their derived HS and MHS) have negative values, which are meaningless [40]. For the instantaneous frequency to always be positive, the sequence must satisfy two conditions [41]: 1. It is symmetric with respect to a null local average value; and 2. It has the same number of zero-crossing and extrema.
A sequence fulfilling these requirements is called an Intrinsic Mode Function (IMF). Therefore, to obtain meaningful Hilbert Transform, the original sequence ( ) should be first decomposed in several IMFs ( ).
A method to achieve this goal is named Empirical Mode Decomposition (EMD) which is based on two nested loops: An outer iteration to obtain the -th IMF ( ); and an inner iteration to obtain ℎ , ( ), the -th trial for the -th IMF. The EMD algorithm (also called a sifting process) can be summarized in the following steps [42] 8. Obtain the -th trial for the -th proto-IMF subtracting the mean value from the original data: Combining the instantaneous amplitude and frequency in a single vector of two functions, the Hilbert Spectrum (HS) is obtained, that is: Sensors 2020, 20, x FOR PEER REVIEW In most applications, the resulting Hilbert Transform ̂( function, obtaining the so-called analytic function, defined as: which is a time-dependent complex sequence. Its module is calle ( ) ≡ |̃( )| . An instantaneous frequency ( ) can also be analytic function: Combining the instantaneous amplitude and frequency in a Hilbert Spectrum (HS) is obtained, that is: For the case of stationary sequences, it is better to obtain the M which is defined as:

Empirical Decomposition and the Hilbert-Huang Transform for Sp
In many practical sequences it is not uncommon that some (and their derived HS and MHS) have negative values, whic instantaneous frequency to always be positive, the sequence must 1. It is symmetric with respect to a null local average value; and 2. It has the same number of zero-crossing and extrema.
A sequence fulfilling these requirements is called an Intrinsi to obtain meaningful Hilbert Transform, the original sequence several IMFs ( ).
A method to achieve this goal is named Emp which is based on two nested loops: An outer iteration to obtai iteration to obtain ℎ , ( ), the -th trial for the -th IMF. The EM process) can be summarized in the following steps [42]  is combined with function, obtaining the so-called analytic function, defined as: which is a time-dependent complex sequence. Its module is called instantaneous amplit ( ) ≡ |̃( )| . An instantaneous frequency ( ) can also be obtained using the ph analytic function: Combining the instantaneous amplitude and frequency in a single vector of two fu Hilbert Spectrum (HS) is obtained, that is: For the case of stationary sequences, it is better to obtain the Marginal Hilbert Spect which is defined as:

Empirical Decomposition and the Hilbert-Huang Transform for Spectral Analysis
In many practical sequences it is not uncommon that some Hilbert instantaneous (and their derived HS and MHS) have negative values, which are meaningless [4 instantaneous frequency to always be positive, the sequence must satisfy two conditions 1. It is symmetric with respect to a null local average value; and 2. It has the same number of zero-crossing and extrema.
A sequence fulfilling these requirements is called an Intrinsic Mode Function (IMF) to obtain meaningful Hilbert Transform, the original sequence ( ) should be first dec several IMFs ( ).
A method to achieve this goal is named Empirical Mode Decompos which is based on two nested loops: An outer iteration to obtain the -th IMF ( ); a iteration to obtain ℎ , ( ), the -th trial for the -th IMF. The EMD algorithm (also call process) can be summarized in the following steps [42]: 1. Consider original data as the first residue: 1 ( ) = ( ) 2. Make = 1 3. Obtaining the -th IMF: ( ) 1. Consider the -th residue as the first trial for the -th proto-IMF: ℎ ,1 ( ) = ( ) 2. Make = 1 3. Consider the -th trial for the -th proto-IMF as the ongoing sequence: , ( ) = 4. Obtain upper local extrema of , ( ) For the case of stationary sequences, it is better to obtain the Marginal Hilbert Spectrum (MHS), which is defined as: In most applications, the resulting Hilbert Transform ̂( ) is combined with the origina function, obtaining the so-called analytic function, defined as: which is a time-dependent complex sequence. Its module is called instantaneous amplitude, that i ( ) ≡ |̃( )| . An instantaneous frequency ( ) can also be obtained using the phase of th analytic function: Combining the instantaneous amplitude and frequency in a single vector of two functions, th Hilbert Spectrum (HS) is obtained, that is: For the case of stationary sequences, it is better to obtain the Marginal Hilbert Spectrum (MHS which is defined as: (14

Empirical Decomposition and the Hilbert-Huang Transform for Spectral Analysis
In many practical sequences it is not uncommon that some Hilbert instantaneous frequencie (and their derived HS and MHS) have negative values, which are meaningless [40]. For th instantaneous frequency to always be positive, the sequence must satisfy two conditions [41]: 1. It is symmetric with respect to a null local average value; and 2. It has the same number of zero-crossing and extrema.
A sequence fulfilling these requirements is called an Intrinsic Mode Function (IMF). Therefor to obtain meaningful Hilbert Transform, the original sequence ( ) should be first decomposed i several IMFs ( ).
A method to achieve this goal is named Empirical Mode Decomposition (EMD which is based on two nested loops: An outer iteration to obtain the -th IMF ( ); and an inne iteration to obtain ℎ , ( ), the -th trial for the -th IMF. The EMD algorithm (also called a siftin process) can be summarized in the following steps [ Combining the instantaneous amplitude and frequency Hilbert Spectrum (HS) is obtained, that is: For the case of stationary sequences, it is better to obtain which is defined as:

Empirical Decomposition and the Hilbert-Huang Transform f
In many practical sequences it is not uncommon that s (and their derived HS and MHS) have negative values, instantaneous frequency to always be positive, the sequence 1. It is symmetric with respect to a null local average value 2. It has the same number of zero-crossing and extrema.
A sequence fulfilling these requirements is called an In to obtain meaningful Hilbert Transform, the original sequen several IMFs ( ).
A method to achieve this goal is named which is based on two nested loops: An outer iteration to iteration to obtain ℎ , ( ), the -th trial for the -th IMF. Th process) can be summarized in the following steps [42]: 1. Consider original data as the first residue:

Empirical Decomposition and the Hilbert-Huang Transform for Spectral Analysis
In many practical sequences it is not uncommon that some Hilbert instantaneous frequencies (and their derived HS and MHS) have negative values, which are meaningless [40]. For the instantaneous frequency to always be positive, the sequence must satisfy two conditions [41]:

1.
It is symmetric with respect to a null local average value; and 2.
It has the same number of zero-crossing and extrema.
A sequence fulfilling these requirements is called an Intrinsic Mode Function (IMF). Therefore, to obtain meaningful Hilbert Transform, the original sequence x(n) should be first decomposed in several IMFs c i (n). A method to achieve this goal is named Empirical Mode Decomposition (EMD) which is based on two nested loops: An outer iteration to obtain the i-th IMF c i (n); and an inner iteration to obtain h i,j (n), the j-th trial for the i-th IMF. The EMD algorithm (also called a sifting process) can be summarized in the following steps [42]:
Obtaining the i-th IMF: c i (n)

1.
Consider the i-th residue as the first trial for the i-th proto-IMF: h i,1 (n) = r i (n) 2.
Consider the j-th trial for the i-th proto-IMF as the ongoing sequence: Obtain upper local extrema of x i,j (n) Sensors 2020, 20, 2912 6 of 25
Obtain lower local extrema and lower envelope: l i,j (n) 7.
Obtain mean value of upper and lower envelopes: m i,j (n) = u i,j (n) + l i,j (n) /2

8.
Obtain the j-th trial for the i-th proto-IMF subtracting the mean value from the original data: Repeat inner steps 3 to 7 (increasing j) until the inner loop stop criterion is fulfilled, which occurs in the k-th inner iteration 10. Consider the last proto-IMF h i,k (n) as the i-th IMF: c i (n) = h i,k (n) 11. Inner loop stop criterion: During S consecutives iterations the number of upper extrema (U e ), lower extrema (L e ), and zero-crossing (Z c ) of the ongoing sequence x i,j (n) satisfy the equation |U e + L e − Z e |≤ 1 , where S is a predefined constant (usually 4 ≤ S ≤ 8)

4.
Obtain the i-th residue subtracting i-th IMF from the first trial of the ongoing sequence: Repeat outer steps 3 and 4 (increasing i) until the outer loop stop criterion is fulfilled, which occurs in the q-th outer iteration 6.
Outer loop stop criterion: Residue is a monotonic function; or IMF or residue are smaller than a predetermined value 7.
Consider the last residue as the overall residue: r(n) = r q (n).
Finally, using EMD, the sequence x(n) can be decomposed into a sum of q IMFs and a residue, according to the following expression: Considering the way of obtaining the IMFs, they satisfy the conditions so that their HSs do not have meaningless negative instantaneous frequencies. Calling Combining the instantaneous amplitude and f Hilbert Spectrum (HS) is obtained, that is: For the case of stationary sequences, it is better which is defined as:

Empirical Decomposition and the Hilbert-Huang T
In many practical sequences it is not uncomm (and their derived HS and MHS) have negativ instantaneous frequency to always be positive, the 1. It is symmetric with respect to a null local aver 2. It has the same number of zero-crossing and e A sequence fulfilling these requirements is cal to obtain meaningful Hilbert Transform, the origin several IMFs ( ). A method to achieve this goal which is based on two nested loops: An outer ite iteration to obtain ℎ , ( ), the -th trial for the -th i (n) and Sensors 2020, 20, x FOR PEER REVIEW In most applications, the resulting Hilbert Transform ̂( ) i function, obtaining the so-called analytic function, defined as: which is a time-dependent complex sequence. Its module is called in ( ) ≡ |̃( )| . An instantaneous frequency ( ) can also be obta analytic function: Combining the instantaneous amplitude and frequency in a sing Hilbert Spectrum (HS) is obtained, that is: For the case of stationary sequences, it is better to obtain the Mar which is defined as:

Empirical Decomposition and the Hilbert-Huang Transform for Spectr
In many practical sequences it is not uncommon that some Hil (and their derived HS and MHS) have negative values, which a instantaneous frequency to always be positive, the sequence must sat 1. It is symmetric with respect to a null local average value; and 2. It has the same number of zero-crossing and extrema.
A sequence fulfilling these requirements is called an Intrinsic M to obtain meaningful Hilbert Transform, the original sequence ( ) several IMFs ( ). A method to achieve this goal is named Empirica which is based on two nested loops: An outer iteration to obtain th iteration to obtain ℎ , ( ), the -th trial for the -th IMF. The EMD process) can be summarized in the following steps [42]: 1. Consider original data as the first residue: 1 ( ) = ( ) i (n) the instantaneous amplitude and frequency of the i-th IMF, the Hilbert-Huang Transform (HHT) is defined as the ensemble of the Hilbert Spectrum (HS) of each IMF and the residue, that is, q + 1 vectors each of them of two functions: HHT [x(n)] = Sensors 2020, 20, x FOR PEER REVIEW In most applications, the resulting Hilbert Transform ̂( ) is combined with the function, obtaining the so-called analytic function, defined as: which is a time-dependent complex sequence. Its module is called instantaneous amplitude ( ) ≡ |̃( )| . An instantaneous frequency ( ) can also be obtained using the phase analytic function: Combining the instantaneous amplitude and frequency in a single vector of two funct Hilbert Spectrum (HS) is obtained, that is: For the case of stationary sequences, it is better to obtain the Marginal Hilbert Spectrum which is defined as:

Empirical Decomposition and the Hilbert-Huang Transform for Spectral Analysis
In many practical sequences it is not uncommon that some Hilbert instantaneous freq (and their derived HS and MHS) have negative values, which are meaningless [40]. instantaneous frequency to always be positive, the sequence must satisfy two conditions [41 1. It is symmetric with respect to a null local average value; and In most applications, the resulting Hilbert Transform ̂( ) is combined with the original function, obtaining the so-called analytic function, defined as: which is a time-dependent complex sequence. Its module is called instantaneous amplitude, that is, ( ) ≡ |̃( )| . An instantaneous frequency ( ) can also be obtained using the phase of the analytic function: Combining the instantaneous amplitude and frequency in a single vector of two functions, the Hilbert Spectrum (HS) is obtained, that is: For the case of stationary sequences, it is better to obtain the Marginal Hilbert Spectrum (MHS), which is defined as: (14)

Empirical Decomposition and the Hilbert-Huang Transform for Spectral Analysis
In many practical sequences it is not uncommon that some Hilbert instantaneous frequencies (and their derived HS and MHS) have negative values, which are meaningless [40]. For the instantaneous frequency to always be positive, the sequence must satisfy two conditions [41]: 1. It is symmetric with respect to a null local average value; and 2. It has the same number of zero-crossing and extrema.
A sequence fulfilling these requirements is called an Intrinsic Mode Function (IMF). Therefore, Combining the instantaneous amplitude and frequency in a single vecto Hilbert Spectrum (HS) is obtained, that is: For the case of stationary sequences, it is better to obtain the Marginal Hi which is defined as:

Empirical Decomposition and the Hilbert-Huang Transform for Spectral Analys
In many practical sequences it is not uncommon that some Hilbert inst (and their derived HS and MHS) have negative values, which are mea instantaneous frequency to always be positive, the sequence must satisfy two 1. It is symmetric with respect to a null local average value; and In most applications, the resulting Hilbert Transform ̂( ) is combined with the origina function, obtaining the so-called analytic function, defined as: ̃( ) ≡ ( ) + ̂( ), (11 which is a time-dependent complex sequence. Its module is called instantaneous amplitude, that is ( ) ≡ |̃( )| . An instantaneous frequency ( ) can also be obtained using the phase of th analytic function:

(12
Combining the instantaneous amplitude and frequency in a single vector of two functions, th Hilbert Spectrum (HS) is obtained, that is: For the case of stationary sequences, it is better to obtain the Marginal Hilbert Spectrum (MHS) which is defined as: (14

Empirical Decomposition and the Hilbert-Huang Transform for Spectral Analysis
In many practical sequences it is not uncommon that some Hilbert instantaneous frequencie (and their derived HS and MHS) have negative values, which are meaningless [40]. For th instantaneous frequency to always be positive, the sequence must satisfy two conditions [41]: 1. It is symmetric with respect to a null local average value; and 2. It has the same number of zero-crossing and extrema.
A sequence fulfilling these requirements is called an Intrinsic Mode Function (IMF). Therefore Combining the instantaneous amplitude and frequenc Hilbert Spectrum (HS) is obtained, that is: For the case of stationary sequences, it is better to obta which is defined as:

Empirical Decomposition and the Hilbert-Huang Transform
In many practical sequences it is not uncommon that (and their derived HS and MHS) have negative values instantaneous frequency to always be positive, the sequenc 1. It is symmetric with respect to a null local average valu q (n), Sensors 2020, 20, x FOR PEER REVIEW In most applications, the resulting Hilbert Transform ̂( ) is combined function, obtaining the so-called analytic function, defined as: which is a time-dependent complex sequence. Its module is called instantaneous ( ) ≡ |̃( )| . An instantaneous frequency ( ) can also be obtained using analytic function: Combining the instantaneous amplitude and frequency in a single vector of Hilbert Spectrum (HS) is obtained, that is: For the case of stationary sequences, it is better to obtain the Marginal Hilber which is defined as:

Empirical Decomposition and the Hilbert-Huang Transform for Spectral Analysis
In many practical sequences it is not uncommon that some Hilbert instant (and their derived HS and MHS) have negative values, which are meaning instantaneous frequency to always be positive, the sequence must satisfy two con 1. It is symmetric with respect to a null local average value; and 2. It has the same number of zero-crossing and extrema.
A sequence fulfilling these requirements is called an Intrinsic Mode Functio . An instantaneous frequency analytic function: Combining the instantaneous amplitude Hilbert Spectrum (HS) is obtained, that is: For the case of stationary sequences, it is which is defined as: (

Empirical Decomposition and the Hilbert-Hu
In many practical sequences it is not un (and their derived HS and MHS) have ne instantaneous frequency to always be positiv 1. It is symmetric with respect to a null loc q+1 (n), Sensors 2020, 20, x FOR PEER REVIEW In most applications, the resulting Hilbert Transform ̂( function, obtaining the so-called analytic function, defined as: which is a time-dependent complex sequence. Its module is call ( ) ≡ |̃( )| . An instantaneous frequency ( ) can also be analytic function: Combining the instantaneous amplitude and frequency in a Hilbert Spectrum (HS) is obtained, that is: For the case of stationary sequences, it is better to obtain the which is defined as:

Empirical Decomposition and the Hilbert-Huang Transform for S
In many practical sequences it is not uncommon that some (and their derived HS and MHS) have negative values, wh instantaneous frequency to always be positive, the sequence mus 1. It is symmetric with respect to a null local average value; an 2. It has the same number of zero-crossing and extrema.
A sequence fulfilling these requirements is called an Intrins In these equations the residue is understood as the (q + 1)-th IMF. This function, which is discrete in time but continuous in frequency and amplitude, is also known as the Hilbert-Huang Spectrum (HHS).
To manage HHS as an ensemble of q + 1 pairs of vectors is often quite burdensome and leads to graphics that are not easy to read. To simplify its representation it is common to discretize the frequencies using a quantizing ∆ f step, building up a grid in the time-frequency plane. Then the amplitudes of all IMFs for every single grid spot are added, obtaining a Discrete Hilbert-Huang Spectrum (DHHS) defined by: In most applications, the resulting Hilbert Transform ̂( ) is combined with the function, obtaining the so-called analytic function, defined as: which is a time-dependent complex sequence. Its module is called instantaneous amplitud ( ) ≡ |̃( )| . An instantaneous frequency ( ) can also be obtained using the phas analytic function: Combining the instantaneous amplitude and frequency in a single vector of two func Hilbert Spectrum (HS) is obtained, that is: For the case of stationary sequences, it is better to obtain the Marginal Hilbert Spectrum which is defined as: . An instantaneous frequency analytic function: Combining the instantaneous amplitud Hilbert Spectrum (HS) is obtained, that is: For the case of stationary sequences, it i which is defined as: Sensors 2020, 20, 2912 For the case of stationary sequences, it is sometimes better to obtain the Marginal Hilbert-Huang Spectrum (MHHS), which is derived from DHHS and defined as: In most applications, the resulting Hilbert Transform ̂( ) is combined with the function, obtaining the so-called analytic function, defined as: which is a time-dependent complex sequence. Its module is called instantaneous amplitude ( ) ≡ |̃( )| . An instantaneous frequency ( ) can also be obtained using the phas analytic function: Combining the instantaneous amplitude and frequency in a single vector of two funct Hilbert Spectrum (HS) is obtained, that is: For the case of stationary sequences, it is better to obtain the Marginal Hilbert Spectrum which is defined as:

Empirical Decomposition and the Hilbert-Huang Transform for Spectral Analysis
In many practical sequences it is not uncommon that some Hilbert instantaneous fre (and their derived HS and MHS) have negative values, which are meaningless [40]. instantaneous frequency to always be positive, the sequence must satisfy two conditions [41 1. It is symmetric with respect to a null local average value; and 2. It has the same number of zero-crossing and extrema.
A sequence fulfilling these requirements is called an Intrinsic Mode Function (IMF). T to obtain meaningful Hilbert Transform, the original sequence ( ) should be first decom several IMFs ( ). A method to achieve this goal is named Empirical Mode Decompositio which is based on two nested loops: An outer iteration to obtain the -th IMF ( ); and iteration to obtain ℎ , ( ), the -th trial for the -th IMF. The EMD algorithm (also called process) can be summarized in the following steps [42] Combining the instantaneous amplitud Hilbert Spectrum (HS) is obtained, that is: For the case of stationary sequences, it i which is defined as:

Empirical Decomposition and the Hilbert-H
In many practical sequences it is not u (and their derived HS and MHS) have n instantaneous frequency to always be positiv 1. It is symmetric with respect to a null loc 2. It has the same number of zero-crossing A sequence fulfilling these requirement to obtain meaningful Hilbert Transform, the several IMFs ( ). A method to achieve thi which is based on two nested loops: An ou iteration to obtain ℎ , ( ), the -th trial for process) can be summarized in the following 1. Consider original data as the first residu 2. Make = 1 3. Obtaining the -th IMF: ( ) 1. Consider the -th residue as the firs 2. Make = 1 3. Consider the -th trial for the -th p 4. Obtain upper local extrema of , ( 5. Obtain upper envelope (joining upp 6. Obtain lower local extrema and low 7. Obtain mean value of upper and lo 8. Obtain the -th trial for the -th pr data: ℎ , ( ) = , ( ) − , ( ) 9. Repeat inner steps 3 to 7 (increasing occurs in the -th inner iteration 10. Consider the last proto-IMF ℎ , ( )

Spectral Analysis of Electricity Demand
In this subsection the spectral analysis methods previously explained will be applied to the MHLV load curve for Spain.

Fourier Transform of Electricity Demand
Let us begin applying FT to the dataset, resulting in the spectrum depicted in Figure 2a. As the sampling frequency is f s = 8760 year −1 , full spectrum contains meaningful frequencies up to f s /2 = 4380 year −1 . Regardless of the Direct Current (DC) component (corresponding to the energy demand mean value), two main harmonics can be seen: At 365 (times per year) corresponding to once per day; and at 730 corresponding to twice per day. These two harmonics reveal the main periodic behavior of the load curve featured by a daily evolution with two humps at noon and mid-afternoon. Let us begin applying FT to the dataset, resulting in the spectrum depicted in Figure 2a. As the sampling frequency is = 8760 year −1 , full spectrum contains meaningful frequencies up to /2 = 4380 year −1 . Regardless of the Direct Current (DC) component (corresponding to the energy demand mean value), two main harmonics can be seen: At 365 (times per year) corresponding to once per day; and at 730 corresponding to twice per day. These two harmonics reveal the main periodic behavior of the load curve featured by a daily evolution with two humps at noon and mid-afternoon.
Zooming in on the spectrum to discard higher frequencies (Figure 2b), a third harmonic stands out at the frequency of 52 (times per year), revealing a weekly periodicity characterized by a higher demand on labor days and lower consumption on weekends. A much closer zoom centered on the low frequency spectrum ( Figure 2c) shows a harmonic peak at frequency 2 (twice per year) and a lower and not very well defined peak at 1 (once per year), corresponding to a yearly periodicity with two seasonal peaks of high demand (due to heating in winter and air conditioning in summer). It has to be noted that this graph has a frequency resolution which can be computed by:

Short-Time Fourier Transform of Electricity Demand
Trying to discover transient behaviors, an STFT with time resolution = 1 has been applied to the dataset's sequence, obtaining the result depicted in Figure 3a. It can be seen that the spectrogram is quite stationary, that is, it has approximately the same values along time, although a small decrease during summer can be noted (mainly in higher frequencies). The most significant harmonics (dark red in the graph) appear in the low-frequency band, and a zoom-in of the frequency axis is shown in Figure 3b, where several harmonics clearly appear (around 50, 400, and 750 −1 ). However, their values are not precisely determined as the corresponding frequency resolution is, according to Equation (5), Δ = 52 −1 . Zooming in on the spectrum to discard higher frequencies (Figure 2b), a third harmonic stands out at the frequency of 52 (times per year), revealing a weekly periodicity characterized by a higher demand on labor days and lower consumption on weekends. A much closer zoom centered on the low frequency spectrum (Figure 2c) shows a harmonic peak at frequency 2 (twice per year) and a lower and not very well defined peak at 1 (once per year), corresponding to a yearly periodicity with two seasonal peaks of high demand (due to heating in winter and air conditioning in summer). It has to be noted that this graph has a frequency resolution which can be computed by:

Short-Time Fourier Transform of Electricity Demand
Trying to discover transient behaviors, an STFT with time resolution ∆t = 1 week has been applied to the dataset's sequence, obtaining the result depicted in Figure 3a. It can be seen that the spectrogram is quite stationary, that is, it has approximately the same values along time, although a small decrease during summer can be noted (mainly in higher frequencies). The most significant harmonics (dark red in the graph) appear in the low-frequency band, and a zoom-in of the frequency axis is shown in Figure 3b, where several harmonics clearly appear (around 50, 400, and 750 years −1 ). However, their values are not precisely determined as the corresponding frequency resolution is, according to Equation (5), ∆ f = 52 years −1 . spectrum discarding higher frequencies; (c) low frequency spectrum.

Short-Time Fourier Transform of Electricity Demand
Trying to discover transient behaviors, an STFT with time resolution = 1 has been applied to the dataset's sequence, obtaining the result depicted in Figure 3a. It can be seen that the spectrogram is quite stationary, that is, it has approximately the same values along time, although a small decrease during summer can be noted (mainly in higher frequencies). The most significant harmonics (dark red in the graph) appear in the low-frequency band, and a zoom-in of the frequency axis is shown in Figure 3b, where several harmonics clearly appear (around 50, 400, and 750 ). However, their values are not precisely determined as the corresponding frequency resolution is, according to Equation (5), Δ = 52 . To increase frequency resolution, the width of the window function has to be enlarged and, consequently, the time resolution decreases. Figure 3c depicts the result when a ∆t = 0.5 years is applied, resulting in a frequency resolution of ∆ f = 2 years −1 .

Wavelet Transform of Electricity Demand
To overcome the limitations of STFT regarding the time-frequency resolution, a Morlet WT is applied to the sequence according to Equation (6), obtaining the scalogram shown in Figure 4a. For an easier reading, the vertical axis has been converted from scale s to frequency. Again, the electricity demand frequency representation shows a quite stationary behavior, although a small decrease during summer can be noted in the highest frequencies. Several harmonics clearly appear in the plot, centered around 50, 400, and 750 years −1 . Figure 4b presents a detail of the scalogram for the lower frequencies.
Sensors 2020, 20, x FOR PEER REVIEW 8 of 25 To increase frequency resolution, the width of the window function has to be enlarged and, consequently, the time resolution decreases. Figure 3c depicts the result when a = 0.5 is applied, resulting in a frequency resolution of Δ = 2 .

Wavelet Transform of Electricity Demand
To overcome the limitations of STFT regarding the time-frequency resolution, a Morlet WT is applied to the sequence according to Equation (6), obtaining the scalogram shown in Figure 4a. For an easier reading, the vertical axis has been converted from scale to frequency. Again, the electricity demand frequency representation shows a quite stationary behavior, although a small decrease during summer can be noted in the highest frequencies. Several harmonics clearly appear in the plot, centered around 50, 400, and 750 . Figure 4b presents a detail of the scalogram for the lower frequencies.

Marginal Wavelet Transform of Electricity Demand
If the sequence analysis considers only its stationary behavior, it is better to obtain the marginal spectrum applying Equation (8). The result is shown in Figure 5 using different graphic scales and with different frequency resolution. The accumulated over time spectral amplitude (vertical axis) versus frequency (horizontal axis) is plotted. Depicted in red are those harmonics already identified in Figure 2, while the values in green reveal new harmonics.

Marginal Wavelet Transform of Electricity Demand
If the sequence analysis considers only its stationary behavior, it is better to obtain the marginal spectrum applying Equation (8). The result is shown in Figure 5 using different graphic scales and with different frequency resolution. The accumulated over time spectral amplitude (vertical axis) versus frequency (horizontal axis) is plotted. Depicted in red are those harmonics already identified in Figure 2, while the values in green reveal new harmonics.
If the sequence analysis considers only its stationary behavior, it is better to obtain the marginal spectrum applying Equation (8). The result is shown in Figure 5 using different graphic scales and with different frequency resolution. The accumulated over time spectral amplitude (vertical axis) versus frequency (horizontal axis) is plotted. Depicted in red are those harmonics already identified in Figure 2, while the values in green reveal new harmonics.

Hilbert Spectrum of Electricity Demand
Now that the more conventional spectral analysis techniques have been employed, let us focus on the application of Hilbert Transform to our dataset. In Figure 6a, Hilbert Spectrum is depicted in accordance with the definitions in Equation (13).
Sensors 2020, 20, x FOR PEER REVIEW 9 of 25 Now that the more conventional spectral analysis techniques have been employed, let us focus on the application of Hilbert Transform to our dataset. In Figure 6a, Hilbert Spectrum is depicted in accordance with the definitions in Equation (13). The plot represents the instantaneous frequency (vertical axis) versus time (horizontal axis), while the instantaneous amplitude corresponding to each time is also color-coded. As there is an instantaneous amplitude-frequency pair at every hour, the total number of points (pairs) in the full-scale graph is quite high ( = 29,183). Therefore, the figure is easier to understand if the pairs are drawn as scattered points, that is, not forming a line.
On the other hand, Figure 6b shows the Marginal Hilbert Spectrum as it is defined in Equation (14). The accumulated-over-time amplitude (vertical axis) versus instantaneous frequency (horizontal axis) is plotted. In both figures, meaningless negative frequencies clearly appear, making a physical interpretation of these spectra difficult.

Hilbert-Huang Spectrum of Electricity Demand
To tackle the Hilbert Spectrum's lack of meaning, the original electricity demand sequence is split into several Intrinsic Mode Functions (IMF) using the Empirical Mode Decomposition (EMD) method described in Section 2.4. The result of this process is a set of 12 IMFs and a residue, which are depicted in Figure 7 using a monthly timescale. The corresponding original signal is shown in the first panel and in Figure 1a. The last IMFs (those with higher indexes) contain the electricity demand slow-rate evolution. Therefore, IMFs 12, 11, 10, and 9 show an approximate period of 30, 12, 6, and 2 months, respectively. It can be seen that, for instance, the higher demands in winter and summer are neatly uncovered by the peaks in IMF 10 (6-month period). The plot represents the instantaneous frequency (vertical axis) versus time (horizontal axis), while the instantaneous amplitude corresponding to each time is also color-coded. As there is an instantaneous amplitude-frequency pair at every hour, the total number of points (pairs) in the full-scale graph is quite high (N = 29,183) . Therefore, the figure is easier to understand if the pairs are drawn as scattered points, that is, not forming a line.
On the other hand, Figure 6b shows the Marginal Hilbert Spectrum as it is defined in Equation (14). The accumulated-over-time amplitude (vertical axis) versus instantaneous frequency (horizontal axis) is plotted. In both figures, meaningless negative frequencies clearly appear, making a physical interpretation of these spectra difficult.

Hilbert-Huang Spectrum of Electricity Demand
To tackle the Hilbert Spectrum's lack of meaning, the original electricity demand sequence is split into several Intrinsic Mode Functions (IMF) using the Empirical Mode Decomposition (EMD) method described in Section 2.4. The result of this process is a set of 12 IMFs and a residue, which are depicted in Figure 7 using a monthly timescale. The corresponding original signal is shown in the first panel and in Figure 1a. The last IMFs (those with higher indexes) contain the electricity demand slow-rate evolution. Therefore, IMFs 12, 11, 10, and 9 show an approximate period of 30, 12, 6, and 2 months, respectively. It can be seen that, for instance, the higher demands in winter and summer are neatly uncovered by the peaks in IMF 10 (6-month period). method described in Section 2.4. The result of this process is a set of 12 IMFs and a residue, which are depicted in Figure 7 using a monthly timescale. The corresponding original signal is shown in the first panel and in Figure 1a. The last IMFs (those with higher indexes) contain the electricity demand slow-rate evolution. Therefore, IMFs 12, 11, 10, and 9 show an approximate period of 30, 12, 6, and 2 months, respectively. It can be seen that, for instance, the higher demands in winter and summer are neatly uncovered by the peaks in IMF 10 (6-month period). For a better understanding of the IMFs behavior, they are also represented in a weekly timescale. Figure 8 depicts a 10-week sequence corresponding to the time-domain representation of For a better understanding of the IMFs behavior, they are also represented in a weekly timescale.  Finally, Figure 9 depicts a 15-day sequence corresponding to the time-domain representation of the original signal shown in the first panel and in Figure 1c. It reveals approximate periods of 4, 2, 1, and 0.5 days for IMFs 4, 3, 2, and 1, respectively. Again, connections between original signal and IMFs can be easily established. For instance, IMF 1 (0.5-day period) presents very low amplitude corresponding to two low energy demand periods during weekend days 16 and 23, both clearly visible in Figure 1b. These weekend days can also be detected by low amplitude in IMF 2 (1-day period).   Figure 1c. It reveals approximate periods of 4, 2, 1, and 0.5 days for IMFs 4, 3, 2, and 1, respectively. Again, connections between original signal and IMFs can be easily established. For instance, IMF 1 (0.5-day period) presents very low amplitude corresponding to two low energy demand periods during weekend days 16 and 23, both clearly visible in Figure 1b. These weekend days can also be detected by low amplitude in IMF 2 (1-day period).
After the EMD process, every resulting IMF does satisfy the conditions expressed in Section 2.4 for having a meaningful Hilbert Transform (HT). The electricity demand Hilbert-Huang Transform (HHT) is then no more than the ensemble of the HT of each IMF, as reflected in Equation (17), which is depicted in Figure 10a. The plot represents, for each IMF, the instantaneous frequency (vertical axis) versus time (horizontal axis), while the instantaneous amplitude corresponding to each point is also color-coded. This plot can be seen as 13 overlapping Hilbert Spectra (corresponding to each of the 12 IMFs plus the residue). As the graph covers the full timescale, it is drawn as a scatter plot for a better interpretation. Finally, Figure 9 depicts a 15-day sequence corresponding to the time-domain representation of the original signal shown in the first panel and in Figure 1c. It reveals approximate periods of 4, 2, 1, and 0.5 days for IMFs 4, 3, 2, and 1, respectively. Again, connections between original signal and IMFs can be easily established. For instance, IMF 1 (0.5-day period) presents very low amplitude corresponding to two low energy demand periods during weekend days 16 and 23, both clearly visible in Figure 1b. These weekend days can also be detected by low amplitude in IMF 2 (1-day period). After the EMD process, every resulting IMF does satisfy the conditions expressed in Section 2.4 for having a meaningful Hilbert Transform (HT). The electricity demand Hilbert-Huang Transform (HHT) is then no more than the ensemble of the HT of each IMF, as reflected in Equation (17), which is depicted in Figure 10a. The plot represents, for each IMF, the instantaneous frequency (vertical axis) versus time (horizontal axis), while the instantaneous amplitude corresponding to each point is also color-coded. This plot can be seen as 13 overlapping Hilbert Spectra (corresponding to each of the 12 IMFs plus the residue). As the graph covers the full timescale, it is drawn as a scatter plot for a better interpretation.
On the other hand, Figure 10b shows Marginal Hilbert-Huang Spectrum as it is defined in Equation (19). The accumulated-over-time amplitude (vertical axis) versus instantaneous frequency (horizontal axis) is plotted, with all the IMFs added. In both figures, meaningless negative frequencies have almost disappeared, presenting only some negligible values which will be discarded from now on. When the timescale is smaller, the points in the Hilbert-Huang Spectrum (HHS) can be joined in a line plot without affecting its readability. That is the case in Figure 11a where 13 lines are drawn, one per IMF. Each line represents the instantaneous frequency versus time for a single IMF, with color representing the instantaneous amplitude. It can be seen that the most significant components are in the low frequency range with two outstanding regions, around twice (IMF 1) and once (IMF 2) per day (365 and 730 times per year). It can also be noted that higher frequencies (but with low amplitudes) exist on IMF 1 during weekends.
Reading an HHS, either in its scatter or set-of-lines plot versions, is sometimes difficult due to the fact that several IMFs overlap in the graphic representation (as in Figure 10a and 11a). Therefore, it is a common practice to discretize frequencies to obtain the Discrete Hilbert-Huang Spectrum (DHHS) as it is defined in Equation (18). The result is depicted in Figure 11b where again the once and twice per day frequencies can be seen. The MHHS can now easily be extended to the whole dataset obtaining the results depicted in Figure 12a, where the red arrows on the right axis point to a stationary and well-defined 365 times per year frequency (once per day), as well as a somehow fluctuant 730 frequency (twice per day). As On the other hand, Figure 10b shows Marginal Hilbert-Huang Spectrum as it is defined in Equation (19). The accumulated-over-time amplitude (vertical axis) versus instantaneous frequency (horizontal axis) is plotted, with all the IMFs added.
In both figures, meaningless negative frequencies have almost disappeared, presenting only some negligible values which will be discarded from now on.
When the timescale is smaller, the points in the Hilbert-Huang Spectrum (HHS) can be joined in a line plot without affecting its readability. That is the case in Figure 11a where 13 lines are drawn, one per IMF. Each line represents the instantaneous frequency versus time for a single IMF, with color representing the instantaneous amplitude. It can be seen that the most significant components are in the low frequency range with two outstanding regions, around twice (IMF 1) and once (IMF 2) per day (365 and 730 times per year). It can also be noted that higher frequencies (but with low amplitudes) exist on IMF 1 during weekends. In both figures, meaningless negative frequencies have almost disappeared, presenting only some negligible values which will be discarded from now on. When the timescale is smaller, the points in the Hilbert-Huang Spectrum (HHS) can be joined in a line plot without affecting its readability. That is the case in Figure 11a where 13 lines are drawn, one per IMF. Each line represents the instantaneous frequency versus time for a single IMF, with color representing the instantaneous amplitude. It can be seen that the most significant components are in the low frequency range with two outstanding regions, around twice (IMF 1) and once (IMF 2) per day (365 and 730 times per year). It can also be noted that higher frequencies (but with low amplitudes) exist on IMF 1 during weekends.
Reading an HHS, either in its scatter or set-of-lines plot versions, is sometimes difficult due to the fact that several IMFs overlap in the graphic representation (as in Figure 10a and 11a). Therefore, it is a common practice to discretize frequencies to obtain the Discrete Hilbert-Huang Spectrum (DHHS) as it is defined in Equation (18). The result is depicted in Figure 11b where again the once and twice per day frequencies can be seen. The MHHS can now easily be extended to the whole dataset obtaining the results depicted in Figure 12a, where the red arrows on the right axis point to a stationary and well-defined 365 times per year frequency (once per day), as well as a somehow fluctuant 730 frequency (twice per day). As most of the harmonics are in the low frequency region, Figure 12 (panels b to f) shows the MHHS at Reading an HHS, either in its scatter or set-of-lines plot versions, is sometimes difficult due to the fact that several IMFs overlap in the graphic representation (as in Figures 10a and 11a). Therefore, it is a common practice to discretize frequencies to obtain the Discrete Hilbert-Huang Spectrum (DHHS) as it is defined in Equation (18). The result is depicted in Figure 11b where again the once and twice per day frequencies can be seen.
The MHHS can now easily be extended to the whole dataset obtaining the results depicted in Figure 12a, where the red arrows on the right axis point to a stationary and well-defined 365 times per year frequency (once per day), as well as a somehow fluctuant 730 frequency (twice per day). As most of the harmonics are in the low frequency region, Figure 12  If the electricity demand can be considered stationary, or if we are not interested in its transient behavior, the temporal axis of DHHS can be squeezed to obtain the Marginal Hilbert-Huang Spectrum (MHHS) according to Equation (19). Results obtained using different frequency scales are depicted in Figure 13. Besides the harmonics obtained by FFT (marked in red), the Hilbert-Huang analysis reveals new harmonics pointed in green in the graphics. The new values, in times per year, are the following: 183 (2-day period); 91 (4-day period); 12 (1-month period); 4.35 (12-week period) and a blurry region till 6 (2-month period); and 0.3 covering the full 40-month dataset. If the electricity demand can be considered stationary, or if we are not interested in its transient behavior, the temporal axis of DHHS can be squeezed to obtain the Marginal Hilbert-Huang Spectrum (MHHS) according to Equation (19). Results obtained using different frequency scales are depicted in Figure 13. Besides the harmonics obtained by FFT (marked in red), the Hilbert-Huang analysis reveals new harmonics pointed in green in the graphics. The new values, in times per year, are the following: 183 (2-day period); 91 (4-day period); 12 (1-month period); 4.35 (12-week period) and a blurry region till 6 (2-month period); and 0.3 covering the full 40-month dataset.
It has to be noted that although DHHS and its derived MHHS are quantized in frequency (also in time), the frequency step ∆ f can be arbitrary selected. To maintain the details at each frequency scale, the number of frequency quanta has been fixed to 100 in every panel of Figure 13 (also in Figure 12). On the other hand, frequency resolution in FT is fixed and depends on the number of values in the dataset (∆ f = 0.3 in our case). Then FT leads to an over-detailed spectra for full range frequencies (Figure 2a), and to an under-detailed spectra for small range frequencies (Figure 2c). However, HHT produces spectra with a selectable frequency granularity which can be adjusted to obtain similar detail levels at each frequency range, as can be seen in Figure 13. behavior, the temporal axis of DHHS can be squeezed to obtain the Marginal Hilbert-Huang Spectrum (MHHS) according to Equation (19). Results obtained using different frequency scales are depicted in Figure 13. Besides the harmonics obtained by FFT (marked in red), the Hilbert-Huang analysis reveals new harmonics pointed in green in the graphics. The new values, in times per year, are the following: 183 (2-day period); 91 (4-day period); 12 (1-month period); 4.35 (12-week period) and a blurry region till 6 (2-month period); and 0.3 covering the full 40-month dataset. It has to be noted that although DHHS and its derived MHHS are quantized in frequency (also in time), the frequency step Δ can be arbitrary selected. To maintain the details at each frequency scale, the number of frequency quanta has been fixed to 100 in every panel of Figure 13 (also in Figure 12). On the other hand, frequency resolution in FT is fixed and depends on the number of values in the dataset (Δ = 0.3 in our case). Then FT leads to an over-detailed spectra for full range frequencies (Figure 2a), and to an under-detailed spectra for small range frequencies (Figure 2c).

Spectral Analysis of the Intrinsic Mode Functions
One of the key elements of Hilbert-Huang Spectrum is the decomposition of the original sequence (the electricity demand in our case) in an ensemble of Intrinsic Mode Functions (IMFs). Although the main goal of this decomposition is to obtain a meaningful Hilbert Spectrum, the IMFs have interest on their own, as any of them convey information about the electricity demand evolution for a certain periodicity (or frequency) scale.
For this reason, it is worthy to study the spectral representation of every individual IMF. The first approach will be to apply Fourier Transform (FT), getting the results depicted in Figure 14, where the spectrum for the original sequence, the 13 IMFs, and the residue, can be compared. However, HHT produces spectra with a selectable frequency granularity which can be adjusted to obtain similar detail levels at each frequency range, as can be seen in Figure 13.

Spectral Analysis of the Intrinsic Mode Functions
One of the key elements of Hilbert-Huang Spectrum is the decomposition of the original sequence (the electricity demand in our case) in an ensemble of Intrinsic Mode Functions (IMFs). Although the main goal of this decomposition is to obtain a meaningful Hilbert Spectrum, the IMFs have interest on their own, as any of them convey information about the electricity demand evolution for a certain periodicity (or frequency) scale.
For this reason, it is worthy to study the spectral representation of every individual IMF. The first approach will be to apply Fourier Transform (FT), getting the results depicted in Figure 14, where the spectrum for the original sequence, the 13 IMFs, and the residue, can be compared. In each panel, the red circle marks the peak frequency, and the dashed green line represents the spectrum centroid defined as: where is the sampling ratio and ( ) is the DFT of the IMF obtained applying Equation (2). It can be seen that centroid is in all cases at a higher frequency than the peak, which indicates that the spectrum contains a significant right-side tail. Our second approach for the spectral characterization of IMFs will be to apply Hilbert Transform recalling that, by definition, IMFs verify the conditions to avoid meaningless negative frequencies. Then, the Marginal Hilbert Spectrum (MHS) is obtained according to Equation (14) and the results are shown in Figure 15. In each panel, the red circle marks the peak frequency, and the dashed green line represents the spectrum centroid defined as: where f s is the sampling ratio and X(k) is the DFT of the IMF obtained applying Equation (2). It can be seen that centroid is in all cases at a higher frequency than the peak, which indicates that the spectrum contains a significant right-side tail. Our second approach for the spectral characterization of IMFs will be to apply Hilbert Transform recalling that, by definition, IMFs verify the conditions to avoid meaningless negative frequencies.
Then, the Marginal Hilbert Spectrum (MHS) is obtained according to Equation (14) and the results are shown in Figure 15.
be seen that centroid is in all cases at a higher frequency than the peak, which indicates that the spectrum contains a significant right-side tail.
Our second approach for the spectral characterization of IMFs will be to apply Hilbert Transform recalling that, by definition, IMFs verify the conditions to avoid meaningless negative frequencies. Then, the Marginal Hilbert Spectrum (MHS) is obtained according to Equation (14) and the results are shown in Figure 15.  Spectra for every IMF can also be plotted as a color map with frequency (vertical axis) versus IMF index (horizontal axis) color-coding the spectrum value. Figure 16a shows the results applying DFT, while MHT results are depicted in Figure 16b. Peaks can be identified as the highest colored regions for each IMF, while the centroid is plotted as a black line. Spectra for every IMF can also be plotted as a color map with frequency (vertical axis) versus IMF index (horizontal axis) color-coding the spectrum value. Figure 16a shows the results applying DFT, while MHT results are depicted in Figure 16b. Peaks can be identified as the highest colored regions for each IMF, while the centroid is plotted as a black line. In the effort to obtain the frequency (or period) that best defines every IMF, the use of the autocorrelation will finally be explored, defined as: where ( ) denotes the electricity demand series, is its mean value, and ( ) is the Pearson autocorrelation coefficient for a certain lag (an integer number). Then, the main periodicity is defined by Δ , where is the lag corresponding to the autocorrelation first peak. Its inverse can be considered the dominant frequency according to: The autocorrelation corresponding to the first IMF is shown in Figure 17a, where the resulting 12 h main periodicity is marked with a red circle. In the effort to obtain the frequency (or period) that best defines every IMF, the use of the autocorrelation will finally be explored, defined as: where x(n) denotes the electricity demand series, µ x is its mean value, and (L) is the Pearson autocorrelation coefficient for a certain lag L (an integer number). Then, the main periodicity is defined by L p ∆t, where L p is the lag corresponding to the autocorrelation first peak. Its inverse can be considered the dominant frequency according to: The autocorrelation corresponding to the first IMF is shown in Figure 17a, where the resulting 12 h main periodicity is marked with a red circle. be considered the dominant frequency according to: The autocorrelation corresponding to the first IMF is shown in Figure 17a, where the resulting 12 h main periodicity is marked with a red circle. In this subsection, three methods to obtain the featuring frequency of every IMF have been presented. The results are summarized and compared in Figure 17b. It can be seen that all of them offer very similar results, except the DFT-obtained centroids. A more detailed analysis shows that In this subsection, three methods to obtain the featuring frequency of every IMF have been presented. The results are summarized and compared in Figure 17b. It can be seen that all of them offer very similar results, except the DFT-obtained centroids. A more detailed analysis shows that the IMF main frequency is approximately halved in every step, as is clearly indicated by the dashed red line in the plot.

Statistical Analysis of the Intrinsic Mode Functions
To get a deeper insight on IMFs, a statistical analysis of its Hilbert Transforms will now be developed. Let us consider the electricity demand sequence x(n) and its Empirical Mode Decomposition (EMD). Let us call z i (n) the resulting i-th IMF and let us obtain its Hilbert Spectrum S i (n) = HS[z i (n)] ≡ Sensors 2020, 20, x FOR PEER REVIEW 5 of 25 In most applications, the resulting Hilbert Transform ̂( ) is combined with the original function, obtaining the so-called analytic function, defined as: which is a time-dependent complex sequence. Its module is called instantaneous amplitude, that is, ( ) ≡ |̃( )| . An instantaneous frequency ( ) can also be obtained using the phase of the analytic function: Combining the instantaneous amplitude and frequency in a single vector of two functions, the Hilbert Spectrum (HS) is obtained, that is: For the case of stationary sequences, it is better to obtain the Marginal Hilbert Spectrum (MHS), which is defined as: (14)

Empirical Decomposition and the Hilbert-Huang Transform for Spectral Analysis
In many practical sequences it is not uncommon that some Hilbert instantaneous frequencies (and their derived HS and MHS) have negative values, which are meaningless [40]. For the instantaneous frequency to always be positive, the sequence must satisfy two conditions [41]: 1. It is symmetric with respect to a null local average value; and 2. It has the same number of zero-crossing and extrema.
A sequence fulfilling these requirements is called an Intrinsic Mode Function (IMF). Therefore, to obtain meaningful Hilbert Transform, the original sequence ( ) should be first decomposed in several IMFs ( ). A method to achieve this goal is named Empirical Mode Decomposition (EMD) which is based on two nested loops: An outer iteration to obtain the -th IMF ( ); and an inner iteration to obtain ℎ , ( ), the -th trial for the -th IMF. The EMD algorithm (also called a sifting process) can be summarized in the following steps [42]: 1. Consider original data as the first residue: 1 ( ) = ( ) In most applications, the resulting Hilbert Transform ̂( ) is combined with the original tion, obtaining the so-called analytic function, defined as: h is a time-dependent complex sequence. Its module is called instantaneous amplitude, that is, ≡ |̃( )| . An instantaneous frequency ( ) can also be obtained using the phase of the tic function: Combining the instantaneous amplitude and frequency in a single vector of two functions, the ert Spectrum (HS) is obtained, that is: For the case of stationary sequences, it is better to obtain the Marginal Hilbert Spectrum (MHS), h is defined as: (14)

mpirical Decomposition and the Hilbert-Huang Transform for Spectral Analysis
In many practical sequences it is not uncommon that some Hilbert instantaneous frequencies their derived HS and MHS) have negative values, which are meaningless [40]. For the ntaneous frequency to always be positive, the sequence must satisfy two conditions [41]: It is symmetric with respect to a null local average value; and It has the same number of zero-crossing and extrema.
A sequence fulfilling these requirements is called an Intrinsic Mode Function (IMF). Therefore, tain meaningful Hilbert Transform, the original sequence ( ) should be first decomposed in ral IMFs ( ). A method to achieve this goal is named Empirical Mode Decomposition (EMD) h is based on two nested loops: An outer iteration to obtain the -th IMF ( ); and an inner tion to obtain ℎ , ( ), the -th trial for the -th IMF. The EMD algorithm (also called a sifting ess) can be summarized in the following steps [42]: Consider original data as the first residue: 1 ( ) = ( ) Obtaining the -th IMF: ( ) 1. Consider the -th residue as the first trial for the -th proto-IMF: ℎ ,1 ( ) = ( ) i (n) according to Equation (13).

If the instantaneous frequency
Sensors 2020, 20, x FOR PEER REVIEW In most applications, the resulting Hilbert Transform ̂( ) is com function, obtaining the so-called analytic function, defined as: which is a time-dependent complex sequence. Its module is called instanta ( ) ≡ |̃( )| . An instantaneous frequency ( ) can also be obtained analytic function: Combining the instantaneous amplitude and frequency in a single ve Hilbert Spectrum (HS) is obtained, that is: For the case of stationary sequences, it is better to obtain the Marginal which is defined as:

Empirical Decomposition and the Hilbert-Huang Transform for Spectral Ana
In many practical sequences it is not uncommon that some Hilbert i (and their derived HS and MHS) have negative values, which are m instantaneous frequency to always be positive, the sequence must satisfy tw 1. It is symmetric with respect to a null local average value; and 2. It has the same number of zero-crossing and extrema.
A sequence fulfilling these requirements is called an Intrinsic Mode F to obtain meaningful Hilbert Transform, the original sequence ( ) shou several IMFs ( ).
A method to achieve this goal is named Empirical Mo which is based on two nested loops: An outer iteration to obtain the -th iteration to obtain ℎ , ( ), the -th trial for the -th IMF. The EMD algori process) can be summarized in the following steps [42] In most applications, the resulting Hilbert Transform ̂( ) is combined with the original function, obtaining the so-called analytic function, defined as: which is a time-dependent complex sequence. Its module is called instantaneous amplitude, that is, ( ) ≡ |̃( )| . An instantaneous frequency ( ) can also be obtained using the phase of the analytic function: Combining the instantaneous amplitude and frequency in a single vector of two functions, the Hilbert Spectrum (HS) is obtained, that is: For the case of stationary sequences, it is better to obtain the Marginal Hilbert Spectrum (MHS), which is defined as: (14)

Empirical Decomposition and the Hilbert-Huang Transform for Spectral Analysis
In many practical sequences it is not uncommon that some Hilbert instantaneous frequencies (and their derived HS and MHS) have negative values, which are meaningless [40]. For the instantaneous frequency to always be positive, the sequence must satisfy two conditions [41]: 1. It is symmetric with respect to a null local average value; and 2. It has the same number of zero-crossing and extrema.
A sequence fulfilling these requirements is called an Intrinsic Mode Function (IMF). Therefore, to obtain meaningful Hilbert Transform, the original sequence ( ) should be first decomposed in several IMFs ( ).
A method to achieve this goal is named Empirical Mode Decomposition (EMD) which is based on two nested loops: An outer iteration to obtain the -th IMF ( ); and an inner iteration to obtain ℎ , ( ), the -th trial for the -th IMF. The EMD algorithm (also called a sifting process) can be summarized in the following steps [42]: 1. Consider original data as the first residue: Sensors 2020, 20, x FOR PEER REVIEW 5 of 25 In most applications, the resulting Hilbert Transform ̂( ) is combined with the original function, obtaining the so-called analytic function, defined as: which is a time-dependent complex sequence. Its module is called instantaneous amplitude, that is, ( ) ≡ |̃( )| . An instantaneous frequency ( ) can also be obtained using the phase of the analytic function: Combining the instantaneous amplitude and frequency in a single vector of two functions, the Hilbert Spectrum (HS) is obtained, that is: For the case of stationary sequences, it is better to obtain the Marginal Hilbert Spectrum (MHS), which is defined as: (14)

Empirical Decomposition and the Hilbert-Huang Transform for Spectral Analysis
In many practical sequences it is not uncommon that some Hilbert instantaneous frequencies (and their derived HS and MHS) have negative values, which are meaningless [40]. For the instantaneous frequency to always be positive, the sequence must satisfy two conditions [41]: 1. It is symmetric with respect to a null local average value; and 2. It has the same number of zero-crossing and extrema.
A sequence fulfilling these requirements is called an Intrinsic Mode Function (IMF). Therefore, to obtain meaningful Hilbert Transform, the original sequence ( ) should be first decomposed in several IMFs ( ).
A method to achieve this goal is named Empirical Mode Decomposition (EMD) which is based on two nested loops: An outer iteration to obtain the -th IMF ( ); and an inner iteration to obtain ℎ , ( ), the -th trial for the -th IMF. The EMD algorithm (also called a sifting process) can be summarized in the following steps [42]: 1. Consider original data as the first residue: Sensors 2020, 20, x FOR PEER REVIEW In most applications, the resulting Hilbert Transform ̂( ) is combined with function, obtaining the so-called analytic function, defined as: which is a time-dependent complex sequence. Its module is called instantaneous amp ( ) ≡ |̃( )| . An instantaneous frequency ( ) can also be obtained using the analytic function: Combining the instantaneous amplitude and frequency in a single vector of two Hilbert Spectrum (HS) is obtained, that is: For the case of stationary sequences, it is better to obtain the Marginal Hilbert Spe which is defined as:

Empirical Decomposition and the Hilbert-Huang Transform for Spectral Analysis
In many practical sequences it is not uncommon that some Hilbert instantaneou (and their derived HS and MHS) have negative values, which are meaningless instantaneous frequency to always be positive, the sequence must satisfy two conditio 1. It is symmetric with respect to a null local average value; and 2. It has the same number of zero-crossing and extrema.
A sequence fulfilling these requirements is called an Intrinsic Mode Function (IM to obtain meaningful Hilbert Transform, the original sequence ( ) should be first d several IMFs ( ).
A method to achieve this goal is named Empirical Mode Decomp which is based on two nested loops: An outer iteration to obtain the -th IMF ( ); iteration to obtain ℎ , ( ), the -th trial for the -th IMF. The EMD algorithm (also c process) can be summarized in the following steps [42]: Then, its derivative is called the Probability Density Function (PDF), commonly approximated by a histogram, and defined as: In most applications, the resulting Hilbert Transform ̂( ) is combined with the original nction, obtaining the so-called analytic function, defined as: hich is a time-dependent complex sequence. Its module is called instantaneous amplitude, that is, ( ) ≡ |̃( )| . An instantaneous frequency ( ) can also be obtained using the phase of the alytic function: Combining the instantaneous amplitude and frequency in a single vector of two functions, the ilbert Spectrum (HS) is obtained, that is: For the case of stationary sequences, it is better to obtain the Marginal Hilbert Spectrum (MHS), hich is defined as: . 4

. Empirical Decomposition and the Hilbert-Huang Transform for Spectral Analysis
In many practical sequences it is not uncommon that some Hilbert instantaneous frequencies nd their derived HS and MHS) have negative values, which are meaningless [40]. For the stantaneous frequency to always be positive, the sequence must satisfy two conditions [41]: . It is symmetric with respect to a null local average value; and . It has the same number of zero-crossing and extrema. In most applications, the resulting Hilbert Transform ̂( ) is combined with the original function, obtaining the so-called analytic function, defined as: which is a time-dependent complex sequence. Its module is called instantaneous amplitude, that is, ( ) ≡ |̃( )| . An instantaneous frequency ( ) can also be obtained using the phase of the analytic function: Combining the instantaneous amplitude and frequency in a single vector of two functions, the Hilbert Spectrum (HS) is obtained, that is: For the case of stationary sequences, it is better to obtain the Marginal Hilbert Spectrum (MHS), which is defined as: (14)

Empirical Decomposition and the Hilbert-Huang Transform for Spectral Analysis
In many practical sequences it is not uncommon that some Hilbert instantaneous frequencies (and their derived HS and MHS) have negative values, which are meaningless [40]. For the instantaneous frequency to always be positive, the sequence must satisfy two conditions [41]: 1. It is symmetric with respect to a null local average value; and 2. It has the same number of zero-crossing and extrema.
A sequence fulfilling these requirements is called an Intrinsic Mode Function (IMF). Therefore, In most applications, the resulting Hilbert Transform ̂( ) is combined with the original function, obtaining the so-called analytic function, defined as: which is a time-dependent complex sequence. Its module is called instantaneous amplitude, that is, ( ) ≡ |̃( )| . An instantaneous frequency ( ) can also be obtained using the phase of the analytic function: Combining the instantaneous amplitude and frequency in a single vector of two functions, the Hilbert Spectrum (HS) is obtained, that is: For the case of stationary sequences, it is better to obtain the Marginal Hilbert Spectrum (MHS), which is defined as: (14)

Empirical Decomposition and the Hilbert-Huang Transform for Spectral Analysis
In many practical sequences it is not uncommon that some Hilbert instantaneous frequencies (and their derived HS and MHS) have negative values, which are meaningless [40]. For the instantaneous frequency to always be positive, the sequence must satisfy two conditions [41]: 1. It is symmetric with respect to a null local average value; and 2. It has the same number of zero-crossing and extrema.
A sequence fulfilling these requirements is called an Intrinsic Mode Function (IMF). Therefore, In most applications, the resulting Hilbert Transform ̂( ) is function, obtaining the so-called analytic function, defined as: which is a time-dependent complex sequence. Its module is called inst ( ) ≡ |̃( )| . An instantaneous frequency ( ) can also be obtain analytic function: Combining the instantaneous amplitude and frequency in a single Hilbert Spectrum (HS) is obtained, that is: For the case of stationary sequences, it is better to obtain the Margi which is defined as:

Empirical Decomposition and the Hilbert-Huang Transform for Spectral
In many practical sequences it is not uncommon that some Hilbe (and their derived HS and MHS) have negative values, which ar instantaneous frequency to always be positive, the sequence must satis 1. It is symmetric with respect to a null local average value; and The resulting PDF for the instantaneous frequencies of every IMF is depicted in Figure 18. In each panel, the red circle and the dashed green line marks mode and mean, respective of the instantaneous frequency. These results are similar to those obtained through the spectral analysis shown in Figure 15.
Besides instantaneous frequency, the instantaneous amplitude In most applications, the resulting Hilbert Transform ̂( ) is combined with the original function, obtaining the so-called analytic function, defined as: which is a time-dependent complex sequence. Its module is called instantaneous amplitude, that is, ( ) ≡ |̃( )| . An instantaneous frequency ( ) can also be obtained using the phase of the analytic function: Combining the instantaneous amplitude and frequency in a single vector of two functions, the Hilbert Spectrum (HS) is obtained, that is: For the case of stationary sequences, it is better to obtain the Marginal Hilbert Spectrum (MHS), which is defined as: i (horizontal axis) versus instantaneous frequency (vertical axis).
The resulting PDF for the instantaneous frequencies of every IMF is depicted in Figure 18. In each panel, the red circle and the dashed green line marks mode and mean, respective of the instantaneous frequency. These results are similar to those obtained through the spectral analysis shown in Figure 15. Besides instantaneous frequency, the instantaneous amplitude ( ) can also be analyzed as a random variable, and its probability density function can be derived in a similar way. The resulting amplitude PDFs (vertical axis) versus instantaneous amplitude (horizontal axis) are depicted in Figure 19 (blue line). For the sake of completeness, the same figure also depicts in a rotated way (orange line), the frequency PDFs (horizontal axis) versus instantaneous frequency (vertical axis). Probability Density Functions for every IMF can also be plotted as a color map with instantaneous amplitude (vertical axis) versus IMF index (horizontal axis), color-coding the PDF value (normalized in the 0-1 range), as is shown in Figure 20a. A similar plot for instantaneous frequency PDF is depicted in Figure 20b. Peaks (mode) can be identified as the highest colored regions for each IMF, while mean value is plotted as a black line. Previous statistical analyses consider instantaneous amplitude and frequency as two independent variables, but this is not really the case. Extending the PDF concept to the full Hilbert Spectrum considered as a 2-dimensional random process ( ) = ( ), ( ) , a joint PDF can be obtained, and the results for each IMF are shown in Figure 21. They are plotted color-coding the logarithm of the joint PDF (normalized in the 0-1 range) as a function of the instantaneous amplitude (horizontal axis) and frequency (vertical axis). The mean value , is marked with a black cross.
It can be seen that some IMFs present a very well defined frequency, for example, IMF 3 at 365 times per year, IMF 6 at 52, or IMF 11 at once per year. Some others, for instance IMF 2 and 5, show this defined frequency at 365 and 52, plus a secondary lower-amplitude higher frequency region at around 730 and 91, respectively. Probability Density Functions for every IMF can also be plotted as a color map with instantaneous amplitude (vertical axis) versus IMF index (horizontal axis), color-coding the PDF value (normalized in the 0-1 range), as is shown in Figure 20a. A similar plot for instantaneous frequency PDF is depicted in Figure 20b. Peaks (mode) can be identified as the highest colored regions for each IMF, while mean value is plotted as a black line. Probability Density Functions for every IMF can also be plotted as a color map with instantaneous amplitude (vertical axis) versus IMF index (horizontal axis), color-coding the PDF value (normalized in the 0-1 range), as is shown in Figure 20a. A similar plot for instantaneous frequency PDF is depicted in Figure 20b. Peaks (mode) can be identified as the highest colored regions for each IMF, while mean value is plotted as a black line. Previous statistical analyses consider instantaneous amplitude and frequency as two independent variables, but this is not really the case. Extending the PDF concept to the full Hilbert Spectrum considered as a 2-dimensional random process ( ) = ( ), ( ) , a joint PDF can be obtained, and the results for each IMF are shown in Figure 21. They are plotted color-coding the logarithm of the joint PDF (normalized in the 0-1 range) as a function of the instantaneous amplitude (horizontal axis) and frequency (vertical axis). The mean value , is marked with a black cross.
It can be seen that some IMFs present a very well defined frequency, for example, IMF 3 at 365 times per year, IMF 6 at 52, or IMF 11 at once per year. Some others, for instance IMF 2 and 5, show this defined frequency at 365 and 52, plus a secondary lower-amplitude higher frequency region at around 730 and 91, respectively. Combining the instantaneous amplitude and frequency in a single vector of tw Hilbert Spectrum (HS) is obtained, that is: For the case of stationary sequences, it is better to obtain the Marginal Hilbert Sp which is defined as: −1 (n) , a joint PDF can be obtained, and the results for each IMF are shown in Figure 21. They are plotted color-coding the logarithm of the joint PDF (normalized in the 0-1 range) as a function of the instantaneous amplitude (horizontal axis) and frequency (vertical axis). The mean value µ Sensors 2020, 20, x FOR PEER REVIEW In most applications, the resulting Hilbert Transform ̂( ) is function, obtaining the so-called analytic function, defined as: which is a time-dependent complex sequence. Its module is called in ( ) ≡ |̃( )| . An instantaneous frequency ( ) can also be obta analytic function: Combining the instantaneous amplitude and frequency in a sing

, µ
Sensors 2020, 20, x FOR PEER REVIEW In most applications, the resulting Hilbert Transform ̂( ) is combined with the or function, obtaining the so-called analytic function, defined as: which is a time-dependent complex sequence. Its module is called instantaneous amplitude, t ( ) ≡ |̃( )| . An instantaneous frequency ( ) can also be obtained using the phase analytic function: Combining the instantaneous amplitude and frequency in a single vector of two function Hilbert Spectrum (HS) is obtained, that is: is marked with a black cross.

Discussion
This section will provide an interpretation of the previous results and a comparison of HHT with more conventional spectral analysis techniques, such as Fourier or Wavelet Transforms. We will focus on three aspects: The ability of HHT to provide more physically accurate frequency detection; its application to abnormal behavior; and its power as an information compression technique.

Reading HHT Spectrum
Comparing spectral analysis by using HHT, as is depicted in Figure 13, and results obtained using conventional FFT, as in Figure 2, clearly shows that HHT provides more clear and meaningful spectra. HHT is also capable of discovering new underlying frequencies, such as 91 and 183 times per year.
Moreover, FFT-based spectrum frequency resolution is limited by the sampling frequency and the number of values, as it is described in Equation (20), while HHT-based spectra can arbitrarily choose frequency resolution by freely defining the quantizing Δ step in Equation (18). This HHT feature is especially powerful to detect low frequency components (long-term periods) of electricity demand.
Wavelet Transform (WT), as it is shown in Figure 5, also enhances FFT low-frequency resolution. However, compared to HHT, WT is not as powerful at discovering new frequencies, and they are less accurate [43].
Therefore, HHT-based spectral analysis results are clearer, more intuitive, and more physically meaningful than FFT and even WT. Other authors have also obtained similar results studying electrical users [44], and other electrical or mechanical sequences [45,46].

Detecting Abnormal Energy Demand Behavior
Although full research on abnormal detection is out of the scope of the present paper, some of the findings in our research also have significant implications in detecting electricity demand abnormalities.
The full sequence decomposition using EMD provides an ensemble of IMFs of different periodicities (or frequencies). Therefore, by choosing the proper IMF, the abnormal behavior at the IMF periodicity scale can be detected. For instance, IMF 6 has a periodicity of about one week (see Figure 15); thus, examining this IMF, it would be possible to detect weeks with abnormal electricity demand. Similarly, as IMF 3 has a periodicity of one day, abnormal daily consumptions would be detected based on its analysis. EMD thus provides a tool for searching for abnormalities at a determined timescale.
On the other hand, the Hilbert Spectrum of each IMF provides instantaneous amplitude and frequency which can be labeled as normal or abnormal using different techniques. As a first approach, amplitude (or frequency or joint) Probability Density Function (PDF) can be used as the likelihood of normal behavior. Detecting low likelihood is therefore equivalent to detecting abnormal demand. Figure 22 shows the instantaneous amplitude logarithmic likelihood for every IMF. It can be seen that some IMFs present a very well defined frequency, for example, IMF 3 at 365 times per year, IMF 6 at 52, or IMF 11 at once per year. Some others, for instance IMF 2 and 5, show this defined frequency at 365 and 52, plus a secondary lower-amplitude higher frequency region at around 730 and 91, respectively.

Discussion
This section will provide an interpretation of the previous results and a comparison of HHT with more conventional spectral analysis techniques, such as Fourier or Wavelet Transforms. We will focus on three aspects: The ability of HHT to provide more physically accurate frequency detection; its application to abnormal behavior; and its power as an information compression technique.

Reading HHT Spectrum
Comparing spectral analysis by using HHT, as is depicted in Figure 13, and results obtained using conventional FFT, as in Figure 2, clearly shows that HHT provides more clear and meaningful spectra. HHT is also capable of discovering new underlying frequencies, such as 91 and 183 times per year.
Moreover, FFT-based spectrum frequency resolution is limited by the sampling frequency and the number of values, as it is described in Equation (20), while HHT-based spectra can arbitrarily choose frequency resolution by freely defining the quantizing ∆ f step in Equation (18). This HHT feature is especially powerful to detect low frequency components (long-term periods) of electricity demand.
Wavelet Transform (WT), as it is shown in Figure 5, also enhances FFT low-frequency resolution. However, compared to HHT, WT is not as powerful at discovering new frequencies, and they are less accurate [43].
Therefore, HHT-based spectral analysis results are clearer, more intuitive, and more physically meaningful than FFT and even WT. Other authors have also obtained similar results studying electrical users [44], and other electrical or mechanical sequences [45,46].

Detecting Abnormal Energy Demand Behavior
Although full research on abnormal detection is out of the scope of the present paper, some of the findings in our research also have significant implications in detecting electricity demand abnormalities.
The full sequence decomposition using EMD provides an ensemble of IMFs of different periodicities (or frequencies). Therefore, by choosing the proper IMF, the abnormal behavior at the IMF periodicity scale can be detected. For instance, IMF 6 has a periodicity of about one week (see Figure 15); thus, examining this IMF, it would be possible to detect weeks with abnormal electricity demand. Similarly, as IMF 3 has a periodicity of one day, abnormal daily consumptions would be detected based on its analysis. EMD thus provides a tool for searching for abnormalities at a determined timescale.
On the other hand, the Hilbert Spectrum of each IMF provides instantaneous amplitude and frequency which can be labeled as normal or abnormal using different techniques. As a first approach, amplitude (or frequency or joint) Probability Density Function (PDF) can be used as the likelihood of normal behavior. Detecting low likelihood is therefore equivalent to detecting abnormal demand. Figure 22 shows the instantaneous amplitude logarithmic likelihood for every IMF. Considering, for instance, IMF 6, a very low likelihood can be appreciated at approximately time 2.0 (years), corresponding to New Year's week. As this IMF has a weekly periodicity, it should correspond to a week with abnormal electricity demand. Let us take a closer look, just as an example, at this IMF and the instantaneous amplitude logarithmic likelihood as it is shown in Figure 23a.  Considering, for instance, IMF 6, a very low likelihood can be appreciated at approximately time 2.0 (years), corresponding to New Year's week. As this IMF has a weekly periodicity, it should correspond to a week with abnormal electricity demand. Let us take a closer look, just as an example, at this IMF and the instantaneous amplitude logarithmic likelihood as it is shown in Figure 23a. Considering, for instance, IMF 6, a very low likelihood can be appreciated at approximately time 2.0 (years), corresponding to New Year's week. As this IMF has a weekly periodicity, it should correspond to a week with abnormal electricity demand. Let us take a closer look, just as an example, at this IMF and the instantaneous amplitude logarithmic likelihood as it is shown in Figure 23a.  This low likelihood value is due to an abnormally high (and wide) value of the instantaneous amplitude of the IMF 6 Hilbert Spectrum as shown in Figure 23b, where the dashed gray line indicates the low likelihood time. The IMF 6 sequence also displays high and wide values at the abnormal time, as is depicted in Figure 23c. Finally, Figure 23d shows the evolution of electricity demand (in blue) and weekly demand (in orange), where it can be seen that, at the abnormal demand time, there is a weekly demand in "V", while one year earlier or later, the weekly demand mainly looks like a "W".
The abnormal behavior of electricity demand detected analyzing IMF 6 is shown in Figure 24, where a normal week is displayed in panel (a), while three New Year's weeks are shown in the remaining panels. A normal week (a) usually has two days of lesser demand: Saturday and Sunday. The New Year's week of 2017 (b) does not very much differ from this pattern. However, 2018 New Year's week (c) has three days (not just two) of lesser demand. On the other hand, 2019 New Year's week (d) also has an abnormal pattern, but with 4 days of lesser demand. This longer abnormality is not detected by IMF 6 but, instead, it is caught by IMF 7, probably due to its longer periodicity. This low likelihood value is due to an abnormally high (and wide) value of the instantaneous amplitude of the IMF 6 Hilbert Spectrum as shown in Figure 23b, where the dashed gray line indicates the low likelihood time. The IMF 6 sequence also displays high and wide values at the abnormal time, as is depicted in Figure 23c. Finally, Figure 23d shows the evolution of electricity demand (in blue) and weekly demand (in orange), where it can be seen that, at the abnormal demand time, there is a weekly demand in "V", while one year earlier or later, the weekly demand mainly looks like a "W".
The abnormal behavior of electricity demand detected analyzing IMF 6 is shown in Figure 24, where a normal week is displayed in panel (a), while three New Year's weeks are shown in the remaining panels. A normal week (a) usually has two days of lesser demand: Saturday and Sunday. The New Year's week of 2017 (b) does not very much differ from this pattern. However, 2018 New Year's week (c) has three days (not just two) of lesser demand. On the other hand, 2019 New Year's week (d) also has an abnormal pattern, but with 4 days of lesser demand. This longer abnormality is not detected by IMF 6 but, instead, it is caught by IMF 7, probably due to its longer periodicity. Abnormal demand detection is an important issue in smart metering when an algorithm is incorporated into the meters to detect abnormal electricity consumption as a result of illegal human intervention [12], or when external events produce anomalies in large systems' energy demand [47,48] or customer-related consumption [49].
Some frequency domain-based methods have been proposed for anomaly detection in power load curves either using Fourier Transform [50] or Wavelet Transform [51]. However, to the best of our knowledge, no solution has been proposed using Hilbert-Huang Transform. As has been revealed throughout the paper, HHT-based anomaly detection provides a multiscale perspective, making possible to detect short-, medium-, and long-term anomalies.
A related issue is load curve characterization in which spectral methods have also been proposed using Fourier [27] or Wavelet Transforms [52,53]. No HHT-based solution for load or client profiling has been found in the technical literature. Again, the multiscale perspective provided for Empirical Mode Decomposition deserves to be further explored.

Electricity Demand Sequence Compression
One of the applications of spectral analysis is its capability of compressing the information required to characterize a signal or sequence. Let us consider a sequence ( ) of values representing, for instance, the hourly electricity demand. Its Fourier Transform (FT) also has harmonics. In fact, only /2 harmonics are required as the resulting spectrum is symmetric, but each harmonic is a complex number defined by 2 values, so values are still required to describe the spectrum. If instead of values, coefficients are kept and − are discarded, the inverse FT recovers an approximation ( ). Therefore, by reducing the amount of information representing Abnormal demand detection is an important issue in smart metering when an algorithm is incorporated into the meters to detect abnormal electricity consumption as a result of illegal human intervention [12], or when external events produce anomalies in large systems' energy demand [47,48] or customer-related consumption [49].
Some frequency domain-based methods have been proposed for anomaly detection in power load curves either using Fourier Transform [50] or Wavelet Transform [51]. However, to the best of our knowledge, no solution has been proposed using Hilbert-Huang Transform. As has been revealed throughout the paper, HHT-based anomaly detection provides a multiscale perspective, making possible to detect short-, medium-, and long-term anomalies.
A related issue is load curve characterization in which spectral methods have also been proposed using Fourier [27] or Wavelet Transforms [52,53]. No HHT-based solution for load or client profiling has been found in the technical literature. Again, the multiscale perspective provided for Empirical Mode Decomposition deserves to be further explored.

Electricity Demand Sequence Compression
One of the applications of spectral analysis is its capability of compressing the information required to characterize a signal or sequence. Let us consider a sequence x(n) of N values representing, for instance, the hourly electricity demand. Its Fourier Transform (FT) also has N harmonics. In fact, only N/2 harmonics are required as the resulting spectrum is symmetric, but each harmonic is a complex number defined by 2 values, so N values are still required to describe the spectrum. If instead of N values, M coefficients are kept and N − M are discarded, the inverse FT recovers an approximation x(n). Therefore, by reducing the amount of information representing the sequence, an error is introduced. A common measure of the overall compressing error is the Root Mean Square Error (RMSE), defined as: An alternative for the spectral representation of a sequence is the use of EMD. By decomposing the sequence x(n), several IMFs are obtained, each of them with N values. However, recalling the procedure to obtain IMFs described in Section 2.4, only the extrema points are required to fully define an IMF as the remaining points are obtained by a standard cubic spline interpolation. As the number of extrema can be considerably lower than N, a reduction in the number of coefficients required to represent x(n) can be expected.
Applying these ideas to our electricity demand dataset, the results shown in Figure 25 are obtained. In panel (a) the number of coefficients required to fully describe each IMF is displayed. This number is approximately halved in every IMF, as is clearly indicated by the dashed red line in the plot. If a compression is required, only extrema representing lower frequencies IMFs are kept.
Sensors 2020, 20, x FOR PEER REVIEW 20 of 25 the sequence, an error is introduced. A common measure of the overall compressing error is the Root Mean Square Error (RMSE), defined as: An alternative for the spectral representation of a sequence is the use of EMD. By decomposing the sequence ( ), several IMFs are obtained, each of them with values. However, recalling the procedure to obtain IMFs described in Section 2.4, only the extrema points are required to fully define an IMF as the remaining points are obtained by a standard cubic spline interpolation. As the number of extrema can be considerably lower than , a reduction in the number of coefficients required to represent ( ) can be expected.
Applying these ideas to our electricity demand dataset, the results shown in Figure 25 are obtained. In panel (a) the number of coefficients required to fully describe each IMF is displayed. This number is approximately halved in every IMF, as is clearly indicated by the dashed red line in the plot. If a compression is required, only extrema representing lower frequencies IMFs are kept. The cumulative number of values required when several IMFs are employed is displayed in Figure 25b. These results, expressed as a percentage of the total number of values in the original sequence ( ), clearly reveal that even when using the full set of IMFs, the information required is 34.85% lower than .
From previous results, it is clear that the electricity demand sequence can be compressed either using a Fourier Decomposition (FD) in its harmonics, or an Empirical Mode Decomposition (EMD) in its IMFs. Comparing the error introduced by both procedures, the results depicted in Figure 26 are obtained. For a very high compression ratio (higher than 90%) both methods offer equivalent performance (FD slightly better). In the medium-range compression ratio, FD obtains smaller errors than EMD. But for zero-error compression EMD clearly outperforms FD. The cumulative number of values required when several IMFs are employed is displayed in Figure 25b. These results, expressed as a percentage of the total number of values in the original sequence (N), clearly reveal that even when using the full set of IMFs, the information required is 34.85% lower than N.
From previous results, it is clear that the electricity demand sequence can be compressed either using a Fourier Decomposition (FD) in its harmonics, or an Empirical Mode Decomposition (EMD) in its IMFs. Comparing the error introduced by both procedures, the results depicted in Figure 26 are obtained. For a very high compression ratio (higher than 90%) both methods offer equivalent performance (FD slightly better). In the medium-range compression ratio, FD obtains smaller errors than EMD. But for zero-error compression EMD clearly outperforms FD.
As a final discussion, let us compare the analysis methods discussed in this section. The first and main goal of any spectral analysis should be identifying a signal's underlying frequencies (or periodicities). Regarding electricity demand, the frequency discovering capabilities of each method are summarized in Table 1.
It can be seen that the Hilbert-Huang Transform discovers more physically meaningful frequencies with an excellent time-frequency resolution. Additionally, a qualitative assessment of the spectral analysis methods is summarized in Table 2.