Theory of the Hilbert Spectrum
Abstract
This paper is a contribution to the old problem of representing a signal in the coordinates of time and frequency. As the starting point, we abandon Gabor’s complex extension and reevaluate fundamental principles of timefrequency analysis. We provide a multicomponent model of a signal that enables rigorous definition of instantaneous frequency on a percomponent basis. Within our framework, we have shifted all uncertainty of the latent signal to its quadrature. In this approach, uncertainty is not a fundamental limitation of analysis, but rather a manifestation of the limited view of the observer. With the appropriate assumptions made on the signal model, the instantaneous amplitude and instantaneous frequency can be obtained exactly, hence exact representation of a signal in the coordinates of time and frequency can be achieved. However, uncertainty now arises in obtaining the correct assumptions, i.e. how to correctly choose the quadrature of the components.
keywords:
Hilbert Space, Signal Analysis, Instantaneous Frequency, Hilbert Spectrum, Latent Signal Analysis, AM–FM Modelingparright #1
url]http://www.HilbertSpectrum.com
Part I. Latent Signal Analysis and the Analytic Signal
\pdfstringdef\BKM@titlePart I. Latent Signal Analysis and the Analytic Signal
Interest in the proper definition of instantaneous frequency first arose with the advent of frequency modulation for radio transmission in the 1920’s …. [T]he “analytic signal procedure” devised by Gabor results in a complex signal that has spectrum identical to that of the real signal for positive frequencies and zero for the negative frequencies …. However, instantaneous frequency is a primitive concept and not a question of mere mathematical definition …. One should keep an open mind regarding the proper definition of the complex signal, that is, the appropriate way to define phase, amplitude, and instantaneous frequency. Probably the last word on the subject has not yet been said.—Cohen cohen1995time ()
In the first part of this paper, we present the Latent Signal Analysis (LSA) problem as a recasting of the classic complex extension problem. Almost universally, the solution approach has been to use the Hilbert Transform (HT) to construct Gabor’s Analytic Signal (AS). This approach relies on Harmonic Correspondence (HC), which may lead to incorrect Instantaneous Amplitude (IA) and Instantaneous Frequency (IF) parameters. We show that by relaxing HC, the resulting complex extension can still be an analytic function and we can arrive at alternate IA/IF parameterizations which may be more accurate at describing the latent signal. Although the existence of other IA/IF parameterizations is not new, Vakman argued that the AS is the only physicallyjustifiable complex extension. We argue that by modifying the differential equation for simple harmonic motion shankar2014fundamentals (); FundAcoust (), our parameterizations are also physically justified.
1 Introduction
Many physical phenomenon are characterized by a complex signal
z(t)=x(t)+jy(t)=\rho(t)e^{j\Theta(t)}  (1) 
where \rho(t) is the signal’s IA, \Theta(t) the signal’s phase, and \Omega(t)=\dfrac{d}{dt}\Theta(t) the signal’s IF. We will call z(t) the latent signal because only the real part, x(t) is observed and the imaginary part, y(t) is hidden, i.e. the act of observation corresponds to the real operator
x(t)=\Re\{z(t)\}.  (2) 
It is often desirable to analyze the latent signal—which is known to completely parameterize the physical phenomenon. Thus the complex extension problem becomes that of determining z(t) from the observation x(t), i.e. recovering the quadrature y(t). In the classical approach, we seek a unique rule \mathcal{L}\{\cdot\} such that the estimate of y(t), \hat{y}(t)=\mathcal{L}\{x(t)\}. These relations and the complex extension problem are illustrated in Fig. 1, where we see many latent signals mapping to a single observation under the real operator. In the LSA problem, given x(t) we seek z(t).
In the context of timefrequency signal analysis, we desire the instantaneous parameters, \rho(t) and \Omega(t). Once we have determined a rule for \hat{y}(t), the instantaneous estimates are given by
\hat{\rho}(t)=\pm\left\hat{z}(t)\right=\pm\mathchoice{{\hbox{$\displaystyle% \sqrt{x^{2}(t)+\hat{y}^{2}(t)\,}$}\lower 0.4pt\hbox{\vrule height 8.959863pt d% epth 7.167891pt width 1px}}}{{\hbox{$\textstyle\sqrt{x^{2}(t)+\hat{y}^{2}(t)% \,}$}\lower 0.4pt\hbox{\vrule height 8.959863pt depth 7.167891pt width 1px}}}% {{\hbox{$\scriptstyle\sqrt{x^{2}(t)+\hat{y}^{2}(t)\,}$}\lower 0.4pt\hbox{% \vrule height 8.399872pt depth 6.719897pt width 1px}}}{{\hbox{$% \scriptscriptstyle\sqrt{x^{2}(t)+\hat{y}^{2}(t)\,}$}\lower 0.4pt\hbox{\vrule h% eight 8.399872pt depth 6.719897pt width 1px}}}  (3) 
and
\hat{\Omega}(t)=\frac{d}{dt}\left[\arctan\left(\frac{\hat{y}(t)}{x(t)}\right)% \right].  (4) 
The very definition of IA and IF hinges on \hat{y}(t) and hence \mathcal{L}\{\cdot\}.
The complex extension problem is wellknown. In 1937, Carson and Fry formally defined the IF based on the phase derivative of a complex FM signal CarsonFry (); VanDerPol1946 (). In Gabor’s seminal 1946 paper gaborTOC (), he introduced a Quadrature Method (QM) as a practical approach for obtaining the complex extension of a real signal and showed its equivalence to the HT. In this context, the rule is given by
\hat{y}(t)=\mathcal{H}\{x(t)\}  (5) 
where \mathcal{H}\{\cdot\} is the HT operator and
\hat{z}(t)=x(t)+j\mathcal{H}\{x(t)\}  (6) 
is termed the AS. This is illustrated in Fig. 1 where we see using the HT to estimate the quadrature leads us to a subset of the latent signals.
However shortly afterward, Shekel pointed out the ambiguity problem in Ville’s work defining the instantaneous parameters of a real signal shekel1953instantaneous (); picinbono1997instantaneous (); BoashashP1 (). As an example, consider
x(t)=\Re\left\{a_{0}(t)e^{j[\int\limits_{\infty}^{t}\omega_{0}(\uptau)d\uptau% +\phi_{0}]}\right\}.  (7) 
There are an infinite set of pairs a_{0}(t) and \omega_{0}(t) for which x(t) may be equivalently described and hence an infinite set of IA/IF parameterizations. Shekel pointed this ambiguity out “with the hope of banishing it [IF] forever from the dictionary of the communication engineer.” Others such as Hupert, suggested that despite the ambiguity, the concept of instantaneous parametrization of real signals was still useful and could be applicable Hupert1953instantaneous (); rihaczek1966hilbert ().
In order to constrain the ambiguity problem so that a unique complex extension can be justified for a real signal, Vakman proposed three conditions tied to physical reality vakman1972 (); vakman1976 (); vakman1979 (): 1) amplitude continuity, 2) phase independence of scaling and homogeneity, and 3) HC. With these conditions, Vakman showed that the unique complex extension is given by the rule \hat{y}(t)=\mathcal{H}\{x(t)\}.
More recently, other authors have proposed other conditions to constrain the ambiguity, such as bounded amplitude and bounded IF variation, leading to alternate IA/IF solutions loughlin1996amplitude (); loughlin1998bounded (). Finally, we point out that a number of other methods have also been introduced, without explicitly stating the conditions. Vakman has shown that all these methods violate one or more of the conditions that he proposed and almost all retain HC vakman1996analytic (); VakmanBook (). Like many authors, we believe Vakman’s conditions justifying Gabor’s use of the HT are reasonable and sensible for several special cases of z(t). On the contrary, we believe that Vakman’s conditions, in particular HC, can be too restrictive for more general z(t) and in the worst case, may lead to incorrect results and interpretations.
Part I is organized as follows. In Section 2, we review the HT and its motivations and the AS. In Section 3, we give a proof showing we can still maintain analyticity without HC and also give a proof that no universal rule for the quadrature exists. In Section 4, we provide an example of a LSA problem and several solutions. Finally, in Section 5 we summarize Part I.
2 The Hilbert Transform and Analytic Signal
Although there exist several methods for estimating the instantaneous parameters, use of the HT dominates science and engineering. The HT of x(t) is given by
\mathcal{H}\{x(t)\}\equiv\frac{1}{\pi}\mathchoice{{\vbox{\hbox{$\textstyle$}% }\kern9.374857pt}}{{\vbox{\hbox{$\scriptstyle$}}\kern7.968628pt}}{{\vbox{% \hbox{$\scriptscriptstyle$}}\kern5.624914pt}}{{\vbox{\hbox{$% \scriptscriptstyle$}}\kern4.687428pt}}\!\int\limits_{\infty}^{\infty}\frac{% x(\uptau)}{\uptaut}d\uptau  (8) 
where \mathchoice{{\vbox{\hbox{$\textstyle$}}\kern7.499886pt}}{{\vbox{\hbox{$% \scriptstyle$}}\kern6.374903pt}}{{\vbox{\hbox{$\scriptscriptstyle$}}\kern4% .499931pt}}{{\vbox{\hbox{$\scriptscriptstyle$}}\kern3.749943pt}}\!\int indicates the Cauchy principle value integral titchmarsh (); LeeWienerLegacy (). The three main motivations for use of the HT are: 1) Vakman’s physical conditions, 2) analyticity of the resulting complex extension, and 3) ease of computation via Gabor’s QM. In this section, we review the motivations.
2.1 Vakman’s Physical Conditions
Vakman proposed conditions in order to constrain the ambiguity in choosing the complex extension vakman1972 (); vakman1976 (); vakman1979 (); vakmanVainshtein1977 (); VakmanBook ().
Condition 1 Amplitude Continuity: Simply stated, amplitude continuity requires that the IA, \rho(t) is a continuous function. This implies that the rule, \hat{y}(t)=\mathcal{L}\{x(t)\} must be continuous, i.e.
\mathcal{L}\{x(t)+\epsilon w(t)\}\rightarrow\mathcal{L}\{x(t)\}~{}\mathrm{for}% ~{}\epsilon w(t)\rightarrow 0.  (9) 
Condition 2 Phase Independence of Scaling and Homogeneity: Let x(t) have a complex extension, z(t)=\rho(t)\exp[j\Theta(t)]. Then for a real constant c>0, cx(t) has associated complex extension z_{1}(t)=[c\rho(t)]\exp[j\Theta(t)], i.e. only the IA of the complex representation is affected and \Theta(t) and \Omega(t) remain unchanged. This implies that the rule for performing the complex extension is scalable
\mathcal{L}\{cx(t)\}=c\mathcal{L}\{x(t)\}.  (10) 
Conditions 1 and 2 force the operator to be linear VakmanBook ().
Condition 3 Harmonic Correspondence: Let x(t)=a_{0}\cos(\omega_{0}t+\phi_{0}), then HC forces the complex extension,
\hat{z}(t)=a_{0}e^{j(\omega_{0}t+\phi_{0})},  (11) 
where we note the IA and IF must be constant. This implies that
\mathcal{L}\{a_{0}\cos(\omega_{0}t+\phi_{0})\}=a_{0}\sin(\omega_{0}t+\phi_{0})  (12) 
and \hat{z}(t) is a Simple Harmonic Component (SHC) with positive IF.
As Vakman showed vakman1972 (), the HT is the only operator that satisfies the above conditions and as a result, the HT (or Gabor’s practical QM implementation) is viewed as the correct way to complex extend a signal. With his work, Vakman was able to refute most objections as to whether use of the HT is physically justified.
Condition 4 Phase Continuity: In addition to real signals having an ambiguity in instantaneous parameterization, complex signals also have an ambiguity: although we construct \hat{z}(t) using a rule for y(t), we want the IA/IF pair \rho(t),~{}\Omega(t) for signal analysis, and there can be ambiguity in this coordinate transformation cohenLoughlinVakman1999 (). There exist at least two choices for resolving this ambiguity when it arises:

Condition 4a Positive IA: \rho(t)\geq 0~{}\forall~{}t.

Condition 4b Phase continuity: the phase \Theta(t) is a continuous function.
Although condition 4a is the traditional choice, we advocate condition 4b to ensure the IF is well defined. This is apparent in (3) where the IA may be negative.
2.2 Gabor’s Quadrature Method
The HT is an ideal operator that in practice is not physically realizable. Gabor’s QM provides a frequency domain approach that under certain conditions is equivalent to the HT. The steps in Gabor’s QM are:

decompose x(t) into SHCs, i.e. compute the Fourier spectrum

double the magnitude of the nonnegative frequency components

negate the negative frequency components.
Gabor’s QM results in a complex signal formulated in terms of nonnegative spectral frequencies,
\hat{z}(t)=\sum\limits_{k=0}^{K1}a_{k}\exp\left\{j\left[\omega_{k}t+\phi_{k}% \right]\right\}  (13) 
i.e. we have decomposed the signal into SHCs each having constant IA, a_{k} and (nonnegative) constant IF, \omega_{k}. By comparing the expressions for \hat{z}(t) in (13) and (1), we can effectively collapse K SHCs into a single AMFM component and then obtain an IA/IF pair. Implementing Gabor’s QM using a FT is very convenient and is one reason for the method’s popularity.
2.3 Analyticity of the Analytic Signal
Unfortunately, the word “analytic” has two distinct meanings when used in signal processing and in addition, another meaning in mathematics. A signal is said to be analytic if it consists only of nonnegative frequency components cohen1995time (); boashash2003time (). The analytic signal refers to the complex extension of a real signal using the HT, i.e. Gabor’s method ville (); bedrosian1962analytic (); cohen1995time (); xia1999analytic (); boashash2003time (). If z({\text{}}), where {\text{}}\equiv t+j\tau, is an analytic function then the real and imaginary parts of the complex function,
z({\text{}})=u(t,\tau)+jv(t,\tau)  (14) 
satisfy the CauchyRiemann (CR) conditions
\dfrac{\partial}{\partial t}u\left(t,\tau\right)=\dfrac{\partial}{\partial\tau% }v\left(t,\tau\right)~{}~{}~{}\text{and}~{}~{}~{}\dfrac{\partial}{\partial\tau% }u\left(t,\tau\right)=\dfrac{\partial}{\partial t}v\left(t,\tau\right).  (15) 
For the AS \hat{z}(t)=x(t)+j\mathcal{H}\{x(t)\}, such that t\leftarrow{\text{}}\equiv t+j\tau, the complex function \hat{z}({\text{}})=\hat{u}(t,\tau)+j\hat{v}(t,\tau) is an analytic function brown2009complex (). Hence the reason for calling a HTextended signal the AS ville (); cohen1995time ().
3 Relaxing the Condition of Harmonic Correspondence
3.1 On Harmonic Correspondence
We can understand Vakman’s motivation for tying HC to physical reality by considering the differential equation
\dfrac{d^{2}}{dt^{2}}z(t)+\omega_{0}^{2}z(t)=0  (16) 
which describes many ideal systems, e.g. a LC circuit or mass/spring model. The solution to (16) is the SHC in (11).
Any deviation from this ideal case, requires modification of the differential equation. For example, when the differential equation describes a circuit and resistance is included or describes motion and damping is included, the equation becomes
\dfrac{d^{2}}{dt^{2}}z(t)+c\dfrac{d}{dt}z(t)+\omega_{0}^{2}z(t)=0  (17) 
where c is a constant. In this case, the solution includes an AM term and has the form
z(t)=a_{0}e^{\nu t}e^{j(\omega_{d}t+\phi_{0})}  (18) 
which is not a SHC FundAcoust (). As another example, when the differential equation is further modified to include timevarying coefficients, nonlinearities, or partial derivatives with respect to \tau or spacial variables, the resulting solution may include both AM and FM terms, which is also not a SHC van1934nonlinear (); VakmanBook (); farlow2012partial ().
While most authors believe HC is a reasonable condition to describe physical phenomena, we advocate that this condition is overly constraining. For many analysis problems, this can lead to incorrect interpretations because real physical systems always deviate from the ideal case. By not assuming HC, we gain a degree of freedom in our analysis that allows us to construct other complex extensions that may be better suited to describing the underlying physical phenomena associated with the signal. We believe that any attempt to find a unique rule to infer z(t) of the form in (1) from x(t) in the form of (2) is fundamentally flawed and that no such universal rule can exist.
3.2 On Analyticity of the Complex Extension
We believe that the use of the term “analytic” has in most cases become too restrictive in signal processing. To wit, one may falsely believe that the HT is the only way to complexextend a real signal in order to result in an analytic function, when time is considered complex. Other complex extensions of real signals can be constructed that result in analytic functions. Choosing \hat{y}(t)=\mathcal{H}\{x(t)\} to obtain the AS \hat{z}(t) ensures z({\text{}}) is an analytic function. However, there are other choices for \hat{y}(t)\neq\mathcal{H}\{x(t)\} such that z({\text{}}) is an analytic function.
Theorem If we do not assume HC, there exists at least one choice for the quadrature, y(t)\neq\mathcal{H}\{x(t)\} that results in z({\text{}}) being an analytic function.
Proof.
Let x(t)=a_{0}\cos(\omega_{0}t) and choose the quadrature such that
z(t)=a_{0}\cos(\omega_{0}t)+j\alpha a_{0}\sin(\beta\omega_{0}t)  (19) 
with real \alpha and \beta. The complex function is given by
z({\text{}})=a_{0}\cos(\omega_{0}[t+j\tau])+j\alpha a_{0}\sin(\beta\omega_{0}[% t+j\tau])  (20) 
where
u(t,\tau)=a_{0}\cos(\omega_{0}t)\cosh(\omega_{0}\tau)\alpha a_{0}\cos(\beta% \omega_{0}t)\sinh(\beta\omega_{0}\tau)  (21) 
and
v(t,\tau)=\alpha a_{0}\sin(\beta\omega_{0}t)\cosh(\beta\omega_{0}\tau)a_{0}% \sin(\omega_{0}t)\sinh(\omega_{0}\tau).  (22) 
It can easily be shown that this choice leads to z({\text{}}) satisfying the CR conditions and hence is an analytic function. Any choice of \alpha\neq 1 or \beta\neq 1 does not imply HC. ∎
Although Gabor’s QM provides a simple rule to obtain an IA/IF pair that parameterizes x(t), it may not address the problem of obtaining the latent signal from the observation. The heart of the problem is that no single rule can determine the latent signal from the observation for every possible latent signal because more than one complex signal map to the same real signal under the real operator. Although a unique rule can be constructed to work for a particular z(t), the same rule cannot generally be used.
Corollary No unique rule for the quadrature, \hat{y}(t)=\mathcal{L}\{x(t)\} exists to obtain the latent signal, z(t) from the observation, x(t)=\Re\{z(t)\} for all z(t).
Proof.
Consider the latent signal, z(t) of the form in (19). Assume a unique rule, \hat{y}(t)=\mathcal{L}\{x(t)\}\equiv\alpha_{0}a_{0}\sin(\beta_{0}\omega_{0}t) exists to obtain z(t). This implies a unique \hat{z}(t)=a_{0}\cos(\omega_{0}t)+j\alpha_{0}a_{0}\sin(\beta_{0}\omega_{0}t). If \alpha_{0}\neq\alpha or \beta_{0}\neq\beta then \hat{z}(t)\neq z(t) and \mathcal{L}\{\cdot\} is not unique. ∎
3.3 On Harmonic Conjugate Functions
If z({\text{}})=u(t,\tau)+jv(t,\tau) is an analytic function, then u(t,\tau) and v(t,\tau) are unique and are harmonic conjugates brown2009complex (). Even though x(t)=u(t,0) and y(t)=v(t,0) in (1) this does not imply that u(t,\tau)=x({\text{}}) and v(t,\tau)=y({\text{}}), but rather u(t,\tau) and v(t,\tau) each contain terms from both x({\text{}}) and y({\text{}}):
\displaystyle z({\text{}})  \displaystyle=  \displaystyle x({\text{}})+jy({\text{}})  (23)  
\displaystyle=  \displaystyle\left[x_{R}(t,\tau)+jx_{I}(t,\tau)\right]+j\left[y_{R}(t,\tau)+jy% _{I}(t,\tau)\right]  
\displaystyle=  \displaystyle\left[x_{R}(t,\tau)y_{I}(t,\tau)\right]+j\left[x_{I}(t,\tau)+y_{% R}(t,\tau)\right]  
\displaystyle=  \displaystyle u(t,\tau)+jv(t,\tau) 
where R,~{}I denote the real and imaginary parts, respectively. This is easily seen in (21), where \alpha and \beta are present in u(t,\tau) despite \alpha and \beta appearing only in y(t) in (19). Thus the problem of finding y(t) cannot be solved by finding a harmonic conjugate of u(t,\tau) because y(t) must be known to obtain u(t,\tau).
3.4 On AM–FM Demodulation
The HT is widely used as a demodulation algorithm for AM–FM signals. This is typically justified as a valid approach because of Vakman and also Bedrosian’s theorem.^{1}^{1}1If the product l(t)h(t) consists of a lowfrequency factor, l(t) and a highfrequency factor, h(t) of nonoverlapping spectra, then \mathcal{H}\{l(t)h(t)\}=l(t)\mathcal{H}\{h(t)\} Bedrosian1963 (); Nuttall1966 (); VakmanBook (). The HT can be used to determine the IA/IF for a small subset of latent signals, i.e. those with HC. The HT can also be used to closely approximate the IA/IF for latent signals whenever \mathcal{H}\{x(t)\}\approx y(t) picinbono1997instantaneous (). However, the HT cannot be used in general to obtain the IA/IF of all latent signals.
4 Example
Recall the LSA problem: given the observation x(t)=\Re\{z(t)\}, find z(t) or more strictly the instantaneous parameters, \rho(t) and \Omega(t). In this section, we will give three solutions to a LSA problem where two of these solutions have \hat{y}(t)\neq\mathcal{H}\{x(t)\}, as illustrated in Fig. 1.
Suppose we have three ideal systems: a frequency modulator, an amplitude modulator, and a Linear TimeInvariant (LTI) system in the steady state. We observe a triangle waveform x(t) that could have come from any of the three ideal systems. What is the corresponding latent signal z(t) or more strictly the instantaneous parameters, \rho(t) and \Omega(t)?
Let the observation be the periodic, evensymmetric triangle waveform where one period is given by
x(t)=\left\{\begin{array}[]{ll}2A\omega_{0}t/\pi+A,&0\leq t\leq T/2\\ 2A\omega_{0}t/\pi+A,&T/2\leq t\leq 0\end{array}\right.  (24) 
with amplitude A, period T, and fundamental frequency \omega_{0}, as illustrated in Fig. 2.
4.1 Solution Assuming Harmonic Correspondence
If we assume HC, i.e. the ideal system is LTI in the steady state then
\displaystyle\hat{z}_{1}(t)  \displaystyle=x(t)+j\mathcal{H}\{x(t)\}  (25a)  
\displaystyle=\sum\limits_{k=0}^{\infty}\frac{8A}{\pi^{2}(2k+1)^{2}}\cos\left[% (2k+1)\omega_{0}t\right]+j\sum\limits_{k=0}^{\infty}\frac{8A}{\pi^{2}(2k+1)^{2% }}\sin\left[(2k+1)\omega_{0}t\right]  (25b) 
where the IA is
\hat{\rho}(t)=\frac{8A}{\pi^{2}}\tilde{a}_{0}(t)  (26) 
and the IF is
\hat{\Omega}(t)=\omega_{0}+\frac{d}{dt}M_{0}(t)  (27) 
where \tilde{a}_{0}(t) and M_{0}(t) are related through
\tilde{a}_{0}(t)e^{jM_{0}(t)}=\sum\limits_{k=0}^{\infty}\frac{1}{(2k+1)^{2}}e^% {j2k\omega_{0}t}.  (28) 
4.2 Solution Assuming Constant IF
If we assume constant IF \hat{\Omega}(t)=\omega_{0}, i.e. the ideal system is an amplitude modulator then
\hat{z}_{2}(t)=\hat{\rho}(t)\cos(\omega_{0}t)+j\hat{\rho}(t)\sin(\omega_{0}t)  (29) 
with the IA given by
\hat{\rho}(t)=x(t)/\cos(\omega_{0}t).  (30) 
In this solution, \hat{y}(t)=\hat{\rho}(t)\sin(\omega_{0}t)\neq\mathcal{H}\{\hat{\rho}(t)\cos(% \omega_{0}t)\}, because Bedrosian’s theorem cannot be applied.
4.3 Solution Assuming Constant IA
If we assume constant IA \hat{\rho}(t)=A, i.e. the ideal system is a frequency modulator then
\hat{z}_{3}(t)=A\cos\left[\omega_{0}t+M_{0}(t)\right]+jA\sin\left[\omega_{0}t+% M_{0}(t)\right]  (31) 
with IF given by
\hat{\Omega}(t)=\omega_{0}+\frac{d}{dt}M_{0}(t)  (32) 
where
\displaystyle M_{0}(t)  \displaystyle=  \displaystyle\arccos\left[x(t)/A\right]\omega_{0}t.  (33) 
As in the solution for constant IF, \hat{y}(t)=A\sin\left[\omega_{0}t+M_{0}(t)\right]\neq\mathcal{H}\{A\cos\left[% \omega_{0}t+M_{0}(t)\right]\}.
5 Summary
We have presented the LSA problem and reviewed the motivation for the use of the HT and AS in signal analysis. We proved the HC condition is not necessary for analyticity and that no unique rule for the complex extension exists to obtain the latent signal. It was argued that by relaxing the HC condition, which forces the use of the HT, we gain alternate choices for the latent signal. Although the HT is widely used for demodulation of AM–FM signals, it cannot be used to obtain the IA/IF of all latent signals. In a strict sense, Gabor’s QM can only be used when the latent signal is a superposition of SHCs with nonnegative IF and can approximate the latent signal whenever \mathcal{H}\{x(t)\}\approx y(t), e.g. communication signals where Bedrosian’s theorem is approximately satisfied.
Part II. Hilbert Spectral Analysis Using Latent AM–FM Components
\pdfstringdef\BKM@titlePart II. Hilbert Spectral Analysis Using Latent AM–FM Components
[W]ellknown methods … can no longer be applied when the frequency itself is made a function of the time. Even the concept of “instantaneousfrequency”… is of a somewhat arbitrary but nevertheless highly useful nature, is often misunderstood and even misinterpreted …. The only way at present available to solve these and similar problems is to go back to the very first and fundamental principles. This implies a theoretical treatment beginning with the differential equations of the problem concerned. Unfortunately these equations can seldom be solved in terms of wellknown functions, such as real or complex exponentials [simple harmonic components]. I think it is precisely to this fact that the main difficulties can be traced.—Van der Pol VanDerPol1946 ()
In the second part of the paper, the central theme is Hilbert Spectral Analysis (HSA) using a superposition of latent AM–FM components. We present HSA as a generalized LSA problem. In the general problem, we seek a representation of z(t) consisting of a superposition of latent components, i.e. a multicomponent model. Furthermore, rather than seeking a single IA/IF pair for z(t), we seek a set of IA/IF pairs each associated with the components. Although timefrequency analysis has been extensively studied the use of a generalized AM–FM model for this analysis, without the HC condition, has never been proposed. Using this model leads to nonunique signal decompositions due to the theorem proved in Part I. However, by imposing assumptions on the form of the AM–FM component, a unique parameterization in terms of IA and IF can be obtained. This analysis requires abandoning Gabor’s complex extension as advocated in Part I and instead allowing the assumptions to imply the complex extension. This model enables us to analyze signals with very few restrictions resulting in alternate and possibly more useful decompositions, especially for nonstationary signals. In this part, AM–FM modeling and HSA theory is presented. Examples using HSA are given and a visualization of the Hilbert spectrum is proposed.
1 Introduction
It has long been known that simple harmonic analysis of nonstationary signals may lead to incorrect interpretations of an underlying signal model BoashashP1 (). In Gabor’s seminal paper, “Theory of Communications” he writes:
The greatest part of the theory of communication has been built up on the basis of Fourier’s reciprocal integral relations …. Though mathematically this theorem is beyond reproach … even experts could not at times conceal an uneasy feeling when it came to the physical interpretation of results obtained by the Fourier method. After having for the first time obtained the spectrum of a frequencymodulated sine wave, Carson wrote: ‘The foregoing solutions, though unquestionably mathematically correct, are somewhat difficult to reconcile with our physical intuitions ….’ gaborTOC ()
As Gabor noted, Carson in 1922 was the first to rigorously study a FM signal and realize that the frequency components obtained from such a nonstationary signal, do not describe the physical system in a meaningful way Carson1922Notes (). A similar argument regarding an AM signal was made by Priestly, who wrote, “…a nonstationary process in general cannot be represented in a meaningful way by the simple Fourier expansion” priestley1988non (). His example,
x(t)=Ae^{t^{2}/o^{2}}\cos(\omega_{0}t+\phi_{0})  (34) 
and its Fourier Transform (FT) consists of two Gaussian functions centered at frequencies \pm\omega_{0}. Thus, the FT contains an infinite number of SHCs. This interpretation of the underlying signal model may be incorrect. For example, an alternative signal model consists of just two components at constant frequencies \pm\omega_{0}, with each component having a timevarying amplitude, (A/2)e^{t^{2}/o^{2}}. Mathematically, these two representations are equally valid and correspond to different families of basic functions used for representation BoashashP1 ().
For nonstationary signals, SHCs lose effectiveness and thus the idea of timevarying components^{2}^{2}2Some authors, such as Ville ville (), refer to timevarying components as “instantaneous spectra” or more generally as functions of time that give the structure of a signal at a given instant., i.e. components with timevarying IA and IF has arisen in order to account for the nonstationarity priestley1988non (); BoashashP1 (). However, the concept of IF is not without its own controversy primarily due to four reasons:

There is an apparent paradox in associating the words “instantaneous” and “frequency” because frequency usually defines the number of cycles undergone during one unit of time BoashashP1 ().

Without assumptions, instantaneous parameterizations of a signal are not unique shekel1953instantaneous (); BoashashP1 (); cohen1995time (); loughlin1998bounded ().

The commonly accepted definition of IF, i.e. phase derivative of Gabor’s AS is correct for only a limited class of signals cohen1995time ().

Although different quantities, harmonic frequency and IF are often confused likely due to the term “frequency” attached to both fink1966relations (); mandel1974interpretation (). Harmonic frequency is a special case of IF and is only equivalent under the assumption of SHCs.
Several authors over the previous decades have shown that problems and paradoxes exist that are related to the definition of IF shekel1953instantaneous (); fink1966relations (); loughlin1998bounded (); cohen1995time (); mandel1974interpretation (); gupta1975definition (); vakmanVainshtein1977 (); priestley1988non (); loughlin1997comments (); oliveira1998concept (); oliveira1999instantaneous (); nho1999instantaneous (); chui2015signal (). However on the whole, these problems seem to have been forgotten or ignored. In 1937, Carson and Fry formally defined the IF based on the phase derivative of a complex FM signal assuming the FM message is slowlyvarying CarsonFry (). This assumption, while perfectly valid in communication theory, implies a narrowband component when used in signal analysis.
In order to exploit the mathematical convenience of using complex exponentials in signal analysis, Gabor introduced a QM for obtaining a complex extension of a real signal gaborTOC (). Gabor’s QM assumed positive IF SHCs and was shown to be equivalent to the HT gaborTOC (). Although Gabor’s QM provides a unique method for complex extension, as Cohen pointed out without strict adherence to the assumption, this results in many counterintuitive consequences cohen1995time (). Ville defined the IF of a real signal by using Gabor’s complex extension and then Carson’s definition of IF BoashashP1 (). By defining the IF as the derivative of the phase of Gabor’s AS, Ville was able to show the average harmonic frequency is equal to the time average of the IF nho1999instantaneous (). He then formulated the WignerVille Distribution (WVD) and showed that the first moment of the WVD with respect to frequency yields the IF. Using Gabor’s AS to extend a real signal results in a number of very convenient relations. As a result, Gabor’s QM is almost universally viewed as the correct way to define the complex signal and subsequently the correct way to define IA, IF, and phase for real signals despite the inherent assumption of HC cohen1995time ().
In our research, we have found that the dogmatic use of Gabor’s complex extension does not provide the necessary flexibility for modeling nonstationary signals. For many problems, as we will demonstrate, it can be advantageous to allow the assumptions of an underlying signal model to imply the complex extension. As a result, assumptions must be made on a per problem basis in order to force a unique decomposition and hence, in order to properly estimate the latent signal and its components. This leads to a decomposition into complex AM–FM components rather than SHCs. As a result of using AM–FM components, we revert back to Carson’s original definition of IF and by further relaxing HC and embracing the resulting nonuniqueness, the associated problems and paradoxes related to IF can all be resolved.
Part II of this paper is organized as follows. In Section 2, we propose the AM–FM model. In Section 3, we propose a general Hilbert spectrum based on the AM–FM model and examine several familiar specializations. In Section 4, a frequency domain view of LSA is presented. In Section 5, we discuss subtleties of complex extension of real signals assuming AMFM components with relaxed HC. In Section 6, we provide examples of HSA and in Section 7 propose a visualization of the Hilbert spectrum to aid in the interpretation. In Section 8, we summarize Part II.
2 The AM–FM Model
2.1 Definitions and Assumptions
Our goal in HSA is to decompose a signal into complex AM–FM components, each of which yield IA and IF estimates that are matched to some criteria. That criteria could be some physical or perceptual observation or simply some intuitive notion of an underlying signal model. The decomposition problem is illustrated in Fig. 3 where we have extended the LSA problem shown in Fig. 1 to show models of the latent signal that in general, are composed of a set of latent AMFM components.
We propose to define the AM–FM model for z(t) as a superposition of K (possibly infinite) AM–FM components
z(t)\equiv\sum\limits_{k=0}^{K1}\psi_{k}\left(t;a_{k}(t),{\omega}_{k}(t),\phi% _{k}\right)  (35) 
where the AM–FM component is defined as
\displaystyle\psi_{k}(t;a_{k}(t),{\omega}_{k}(t),\phi_{k})  \displaystyle\equiv a_{k}(t)\exp\left\{j\left[\int\limits_{\infty}^{t}\omega_% {k}(\uptau)d\uptau+\phi_{k}\right]\right\}  (36a)  
\displaystyle=a_{k}(t)e^{j\theta_{k}(t)}  (36b)  
\displaystyle=s_{k}(t)+j\sigma_{k}(t)  (36c) 
parameterized by the IA a_{k}(t), IF \omega_{k}(t), and phase reference \phi_{k}. We assume the observed real signal x(t) is related to the latent signal z(t) by (2).
When convenient, we will drop the k denoting the parameters of the kth component for notational simplification and unspecified references to IA, IF, and phase refer to a_{k}(t), \omega_{k}(t), and \theta_{k}(t) rather than \rho(t), \Omega(t), and \Theta(t), respectively. We note that most of the paradoxes related to IF are a result of not correctly interpreting between these variables in the multicomponent case, because in the monocomponent case they are equivalent. The geometric interpretation of the AM–FM component in (36) is illustrated with the Argand diagram in Fig. 4. The AM–FM component can be visually interpreted as a single rotating vector in the complex plane with timevarying length and timevarying angular velocity.
Carson expressed the phase \theta(t) in terms of a carrier frequency \omega_{c} and FM message m(t) as CarsonFry ()
\theta(t)=\omega_{c}t+\lambda\int\limits_{0}^{t}m(\uptau)d\uptau  (37) 
assuming the modulation index, \lambda\leq\omega_{c} and m(t)\leq 1. We choose to parameterize the phase as
\displaystyle\theta(t)=\omega_{c}t+\int\limits_{\infty}^{t}m(\uptau)d\uptau+\phi  (38) 
which is more general than (37), avoids an arbitrary normalization on m(t), and does not impose an upper bound on the IF. The phase can also be expressed in terms of \omega_{c} and Phase Modulation (PM) message M(t) as
\displaystyle\theta(t)  \displaystyle=  \displaystyle\omega_{c}t+M(t)+\phi  (39) 
where the relationship between the messages is given by
M(t)=\int\limits_{\infty}^{t}{m}(\uptau)d\uptau.  (40) 
In the AM–FM model, we assume nonnegative IF
{\omega}(t)=\frac{d}{dt}\theta(t)=\omega_{c}+m(t)\geq 0\ \forall\ t.  (41) 
In Section 5, we discuss situations in which this assumption may be relaxed. Finally, we assume without loss of generality, that the phase and carrier references are taken at t=0, i.e.
\int\limits_{\infty}^{0}{m}(t)dt=M(0)=0~{}~{}~{}\implies~{}~{}~{}\phi=\theta(% 0)~{}~{}\text{and}~{}~{}\omega_{c}=\omega(0)m(0).  (42) 
2.2 Monocomponents and Narrowband Components
The concept of a signal being composed of one (mono) or more (multi) components and whether those components are narrowband or wideband has caused much confusion in existing literature. There is no clear agreed upon definition for a monocomponent huang1998empirical (), however, a few definitions have been proposed boashash1992time (); cohen1995time (); boashash2003time (). Some authors describe a monocomponent as defined by a single ‘ridge’ in time and harmonic frequency, corresponding to an elongated region of energy concentration boashash1992time (); cohen1992multicomponent (); cohen1995time (). Cohen agrees that if a component is well localized, the crest of the ridge corresponds to the IF and further states that the width of the ridge depends on the energy spread, or instantaneous bandwidth of the component cohen1992multicomponent (); cohen1995time (). A multicomponent signal may then be defined as any signal which is the sum of two or more monocomponents which can only occur if the separation between the ridges is large in comparison to bandwidth of the components boashash1992time (); cohen1992multicomponent (); cohen1995time (); wei1998instantaneous (). It has also been noted that a signal may be monocomponent at some time instances and multicomponent at other time instances cohen1992multicomponent ().
With these descriptions of a monocomponent, we point out that the definitions and concepts of IA and IF of a complex AM–FM component as defined by (36) are by no means justified only for narrowband components. Also, as pointed out by Cohen cohen1995time (), a multicomponent signal is not defined by the ability to express a signal as a sum of parts, because there are an infinite number of ways to write a signal as a sum of parts. Rather, it is the nature of the parts in relation to themselves and the signal which determines whether the decomposition is of interest cohen1988instantaneous (). We strongly believe that the previous description of a monocomponent is unnecessarily restrictive, limiting the monocomponent to be narrowband. Clearly, the IA and IF in (36) are welldefined for wideband components.
Simply put, a perfectlyvalid component can be defined with (36) while foregoing the restrictive narrowband definitions of Carson CarsonFry (), Ville ville (), Boashash BoashashP1 (); boashash1992time (), and Cohen cohen1995time (). Such wideband components may result in multiple ridges which may be broken down into narrowband components. However, for many problems a wideband component, such as the AM–FM component in (36), can be much more appropriate than multiple, narrowband components. Further, the appearance of structure in the Fourier spectrum can be viewed as an indication that a wideband component is present in the signal. For narrowband components the Fourier and Hilbert spectra can be quite similar, but for wideband components these can be quite different.
3 Hilbert Spectral Analysis
Huang’s original definition of the Hilbert spectrum uses Empirical Mode Decomposition (EMD) to determine a set of Intrinsic Mode Functions (IMFs) which are individually demodulated with the HT to obtain \{a_{k}(t)\} and \{\omega_{k}(t)\} huang1998empirical (). This analysis is also known as the HilbertHuang Transform (HHT) huang1998empirical (). This definition can be generalized, by recognizing that IMFs are a class of AM–FM components and that other decomposition and demodulation methods exist for obtaining the instantaneous parameters, \{a_{k}(t)\} and \{\omega_{k}(t)\}. This generalization can lead to a more powerful and useful signal analysis technique as we will demonstrate.
Towards a more general definition, we define the Hilbert spectrum as a representation or characterization of a signal by an instantaneous spectrum which is parameterized by AM–FM components as in (36). With this definition the Hilbert spectrum is not unique, therefore two problems must be addressed in order to obtain a unique solution:

Determine the appropriate complex extension to obtain an AM–FM decomposition with meaningful interpretation.

Estimate \{a_{k}(t)\} and \{\omega_{k}(t)\}, 0\leq k\leq K1.
The instantaneous parameters are defined on a complex signal, therefore to perform analysis on a real signal it must be extended to the complex space. There are cases where the extension is straightforward, i.e. simple harmonic motion, or in communications where assumptions can be matched at the transmitter and receiver, or other cases where the assumptions are implicit, i.e. Fourier analysis. However, in the general case, signal decomposition and component demodulation are ambiguous and due to this, assumptions must always be made in order to arrive at a unique solution. We will address subtleties of complex extension of real signals assuming AMFM components with relaxed HC in Section 5.
The class of the solution obtained is dependent on the nature of the assumptions. Fig. 5 illustrates the various classes of components under particular assumptions including the SHC, AM, FM, and AM–FM components as well as the IMF (discussed further in Part III). The AM–FM component in (36) may be considered a generalization of all of these components. In the next subsections, we illustrate several forms of the AM–FM component and AM–FM model under various assumptions.
3.1 Two Conventions for Relating the Real Observation to the Latent Signal
As noted by Gabor in gaborTOC (), there are two conventional ways to relate a real signal x(t) to a complex signal z(t). The first convention is
\displaystyle x(t)  \displaystyle=\Re\{z(t)\}  (43a)  
\displaystyle=\Re\{x(t)+jy(t)\}  (43b) 
which can be thought of as a single vector as in Fig. 4 and the second convention is
\displaystyle x(t)  \displaystyle=\left[{z}(t)+{z}^{*}(t)\right]/2  (44a)  
\displaystyle=[x(t)+jy(t)+x(t)jy(t)]/2  (44b) 
which can be thought of as two vectors. If z(t) is a superposition of AM–FM components, the first convention can be viewed as a single rotating vector per component, while the second convention can be viewed as a pair of vectors per component rotating in opposite directions.
Equation (44) is implicit in standard Fourier analysis, with the interpretation that the latent signal is realvalued and its spectrum^{3}^{3}3In this work, use of the unspecified word “spectrum” strictly means Fourier spectrum and not “Hilbert spectrum.” consists of Hermitian symmetric, positive and negative frequency components. The convention in (43) is adopted for the AM–FM model in (35), with the interpretation that only the real part of the signal, x(t) is observed. Inherent in both these conventions is the ambiguity associated with the quadrature y(t), i.e. an infinite number of choices exist for y(t) that lead to the same x(t). Thus, we view the quadrature y(t) as a free parameter that can be chosen to achieve alternate interpretations.
We use the term quadrature in a general sense. For mathematical and physical reasons, the quadrature signal is most often chosen offset in phase by onequarter cycle (\pi/2 radians) relative to the real signal. More specifically, if we assume a single SHC, z(t)=\psi_{0}(t;a_{0},\omega_{0},\phi_{0}) then the AM–FM model is given by (11). Although choosing the quadrature signal given the real signal is trivial for the SHC, the same is not true for the component in (36) where the assumptions a_{0}(t)=a_{0}, and \omega_{0}(t)=\omega_{0} are removed, i.e. the concept of a quadrature signal for the AM–FM component is ambiguous—as will be demonstrated in the following sections.
3.2 Simple Harmonic Component
The simplest form of the AM–FM model has K=1, a_{0}(t)=a_{0}, and \omega_{0}(t)=\omega_{0} in which case the AM–FM component in (36) becomes
\psi_{0}(t;a_{0},\omega_{0},\phi_{0})=a_{0}e^{j(\omega_{0}t+\phi_{0})}.  (45) 
We recognize the significance of this form as representing simple harmonic motion as we discussed in Part I as Vakman’s Condition 3.
3.3 Superposition of Simple Harmonic Components
Perhaps the most familiar form of the AM–FM model has K=\infty, a_{k}(t)=a_{k}, and \omega_{k}(t)=k\omega_{0}
z(t)=\sum\limits_{k=0}^{\infty}a_{k}e^{j(k\omega_{0}t+\phi_{k})}  (46) 
and has real observation given by (2) repeated here as
\displaystyle x(t)  \displaystyle=\Re\left\{z(t)\right\}  (47a)  
\displaystyle=\sum\limits_{k=0}^{\infty}a_{k}\cos(k\omega_{0}t+\phi_{k})  (47b)  
\displaystyle=A_{0}+\sum\limits_{k=1}^{\infty}[A_{k}\cos(k\omega_{0}t)+B_{k}% \sin(k\omega_{0}t)].  (47c) 
Equation (47) is recognized as a Fourier Series (FS) with a particular convention used for the complex extension fourier1807memoire (); Fourier1822analytical (). Note that \phi_{k} is required for proper synthesis since without it, A_{k}=a_{k} and B_{k}=0 and the resulting x(t) is always even, which may not be true. The FS can be considered the simplest of AM–FM superposition models corresponding to the assumption of SHCs.
The FT is the limiting form of the FS as the fundamental, \omega_{0}\rightarrow 0 OppenheimWillsky (). When viewing the FT as a special case of the AM–FM model, subtle and important observations can be made:

In the decomposition problem P1, the AM–FM model corresponding to a separable solution with X(\omega) (not timevarying) and e^{j\omega t} (constant frequency), is the FT, i.e. the inverse FT is a superposition (inner product) of X(\omega) and e^{j\omega t}.

The FT yields a solution that can be described as a superposition of nontimevarying components described by constant state parameters, a_{k}(t)=a_{k}, \omega_{k}(t)=k\omega_{0}, and \phi_{k} which corresponds to a constant model state. However, this does not in any practical way constrain the real superposition x(t). Although analysis of a timevarying system can be accomplished using a constant model, like the FT, in most cases this forces K=\infty which may not accurately describe a physical system.

HSA can be considered a generalization of Fourier analysis huang1998empirical (). However, rather than a generalization of Fourier analysis obtained by generalizing a kernel function (Gabor transform, WignerVille distribution, wavelet transform, etc.) Allen1977ShortTime (); cohen1989time (); jones1990instantaneous (); cohen1993anything (); lovell1993relationship (); cohen1995time (); QianJoint (), an AM–FM component can be viewed as allowing for additional degrees of freedom in the bases of Fourier analysis. Although it may be convenient to think of (36) as a basis for HSA because they span the entire Hilbert space, they do not meet the linear independence property of a formal basis strang09 ().
3.4 The AM Component
One way to generalize the SHC in (45) is to relax the constraint of constant amplitude which leads to an AM component
\displaystyle\psi_{0}(t;a_{0}(t),\omega_{0},\phi_{0})  \displaystyle=  \displaystyle a_{0}(t)e^{j(\omega_{0}t+\phi_{0})}.  (48) 
The simple AM component was originally developed in the context of communication theory armstrong1915some (); voelcker1966towardI (); voelcker1966towardII ().
3.5 Superposition of AM Components
Another form of the AM–FM model is a superposition of AM components that have frequencies as integer multiples of a fundamental
z(t)=\sum\limits_{k=0}^{\infty}a_{k}(t)e^{j(k\omega_{0}t+\phi_{k})}  (49) 
with real observation given by (2).
The AM superposition was first conceived by Gabor and lead to the development of the ShortTime Fourier Transform (STFT) gaborTOC ()
\displaystyle X_{t}(\omega)  \displaystyle=\int x(\uptau)h(\uptaut)e^{j\omega\uptau}d\uptau  (50a)  
\displaystyle\updownarrow  
\displaystyle x_{\omega}(t)  \displaystyle=\frac{1}{2\pi}\int X(\text{})H(\omega\text{})e^{j\text{}t}d% \text{}.  (50b) 
where h(t)\leftrightarrow H(\omega) is a window function. There are two possible interpretations of the STFT in context of the AM–FM model:

The window function h(t) is applied to the signal x(t) in order to use the FT, i.e. simple harmonic analysis at each t. This interpretation couples H(\omega) to X(\omega).

The window function h(t) is applied to the e^{j\omega t} and is effectively an AM superposition analysis where the window is an assumption on a_{k}(t) in the AM–FM model. This interpretation couples H(\omega) to e^{j\omega t}.
The second interpretation provides an important and alternate view of the STFT. What is typically viewed as imprecision in the ability to compute timefrequency parameters in the first interpretation can alternatively be viewed in the second interpretation as a precise ability to compute timecomponent parameters using the modified basis, h(t\uptau)e^{j\omega t}.
3.6 FM Component
Another way to generalize the SHC in (45) is to relax the constraint of constant frequency which leads to a simple FM component
\displaystyle\psi_{0}(t;a_{0},\omega_{0}(t),\phi_{0})  \displaystyle=  \displaystyle a_{0}e^{j[\int\limits_{\infty}^{t}\omega_{0}(\uptau)d\uptau+% \phi_{0}]}.  (51) 
The FM component was also developed in the context of communication theory Carson1922Notes (); CarsonFry (). The FM component has been used in signal synthesis where Chowning observed that a single, wideband FM component is perceived by the ear as spectrally rich, i.e. multiple SHCs FMSynthesis (). This FMsynthesis method has been used in commercial audio synthesizers bedaux1974micro (); chowning1977synthesis ().
3.7 Superposition of FM Components
Signal analysis using a superposition of FM components is not usually considered due to the loss of linear independence of the components in the model and the restriction of constant amplitude. For this reason, signal analysis using a superposition of AM–FM components is more useful than a superposition of FM components.
3.8 Other AM–FM Models
Alternatives to the proposed AM–FM model in Section 2 have been considered. However, common to these are restrictive assumptions which limit the utility as we highlight below. Previous AM–FM models for signal analysis/synthesis usually fall into one of three main groups: 1) HT feldman1994non (); rao2000decomposing (); gianfelici2007multicomponent (); feldman2011hilbert (), 2) peak tracking/sinusoidal modeling mcaulay1986speech (); rao1990estimation (); Pantazis2011 (); Boashash2013 (), and 3) Teager energy operator maragos1993energy (); bovik1993fm (); fertig1996instantaneous (); potamianos1999speech (); boudraa2004if (); boudraa2011instantaneous (). However, some models exist that do not fall into any of these groups Quatieriseparation (); fertig1996instantaneous (). A historical summary of AM–FM modeling is presented by Gianfelici gianfelici2007StateOfArt (). A review of algorithms for estimating IF is presented by Boashash boashash1990algorithms (); BoashashP2 ().
The HT permits the unambiguous definition of IA, IF, and phase of any real signal (random or deterministic) vakmanVainshtein1977 (). Thus, it provides a direct means of performing AM–FM analysis, however, it requires that the quadrature signal be defined in a consistent manner in all cases. Implicit in the HT is HC, as a result the more likely this assumption is true, the better the Hilberttransformed signal approximates the true quadrature signal and the more likely the HTbased solution provides an accurate model of the underlying physical synthesis model boashash1992time (). Direct application of the HT does not however address the signal decomposition problem, as a result, methods for using the HT for decomposition have been proposed gianfelici2007multicomponent (); feldman2008compare (); feldman2011hilbert ().
In peak tracking, one accepts the narrowband definition of a component and as a result, each component appears in the timefrequency plane as a single ridge of energy concentration. Thus, a signal can be parameterized by tracking its ridges in location, intensity, and possibly bandwidth. The timefrequency distribution is usually derived from a STFT mcaulay1986speech () but a generalized timefrequency distribution can also be used rao1990estimation (). Regardless of the timefrequency distribution used, the narrowband assumption is inherent. Extensions of the sinusoidal model have been proposed, such as the harmonic plus noise, and adaptive quasiharmonic model Pantazis2011 ().
The Teager Energy Separation Algorithm (ESA) provides a method of estimating the IA and IF of an AM–FM component, assuming HC. The Teager Energy Operator (TEO) is defined as quatieri02 ()
\displaystyle\varPsi\left\{x(t)\right\}  \displaystyle\equiv  \displaystyle\left[\dot{x}(t)\right]^{2}x(t)\ddot{x}(t)  (52) 
where \dot{x}(t) and \ddot{x}(t) are the first and second time derivatives of x(t), respectively. The TEO applied to a single narrowband AM–FM component in (36), results in quatieri02 ()
\displaystyle\varPsi\left\{\psi_{0}(t)\right\}  \displaystyle\approx  \displaystyle\left[a_{0}(t)\omega_{0}(t)\right]^{2}.  (53) 
Assuming the modulating signals a_{0}(t) and m_{0}(t) are bandlimited, it can be shown that
\displaystyle a_{0}(t)  \displaystyle\approx  \displaystyle\dfrac{\varPsi\left\{\psi_{0}(t)\right\}}{\mathchoice{{\hbox{$% \displaystyle\sqrt{\varPsi\left\{\dot{\psi_{0}}(t)\right\}\,}$}\lower 0.4pt% \hbox{\vrule height 6.999893pt depth 5.599915pt}}}{{\hbox{$\textstyle\sqrt{% \varPsi\left\{\dot{\psi_{0}}(t)\right\}\,}$}\lower 0.4pt\hbox{\vrule height 6.% 999893pt depth 5.599915pt}}}{{\hbox{$\scriptstyle\sqrt{\varPsi\left\{\dot{% \psi_{0}}(t)\right\}\,}$}\lower 0.4pt\hbox{\vrule height 6.999893pt depth 5.5% 99915pt}}}{{\hbox{$\scriptscriptstyle\sqrt{\varPsi\left\{\dot{\psi_{0}}(t)% \right\}\,}$}\lower 0.4pt\hbox{\vrule height 6.999893pt depth 5.599915pt}}}}  (54) 
and
\displaystyle\omega_{0}(t)  \displaystyle\approx  \displaystyle\sqrt{\dfrac{\varPsi\left\{\dot{\psi_{0}}(t)\right\}}{\varPsi% \left\{\psi_{0}(t)\right\}}\,}  (55) 
thus providing a method to estimate narrowband monocomponent IA/IF parameters quatieri02 (); potamianos1994comparison (); diop2011joint (); santhanam2000multicomponent ().
The AM–FM model in these groups all rely on a rigid narrowband component. Surprisingly, the use of wideband components is a wellknown means of synthesis moorer1976synthesis (); chowning1977synthesis (); FMSynthesis (); dodge1997computer (). However, analysis using the wideband components in (36) in the general form, has not been considered. It is our belief that the true power of HSA can only be fully recognized with the use of wideband components without the inherent assumption of HC.
4 Frequency Domain View of Latent Signal Analysis
In Part I, we introduced the LSA problem in the time domain. The frequency domain view of the LSA problem is to estimate the latent spectrum Z(\omega) from X(\omega). Because x(t)=\Re\{z(t)\}, this imposes structure on the spectrum of x(t), namely X(\omega)=[Z(\omega)+Z^{*}(\omega)]/2, i.e. Hermitian symmetry. As discussed in Part I, the conventional way to estimate z(t) from x(t) is via the HT as in (6). Equivalently in the frequency domain, the conventional way to estimate Z(\omega) from X(\omega) is with Gabor’s QM. Using Gabor’s QM effectively forces the FT of (5) \hat{Y}(\omega)=j\mathrm{sgn}(\omega)X(\omega), and \hat{Z}(\omega) is always spectrally onesided Boashash2013 (). In the frequency domain, by relaxing HC we gain the freedom to choose \hat{Y}(\omega) appropriately.
As an example of this frequency domain view, consider the observed spectrum
X(\omega)=a_{0}\pi[\delta(\omega+\omega_{0})+\delta(\omega\omega_{0})]  (56) 
where \delta(\cdot) is the Dirac delta function. If we assume HC,
\hat{Y}(\omega)=ja_{0}\pi[\delta(\omega+\omega_{0})+\delta(\omega\omega_{0})]  (57) 
and the corresponding latent spectrum is given by
\begin{split}\displaystyle\hat{Z}(\omega)&\displaystyle=X(\omega)+j\hat{Y}(% \omega)\\ &\displaystyle=2a_{0}\pi\delta(\omega\omega_{0}).\end{split}  (58) 
If we do not assume HC, there are many choices for \hat{Y}(\omega) leading to many possible latent spectra including,
\hat{Z}(\omega)=2a_{0}\pi\delta(\omega+\omega_{0}),  (59) 
\hat{Z}(\omega)=a_{0}\pi[(1\alpha)\delta(\omega+\omega_{0})+(1+\alpha)\delta(% \omega\omega_{0})],  (60) 
and
\hat{Z}(\omega)=a_{0}\pi[\delta(\omega+\omega_{0})+\delta(\omega\omega_{0})+% \alpha\delta(\omega\beta\omega_{0})\alpha\delta(\omega+\beta\omega_{0})]  (61) 
where \alpha,\beta\in\mathbb{R} and (61) corresponds to the FT of (19). Because of the structure imposed on the spectrum X(\omega)=[Z(\omega)+Z^{*}(\omega)]/2, the latent spectra in (58)(61) all yield the same X(\omega).
In Table 1 we summarize the spectral structure of a complex signal depending on whether or not HC is assumed. If z(t) takes on only real values, then the spectrum has Hermitian symmetry (first row). If z(t) takes on complex values and has HC, then the Fourier spectrum is singlesided (second row). If z(t) takes on complex values and does not have HC, the spectrum does not have Hermitian symmetry.
Class of z(t)  FT Structure 

z(t)\in\mathbb{R}  Z(\omega)=Z^{*}(\omega) 
z(t)\in\mathbb{C} w/ HC  Z(\omega)=0,~{}\omega<0 
z(t)\in\mathbb{C} w/o HC  Z(\omega)\neq Z^{*}(\omega) 
In the frequency domain view of LSA, Gabor’s QM can be interpreted as simply a method to distinguish Z(\omega) from Z^{*}(\omega) which works only if the harmonic frequencies of Z(\omega) are all nonnegative. This is a convenient “trick” and a generalization of this concept for the Hilbert spectrum is possible. As an illustration, consider Fig. 6 where the latent spectrum Z(\omega)=Z_{2}(\omega)+Z_{1}(\omega) and the real signal under analysis has a spectrum [Z_{1}^{*}(\omega)+Z_{2}^{*}(\omega)+Z_{2}(\omega)+Z_{1}(\omega)]/2. Applying Gabor’s QM, the complex signal clearly has the spectrum Z(\omega)=Z_{2}(\omega)+Z_{1}(\omega). Next consider Fig. 6 where the latent spectrum Z(\omega)=Z_{2}(\omega)+Z_{1}(\omega) and the real signal under analysis has a spectrum [Z_{1}^{*}(\omega)+Z_{2}(\omega)+Z_{2}^{*}(\omega)+Z_{1}(\omega)]/2. Applying Gabor’s QM, the complex signal has the spectrum Z_{2}^{*}(\omega)+Z_{1}(\omega) which is the incorrect spectral grouping. Finally, consider Fig. 6 where the latent spectrum Z(\omega)=Z_{3}(\omega)+Z_{2}(\omega)+Z_{1}(\omega) and the real signal under analysis has a spectrum [Z_{1}^{*}(\omega)+\cancel{Z_{2}^{*}(\omega)+Z_{3}(\omega)}+\cancel{Z_{3}^{*% }(\omega)+Z_{2}(\omega)}+Z_{1}(\omega)]/2; cancellation is due to symmetry of the latent spectrum. Applying Gabor’s QM, the complex signal has the spectrum Z_{1}(\omega) and not only has Gabor’s QM failed to yield the latent spectrum but terms are missing entirely.
5 Subtleties of the Hilbert Spectrum
In the previous section, we considered the LSA problem where we assumed SHCs but relaxed the HC assumption. By relaxing HC, we also relax the nonnegative IF constraint associated with it. In this section, we discuss the implications of assuming nonnegative IF and relaxing the assumption of SHCs.
First, consider making the assumption of nonnegative IF and a model composed of SHCs. Table 1, second row shows this implies HC on z(t). Next, consider nonnegative IF and a model not composed of SHCs. Contrary to Table 1 rows two and three where SHCs are assumed, there can exist models with nonnegative IF but without HC. We illustrate the existence of such models with the following example. Consider the complex signal,
z(t)=a_{0}\cos(\omega_{0}t)+j\alpha a_{0}\sin(\omega_{0}t)  (62) 
and the spectrum given by (60) which is not onesided. However in terms of a single AM–FM component, we can have
z(t)=a_{0}(t)e^{j\theta_{0}(t)}  (63) 
where
a_{0}(t)=\mathchoice{{\hbox{$\displaystyle\sqrt{a_{0}^{2}\cos^{2}(\omega_{0}t)% +\alpha^{2}a_{0}^{2}\sin^{2}(\omega_{0}t)\,}$}\lower 0.4pt\hbox{\vrule height % 8.959863pt depth 7.167891pt width 1px}}}{{\hbox{$\textstyle\sqrt{a_{0}^{2}% \cos^{2}(\omega_{0}t)+\alpha^{2}a_{0}^{2}\sin^{2}(\omega_{0}t)\,}$}\lower 0.4% pt\hbox{\vrule height 8.959863pt depth 7.167891pt width 1px}}}{{\hbox{$% \scriptstyle\sqrt{a_{0}^{2}\cos^{2}(\omega_{0}t)+\alpha^{2}a_{0}^{2}\sin^{2}(% \omega_{0}t)\,}$}\lower 0.4pt\hbox{\vrule height 6.999893pt depth 5.599915pt % width 1px}}}{{\hbox{$\scriptscriptstyle\sqrt{a_{0}^{2}\cos^{2}(\omega_{0}t)+% \alpha^{2}a_{0}^{2}\sin^{2}(\omega_{0}t)\,}$}\lower 0.4pt\hbox{\vrule height 6% .999893pt depth 5.599915pt width 1px}}}  (64) 
and
\theta_{0}(t)=\arctan\left[\dfrac{\alpha a_{0}\sin(\omega_{0}t)}{a_{0}\cos(% \omega_{0}t)}\right].  (65) 
The IF, \omega_{0}(t)=\dfrac{d}{dt}\theta_{0}(t) is strictly positive for any positive choice of \alpha. Thus the associated Hilbert spectrum is onesided even though the spectrum is twosided. In other words, saying we assume positive IF is not the same as saying we assume a onesided spectrum. Although signals can be designed such that Z(\omega), for practical purposes, has a onesided spectrum, e.g. communications signals, this does not mean that naturallyoccurring signals have, in general, a onesided spectrum. Therefore, making the onesided spectrum assumption (and assuming HC) may lead to incorrect model parameters.
Another incorrect assumption that is often made is that if the signal is narrowband, the HT will yield a meaningful complex extension. Unfortunately, this is not the case because narrowband signals exist that do not have Hermitian symmetry and hence use of the HT will yield incorrect results. This is demonstrated by (60) which could be considered narrowband, yet use of the HT to extend its real observation would not yield the correct latent signal due to the absence of HC.
By changing from SHCs to AM–FM components and relaxing the assumption of HC, we have generalized the spectrum to a Hilbert spectrum and effectively changed the definition of IF to be componentdependent. By reverting back to Carson’s definition and not using SHCs in the analysis, we move away from Gabor and Ville’s specialized definition back towards the generalized definition of IF.
While using AM–FM components to define IF, which can change in time, we have to consider the possibility that a component’s IF changes sign. Such a sign change is not possible under the assumption of SHCs which have constant IF. In this work, we have arbitrarily assumed that all components must have nonnegative IF for all t, although there may be some signal classes for which this is not true. In these classes, the AM–FM parameters may not match the underlying signal model parameters, although the superposition will yield the correct real signal. This is illustrated in Fig. 7 where a component of z(t) with associated IF \omega(t) and the component of z^{*}(t) with associated IF \omega(t) cannot be separated by assuming nonnegative IF because each has both positive and negative instantaneous frequencies. Therefore in these cases, we need to relax the the assumption of nonnegative IF at some time instances in order to properly estimate the latent signal.
6 Examples using the AM–FM Model
The example in Part I illustrated solutions to the LSA problem where we arrived at a single pair of IA/IF parameters, \rho(t) and \Omega(t), for the latent signal. In this section, we illustrate HSA consisting of a set of IA/IF pairs, \{a_{k}(t),\omega_{k}(t)\} for the signal, one pair for each latent component. This is illustrated by the set diagram in Fig. 3. In our closedform solutions, we choose assumptions which lead to various decompositions and interpretations. For convenience, it is assumed that the phase reference \phi_{k}=0, where possible.
6.1 Periodic Triangle Waveform Example
As an example of HSA, we consider the triangle waveform from Part I (24) and the sinusoidal FM signal.
6.1.1 Analysis assuming simple harmonic components
If we assume the component in (36) has the form given in (45), i.e. a SHC, then the AM–FM model can be obtained using Fourier analysis and the triangle waveform may be expressed with an infinite number of components
x(t)=\Re\left\{\sum\limits_{k=0}^{\infty}\frac{8A}{\pi^{2}}\frac{1}{(2k+1)^{2}% }e^{j(2k+1)\omega_{0}t}\right\}  (66) 
where the IA is given by
a_{k}(t)=\frac{8A}{\pi^{2}}\frac{1}{(2k+1)^{2}},  (67) 
the IF is given by
\omega_{k}(t)=(2k+1)\omega_{0},  (68) 
and \omega_{0} is a constant—conventionally known as the fundamental frequency.
6.1.2 Analysis assuming a single AM–FM component with harmonic correspondence
The solution in (66) may be collapsed into a single, wideband AM–FM component (K=1) by rewriting as in Part I, in terms of the AM–FM model as
x(t)=\Re\left\{a_{0}(t)e^{j[\omega_{0}t+M_{0}(t)]}\right\}  (69) 
where the IA a_{0}(t) is given by \hat{\rho}(t) in (26) and the IF \omega(t) is given by \hat{\Omega}(t) in (27).
6.1.3 Analysis assuming a single AM component
If we assume a single component (K=1) with constant frequency \omega_{0}(t)=\omega_{0}, the triangle waveform may be expressed as in Part I, in terms of the AM–FM model as
x(t)=\Re\left\{a_{0}(t)e^{j\omega_{0}t}\right\}  (70) 
where the IA a_{0}(t) is given by \hat{\rho}(t) in (30).
6.1.4 Analysis assuming a single FM component
If we assume a single component (K=1) with constant amplitude a_{0}(t)=A, the triangle waveform may be expressed as in Part I, in terms of the AM–FM model as
x(t)=\Re\left\{Ae^{j\left[\omega_{0}t+M_{0}(t)\right]}\right\}  (71) 
where the IF \omega_{0}(t) is given by \hat{\Omega}(t) in (32).
6.2 Sinusoidal FM Example
A sinusoidal FM signal with carrier frequency \omega_{c}, modulation frequency \omega_{m}, and constant B, is expressed as
x(t)=\Re\left\{e^{j\left[\omega_{c}t+B\sin(\omega_{m}t)\right]}\right\}.  (72) 
6.2.1 Analysis assuming a single FM component
We recognize that (72) is already in the form of the AM–FM model, if we assume a single component (K=1) with IA
a_{0}(t)=1.  (73) 
This leads to the IF
\omega_{0}(t)=\omega_{c}+\frac{d}{dt}B\sin(\omega_{m}t).  (74) 
6.2.2 Analysis assuming simple harmonic components
Alternatively, the sinusoidal FM signal can be expressed in terms of the AM–FM model as
x(t)=\Re\left\{\sum\limits_{k=\infty}^{\infty}J_{k}(2\pi B/\omega_{m})e^{j[(% \omega_{c}+k\omega_{m})t+\phi_{k}]}\right\}  (75) 
where J_{k}(\cdot) denotes the kthorder Bessel function of the first kind titchmarsh (). This yields an infinite number of components with IA given by
a_{k}(t)=J_{k}(2\pi B/\omega_{m})  (76) 
and IF given by
\omega_{k}(t)=\omega_{c}+k\omega_{m}.  (77) 
This is an example of a signal class where the assumption of nonnegative IF components in (41) must be relaxed, as a consequence of the summation on k becoming doublesided. For most practical choices of B and \omega_{m}, the J_{k}(2\pi B/\omega_{m}) associated with the negative frequencies will be vanishingly small, therefore the assumption of nonnegative IF, is for practical purposes, true.
6.3 Remarks
As demonstrated above, we obtain different parameterizations of the same signal by changing assumptions made during analysis. Some parameterizations, such as those obtained using Fourier analysis and the HT, are closely related due to similar assumptions. However, there exist many parameterizations each with different interpretations and without other information, we cannot say that one particular parameterization is any better or worse than any other. This highlights the importance of considering the hidden assumptions when using any particular method (FT, HT, STFT, HHT, etc.) and interpreting meaning from the parameters. If assumptions in analysis do not match the underlying signal model, erroneous interpretations may be made.
The previous solutions make use of different quadratures each corresponding to different signal models. This demonstrates the nonuniqueness of the quadrature for the AM–FM component, and hence our view of the quadrature as a free parameter. Also, note that the solutions in (66), (69), (70), etc., each utilize a complex extension that is implied by the particular assumptions made, rather than utilizing an extension based on a single rigorouslydefined procedure, like the HT. In fact, for any practicallychosen real signals x(t) and y(t), there exist assumptions leading to instantaneous parameters in which y(t) can be viewed as the quadrature of x(t).
The AM–FM model leads to exact solutions for a(t) and \omega(t), so it might appear to violate the uncertainty principle, i.e. exceed the lower limit of the timebandwidth product. However, when AM–FM modeling is viewed as a quantum mechanics problem, our casting of the problem is fundamentally different than Gabor’s. A comparison of the formal correspondences between quantum mechanics and shorttime Fourier analysis and quantum mechanics and Hilbert spectral analysis is given in Table 2. In Gabor’s casting, time and frequency are the complementary variables, i.e. (position, momentum)\leftrightarrow(t, \omega) cohen1995time (); loughlin2004uncertainty () while for the AM–FM model, the real and quadrature signals are the complementary variables, i.e. (position, momentum)\leftrightarrow(x(t), y(t)). Therefore, all uncertainty arises because only the real signal is observed. We believe that our casting is not only more appropriate, but also more useful and powerful.
7 Proposed Visualization of the Hilbert Spectrum
The ability to visualize and interpret model parameters is key to the adoption of any analysis method. Often complex AM–FM signals are plotted as a series of real 2D plots, i.e. s_{k}(t) vs. t, which for AM–FM components, would provide little insight into the underlying signal model. Alternatively, Argand diagrams may be used for AM–FM signal visualization, however, drawbacks include the quadrature signal possibly having no assigned meaning and the revolution rate not intuitively displayed.
A more appropriate visualization plots the Hilbert spectrum in its entirety. Unfortunately, plots of the Hilbert spectrum are often crudely discretized because a clear distinction between instantaneous parameters and spectral parameters is not made huang1998empirical (). We propose a method for visualizing the Hilbert spectrum, which is both complete, intuitive, and avoids the coarse discretization. The proposed visualization can be considered a (pseudo) phasespace plot because every degree of freedom or parameter of each AM–FM component is represented.
By plotting \omega_{k}(t) vs. s_{k}(t) vs. t as a line in a 3D space and coloring the line with respect to a_{k}(t) for each component, the simultaneous visualization of multiple parameters for each component is possible. Further, orthographic projections yield common plots: the timereal plane (the real signal waveforms), the timefrequency plane (Hilbert spectrum), and the realfrequency plane (analogous to the Fourier magnitude spectrum). We have found it beneficial to interpret each component as an “illuminated” particle moving in 3D space. Each particle’s motion in time and frequency is governed by \psi_{k}(t). The proposed method allows one to visualize the assumed underlying signal model.
To illustrate visualization of the Hilbert spectrum, we plot the examples in Section 6. However, the sophisticated nature of the proposed 3D visualization of the Hilbert spectrum is not wellaccommodated with paper media, and as a result we provide associated matlab functions and additional visualizations online HSAlink (). It is our preference that Hilbert spectra not be visualized with paper media. In the plots, we have utilized a perceptuallymotivated colormap in order to improve interpretation over other colormaps borland2007rainbow (); NiccoliPMK ().
For the triangle waveform example, Figs. 8(a),(c),(e) illustrate the Hilbert spectrum plots for the assumptions of SHC, single AM–FM component with HC, and single AM component, respectively. Fig. 8(a) shows the first three SHCs [constant amplitude and constant frequency in (67)] at integer multiples of a fundamental frequency, 50\pi rads/s. Fig. 8(c) shows the single AM–FM component with HC, where we see line color variation indicating a timevarying IA in (26) and a clear timevarying IF in (27). Fig. 8(e) shows the single AM component, where we see a constant IF and color variation indicating the timevarying IA in (30). Figs. 8(b),(d),(f) show the corresponding timefrequency planes of Figs. 8(a),(c),(e) by projecting out the s_{k}(t) dimension.
For the sinusoidal FM example, Figs. 9(a) and (c) illustrate the Hilbert spectrum plots for the assumptions of a single FM component and SHCs respectively. Fig. 9(a) shows the constant IA in (73) indicated by a constant line color and time varying IF in (74). Fig. 9(c) shows SHCs at integer multiples of a fundamental frequency 4\pi rads/s where each component has constant IA in (76). Figs. 9(b) and (d) show the corresponding timefrequency planes of Figs. 9(a) and (c) by projecting out the s_{k}(t) dimension.
8 Summary
In Part II of this paper, we have presented theory for the Hilbert spectrum as a generalized LSA problem. In the general problem, we seek a representation of z(t) consisting of a superposition of latent signals, i.e. multicomponent model consisting of a superposition of complex AM–FM components. We have used the AM–FM model to define the Hilbert spectrum as parameterized by a set of IA/IF pairs (and phase references) each associated with the components. In the LSA problem, relaxation of the HC condition allowed freedom in the choice of the signal’s quadrature and admitted many solutions for the latent signal. In the HSA problem, relaxation of the HC condition allows freedom in the choice of each component’s quadrature thereby allowing even greater freedom in the construction of the signal’s model. We illustrated how assumptions on the form of the component, i.e. constant amplitude and constant frequency, timevarying amplitude and constant frequency, constant amplitude and timevarying frequency, etc., lead to familiar specializations of the model.
We discussed the frequency domain view of LSA including the implications of relaxing HC. From there we discussed the Hilbert spectrum view of LSA, specifically the implications of using a general AM–FM component with relaxed HC and assumed positive IF. Closed form HSA examples were provided in which we assume various forms of the AMFM component and determined the unique corresponding instantaneous parameterization in terms of IA and IF. A novel 3D visualization of the Hilbert spectrum was proposed by plotting \omega(t) vs. s(t) vs. t and coloring with respect to a(t).
Finally, we discussed the analogy of the Hilbert spectrum to quantum mechanics in which our casting of timefrequency analysis is fundamentally different than Gabor’s casting, where the uncertainty in our framework is in the quadrature signal and not the frequency variable. This provides a new and powerful framework for nonstationary signal analysis.
Part III: Numerical HSA Assuming Intrinsic Mode Functions
\pdfstringdef\BKM@titlePart III: Numerical HSA Assuming Intrinsic Mode Functions
[T]here are some crucial restrictions of the Fourier spectral analysis; the system must be linear; and the data must be strictly periodic or stationary; otherwise, the resulting spectrum will make little physical sense.—Huang huang1998empirical ()
In Part II, the AM–FM signal model and HSA theory was presented, analysis was performed on example signals, and a visualization of the Hilbert spectrum was proposed. In Part III, numerical algorithms are given for estimating the IAs and IFs of the components in the AM–FM model where now the component is assumed to be an Intrinsic Mode Function (IMF). These algorithms first decompose the signal into IMFs using an improved version of Huang’s original Empirical Mode Decomposition (EMD) algorithm and second, demodulate the IMFs to obtain the instantaneous parameters. Unlike previous studies, close attention is paid to the assumptions made in the definition of the IMF which are carried forward to the demodulation step, thereby avoiding any ambiguity associated with obtaining the instantaneous parameters. We begin with a comprehensive review of EMD and several variations, and propose an algorithm for the computation of the Hilbert spectrum assuming IMF components. It is important to note that while IMFs can be considered latent AM–FM components, there are other classes of AM–FM components that are not IMFs as illustrated in Fig. 5. Examples using the proposed algorithm are provided that highlight alternative decompositions compared to traditional Fourier analysis and demonstrates the advantages of using the HSA framework.
1 Introduction
In Part II of this paper, we defined the Hilbert spectrum using the AM–FM model in (35) where the AM–FM component is defined in (36) and parameterized by the IA a_{k}(t), IF \omega_{k}(t), and phase reference \phi_{k}. We assume a real observation, x(t) of a latent signal, z(t).
In Part III of this paper, we turn our attention to numerical computation of instantaneous parameters of the AM–FM model. Many methods for computing the parameters of an AM–FM model already exist, each corresponding to a specific set of assumptions. As discussed in Part II, a FT corresponds to the assumption of SHCs and a STFT corresponds to the assumption of AM components. Other AM–FM models have been proposed with alternative assumptions, however, common to these are restrictions which limit the utility of the model.
Practical estimation of the instantaneous parameters of the AM–FM model in (35) is a twostep process. First, the signal must be decomposed into a set of latent AM–FM components and second, the instantaneous parameters \{a_{k}(t),\omega_{k}(t)\} of each component must be estimated. To date, the most flexible decomposition method for AM–FM modeling is Huang’s EMD and its variations huang1998empirical (). In contrast to most timefrequency analysis methods, EMD makes no assumption of constant amplitude or frequency on the component and has mild bandlimiting assumptions. In huang1998empirical (), Huang proposed the original EMD algorithm which sequentially determines a set of AM–FM components, i.e. IMFs via an iterative sifting algorithm. The Ensemble Empirical Mode Decomposition (EEMD) wu2009ensemble () and tone masking DeeringMasking () introduced ensemble averaging in order to address the mode mixing problem. The complete EEMD was proposed to address some of the undesirable features of EEMD by averaging at the componentlevel as each component is estimated rather than averaging at the conclusion of EMD torres2011complete (). Several improvements to the sifting algorithm have also been proposed including those by Rato rato2008hht (). These improvements to the original EMD will be more fully described in this part.
In the second step, instantaneous parameters must be estimated through demodulation of the IMFs. Despite the numerous attempts to demodulate IMFs, most of the proposed methods have fallen short because the assumptions made during decomposition have not been maintained during demodulation. In this paper, we identify and pay close attention to these assumptions as we develop a demodulation technique which adheres to the assumptions. We point out that Rato proposed an AM–FM demodulation procedure in which the IA estimation was consistent with the assumptions but the IF estimation was not rato2008hht (). Also, Huang has examined numerous demodulation methods, including an iterative normalization to obtain an FM signal which is then demodulated to estimate IF using an \arctan approach that unfortunately suffers from numerical instability huang2009instantaneous (). The iterative normalization method is consistent with the decomposition assumptions. Utilizing both Rato’s AM estimation and Huang’s iterative normalization procedure, we propose a mathematically equivalent method to get IF from the FM signal that is more numerically stable. We then incorporate the proposed demodulation and EMD into a single HSA–IMF algorithm which gives very good estimates for the IA and IF parameters of the AM–FM model.
Part III of this paper is organized as follows. In Section 2, we review the EMD algorithm and its variations and improvements by other researchers. By carefully noting the assumptions made in EMD, we propose demodulation of the resulting IMFs and present a complete HSA–IMF algorithm composed of EMD and demodulation. In Section 3, we provide two sets of examples using the HSA–IMF algorithm. The first set of examples uses synthetic signals which illustrate how EMD can estimate the parameters of the underlying, assumed AM–FM signal model. The second set of examples uses realworld audio signals including speech which further illustrates AM–FM model parameterizations with alternate interpretations. These examples clearly demonstrate the alternate decomposition and interpretation offered by HSA. In Section 4, we discuss other issues related to HSA and in Section 5, we discuss future research. In Section 6 we summarize Part III.
2 Empirical Mode Decomposition
EMD consists of an iterative procedure for decomposing a signal into a set of IMFs, \{\varphi_{k}(t)\} huang1998empirical (). In huang1998empirical (), the definition of an IMF is any signal that satisfies two conditions:

In the whole signal segment, the number of extrema and the number of zero crossings must be either equal or differ at most by one.

At any point the mean value of the envelope, defined by the local maxima and the envelope defined by the local minima, is zero.
In the context of LSA and AM–FM modeling, an IMF can be thought of as a latent component,
\varphi_{k}(t)=\Re\{\psi_{k}(t)\}.  (78) 
Thus we view an IMF as an inherently complexvalued component, which is not the conventional interpretation. Furthermore, we argue that the definition of an IMF forces a unique quadrature that in general does not equal \mathcal{H}\{\varphi_{k}(t)\} as will be discussed in Section 2.4.
Unlike traditional methods in timefrequency analysis, EMD is defined by an algorithm rather than transform theory, which has both advantages and disadvantages. One advantage is that EMD utilizes a component that is far less restrictive than other timefrequency methods cohen1995time (). One disadvantage is that EMD is understood primarily through experimentation DeeringMasking (). Empirical experiments using white noise have shown EMD to act as a dyadic filter bank wu2004study (); flandrin2004empirical (); FlandrinHuang2005 (); wu2009ensemble (); mandic2011filter (). Using fractional Gaussian noise as a model for broadband noise, it has been shown that the builtin adaptivity of EMD makes it behave spontaneously as a ‘waveletlike’ filter, i.e. the result of sifting is a ‘detail’ and a ‘trend’ flandrin2004empirical ().
Efforts have also been made to replace the sifting algorithm with alternate formulations which are more mathematically grounded, such as techniques based on optimization meignen2007new (); kopsinis2008investigation (); hou2013data (), machine learning looney2008machine (), PDEs delechelle2005empirical (); sharpley2006analysis (); vatchev2008decomposition (); diop2009pde (); el2010analysis (); el2013pde (), and Fourier analysis niang2010spectral (). We also point out other research on the EMD algorithm including the multivariate EMD damerval2005fast (); rilling2007bivariate (); wu2009multi (); linderhed2009image (); mandic2011filter (); nunes2009empirical (); ThePowerAdaptive () and the complex EMD tanaka2007complex (). As a final note, not all variations of EMD utilize the IMF component or the proposed AM–FM model and thus cannot be considered as a form of HSA. Examples include the Hilbert variational decomposition feldman2006hvd (); feldman2008compare (); feldman2011hilbert (), the timedependent intrinsic correlation chen2010time (), and synchrosqueezed wavelet transforms daubechies2011synchrosqueezed (); wu2011one (); chui2015signal ().
2.1 The Original EMD Algorithm
The original EMD and sifting algorithms proposed by Huang huang1998empirical () are listed in Algorithms 1 and 2. The purpose of the sifting algorithm is to iteratively identify and remove the trend from the signal, effectively acting as a high pass filter. Step 5 of Algorithm 1 removes the high frequency component \varphi_{k}(t), estimated during the sifting process. The process is then repeated to remove additional IMFs from the signal if they exist. The resulting decomposition is complete and sparse huang1998empirical (); hou2011adaptive (); hou2013data (). EMD is formulated with continuoustime signals, however, in practice, EMD is applied to discretetime signals which may result in errant decompositions. The effects of sampling in the context of EMD have been considered by Rilling and it is generally recommended to oversample but not resample before application of EMD, so that EMD effectively behaves like a continuous operator rilling2006sampling ().
2.2 Improving the Sifting Algorithm
If EMD is viewed as an AM–FM decomposition technique, then the sifting algorithm is an iterative way of removing the asymmetry between the upper and lower envelopes in order to transform the input r(t) into an IMF rato2008hht (). By doing so, low frequency content is discarded at every sifting iteration, effectively making the sifting algorithm behave as a high frequency filter or high frequency component tracker. Due to the doubly iterative nature of EMD and termination conditions, numerical imprecision and differing implementations can lead to very different IMFs. To achieve consistency when using EMD, Rato proposes the following constraints rato2008hht ():

IMFs should scale with the signal.

Any signal bias, should only be reflected in the trend.

The EMD of an IMF should be the IMF itself.^{6}^{6}6As Rato points out, this constraint may be relaxed in some cases rato2008hht ().

Timereversal of the signal should timereverse the IMFs.
Several improvements have been made to the sifting algorithm (Algorithm 2) to improve the decomposition accuracy rato2005Modified (); rato2008hht ():
Although several stopping criteria have been proposed huang2003confidence (); rilling2003algorithms (); rato2008hht (), Rato suggested the use of a resolution factor which is the ratio between the energy of the signal at the beginning of sifting r(t) and the energy of the average of the envelopes e(t). If this ratio increases above a predetermined threshold, then the IMF computation terminates. We have found that Rato’s use of a parabolic interpolator rato2008hht () to identify the extrema in Steps 4 and 5 is convenient because it uses as few as three samples and it avoids the classificatory function to determine if a particular sample is a maxima, minima, or neither. End effects appear due to the fact that a given interpolator may not be a good extrapolator rato2008hht (). In order to deal with this problem, Rato suggested to insert artificial minima and maxima in order to control the behavior of the interpolator. In addition, the mean envelope removal is scaled in Step 9,
r(t)\leftarrow r(t)\alpha e(t)  (79) 
where the stepsize 0<\alpha\leq 1 increases the number of sifting iterations but improves stability and robustness of the resulting IMFs rato2008hht ().
We illustrate a single iteration of the sifting algorithm in Fig. 10. In Fig. 10, we show two unknown components, \varphi_{0}(t) ( ) and \varphi_{1}(t) ( ) and in Fig. 10, we show the signal under analysis r(t)=\varphi_{0}(t)+\varphi_{1}(t) ( ) which is the input to the sifting algorithm. Figs. 10 and (d) show the location and interpolation of the extrema as given in Steps 4  7. Interpolations lead to estimates of the upper envelope u(t) ( ) and lower envelope l(t) ( ) of r(t). The average of the upper and lower envelopes e(t)\approx\varphi_{1}(t) is shown in Fig. 10. At the end of the first iteration of sifting, r(t)e(t)\approx\varphi_{0}(t) and is shown in Fig. 10.
2.3 Improving the EMD Algorithm
The major problem in the EMD algorithm is mode mixing, which is defined as a single IMF either consisting of components of disparate scales or components of similar scale residing in the same IMFs wu2009ensemble (). Mode mixing is a consequence of signal intermittency, or more specifically relative component intermittency. As a result, the particular component(s) tracked by the sifting algorithm in a particular IMF at any instant may change as intermittent components begin or end DeeringMasking (). This is illustrated in Fig. 11, where in the left part of the figure, we show components of disparate scales being in the same IMF denoted by ⚪, while in the center part of the figure we show two components of similar scale in the same IMF denoted by ◻. The ability of EMD to resolve two components considering both the relative IAs and IFs of components, was examined and quantified by Rilling rilling2008one (). However as was noted, resolving closelyspaced components may not be the ultimate goal, provided that the decomposition is suitably matched to some meaningful interpretation rilling2008one ().
Two commonly used methods of mitigating mode mixing are EEMD wu2009ensemble () and tone masking DeeringMasking (). EEMD (given in Algorithm 3, where \mathbb{E}[\cdot] denotes the expectation) utilizes zeromean white noise, w^{(i)}(t) to perturb the signal so a component may be tracked properly over an ensemble average. As the illustration in Fig. 11 shows, noise can be used to assist the sifting algorithm. Inserting noise with high enough power gives the sifting algorithm something to track when the highest frequency component is intermittent, then vanishes in the ensemble average. Although injecting noise can help to track components properly, a carefully designed masking signal can result in better performance. A situation where a carefully designed masking signal may be beneficial is illustrated in Fig. 11 denoted by ◻.
EEMD is not without its disadvantages as it is more computationally complex, loses the perfect reconstruction property, propagates IMF estimation error, results in inconsistent numbers of IMFs across the trials, and the resulting set of averaged IMFs \{\bar{\varphi}_{k}(t)\} are not necessarily IMFs torres2011complete (). Torres proposed the “complete EEMD” given in Algorithm 4 to address some of these issues torres2011complete (). Complete EEMD defines a procedure \mathrm{EMD}_{k}(\cdot) which returns the kth IMF obtained using EMD torres2011complete (). This method of EEMD requires fewer sifting iterations, a smaller ensemble size, and recovers the completeness property of the original EMD algorithm to within the numerical precision of the computer torres2011complete ().
Rather than using a noise signal, a deterministic signal v(t) can also be used as a perturbation and then removed after IMF estimation. The “tone masking” technique is given in Algorithm 5 and can be used, for example, in place of Step 4 in Algorithm 1 DeeringMasking (). Advanced forms of tone masking have been proposed and are termed SignalAssisted EMD (SAEMD) FlandrinHuang2005 (); senroy2007two (); guanlei2009time (); li2011sinusoidal (). These methods have additional advantages in their ability to track closelyspaced components, but require careful selection of the masking signal.
2.4 IMF Demodulation
In order to obtain the IA/IF parameters in the AM–FM model, we are required to demodulate the IMFs returned by EMD. One approach, used by Huang, is to apply the HT to each IMF in order to obtain estimates of IA and IF. This analysis is often referred to as the HilbertHuang Transform (HHT) huang1998empirical (). Using the HT for IMF demodulation is inappropriate because as argued, the HT assumes SHCs and HC—which is not true since an IMF is in general an AM–FM component without HC. A similar argument can be made of many of the other demodulation methods that have been considered huang2009instantaneous (). A second approach, as advocated in Part II, is to let the assumed form of the AM–FM component imply a complex extension. When viewed in the context of LSA, the definition of the IMF given in Section 2 forces a unique complex extension to the IMF—justifying our view of the IMF as a latent component specified by a real signal.
Inherent in C2 of the IMF definition, are the following assumptions:

a(t_{p})=s(t_{p}), where s(t_{p}) are the extrema of s(t).

a(t),t\notin\{t_{p}\} is inferred by cubic spline interpolation.
The first assumption, A1 can be viewed as quadrature forced to zero at \{t_{p}\}. To see this, note that from (36), we have
a(t)=\mathchoice{{\hbox{$\displaystyle\sqrt{s^{2}(t)+\sigma^{2}(t)\,}$}% \lower 0.4pt\hbox{\vrule height 8.959863pt depth 7.167891pt width 1px}}}{{% \hbox{$\textstyle\sqrt{s^{2}(t)+\sigma^{2}(t)\,}$}\lower 0.4pt\hbox{\vrule hei% ght 8.959863pt depth 7.167891pt width 1px}}}{{\hbox{$\scriptstyle\sqrt{s^{2}(% t)+\sigma^{2}(t)\,}$}\lower 0.4pt\hbox{\vrule height 6.999893pt depth 5.59991% 5pt width 1px}}}{{\hbox{$\scriptscriptstyle\sqrt{s^{2}(t)+\sigma^{2}(t)\,}$}% \lower 0.4pt\hbox{\vrule height 6.999893pt depth 5.599915pt width 1px}}}  (80) 
and thus, \sigma(t_{p})=0 and a(t_{p})=s(t_{p}). A1 also implies nonnegativity of a(t_{p}), however, nonnegativity of a(t) is not guaranteed for all t due to the interpolation.
The second assumption, A2 can be viewed as a relative bandlimiting condition on a(t) controlled by two factors: the choice of interpolator and the density of extrema. To see this, note that a(t) between extrema of s(t) is defined by an interpolator. This implies that a(t) will only have as much variation between extrema as the interpolator allows. However, with dense extrema (high IF) much less restrictive constraints have been imposed on a(t) than if extrema are sparse (low IF)—thus our view as a relative bandlimiting condition on a(t). Also note, using an interpolator other than the cubic spline is likely to change the IA, effectively changing the resulting IMFs rato2008hht (); hawley2010some ().
In the IMF definition, C1 requires that the number of extrema and the number of zero crossings must be either equal or differ at most by one. A general AM–FM component may not satisfy this condition, for example as a result of sign reversal of \omega(t) or in cases where A1 and A2 are not satisfied. Assumption of positive IF on the AM–FM component in Part II, eliminates the possibility of sign reversal of the IF but may not be true for all signal classes.
Rato proposed an AM demodulation approach given in Algorithm 6, that is consistent with the decomposition assumptions in Algorithm 2 rato2008hht (). Starting with an IMF estimate \hat{\varphi}(t), we obtain an estimate for IA \hat{a}(t) which can then be used to estimate the IF via the FM signal
\hat{s}_{\mathrm{FM}}(t)=\frac{\hat{\varphi}(t)}{\hat{a}(t)}.  (81) 
However, the estimate in (81) may result in \hat{s}_{\mathrm{FM}}(t)>1. Thus, Huang proposed an iterative normalization procedure, given in Algorithm 7, which removes the AM from the signal to obtain a more accurate \hat{s}_{\mathrm{FM}}(t) huang2009instantaneous (). Although the iterative procedure improves demodulation accuracy, we noticed that it can be susceptible to oscillating artifacts introduced by overfitting of the cubic spline interpolator. As a result, we have found that these artifacts can be minimized by limiting the number of iterations to three.
We can directly obtain two estimates of the IF, \pm\hat{\omega}(t) by substituting a(t)=1 in (36b), s(t)=\hat{s}_{\mathrm{FM}}(t) in (36c), and then computing (41). Direct FM demodulation is not straightforward because different approaches exist for obtaining the IF from \hat{s}_{\mathrm{FM}}(t), that although are mathematically equivalent, may differ in numerical stability huang2009instantaneous ().
We propose a FM demodulation method to address the numerical stability issues. We begin by estimating the quadrature in (36c) as
\hat{\sigma}_{\mathrm{FM}}(t)=\text{sgn}\left[\frac{d}{dt}\hat{s}_{\mathrm{FM% }}(t)\right]\mathchoice{{\hbox{$\displaystyle\sqrt{1^{2}\hat{s}_{\mathrm{FM}}% ^{2}(t)\,}$}\lower 0.4pt\hbox{\vrule height 8.959863pt depth 7.167891pt width% 1px}}}{{\hbox{$\textstyle\sqrt{1^{2}\hat{s}_{\mathrm{FM}}^{2}(t)\,}$}\lower 0% .4pt\hbox{\vrule height 8.959863pt depth 7.167891pt width 1px}}}{{\hbox{$% \scriptstyle\sqrt{1^{2}\hat{s}_{\mathrm{FM}}^{2}(t)\,}$}\lower 0.4pt\hbox{% \vrule height 8.399872pt depth 6.719897pt width 1px}}}{{\hbox{$% \scriptscriptstyle\sqrt{1^{2}\hat{s}_{\mathrm{FM}}^{2}(t)\,}$}\lower 0.4pt% \hbox{\vrule height 8.399872pt depth 6.719897pt width 1px}}}  (82) 
where \text{sgn}\left[\frac{d}{dt}\hat{s}_{\mathrm{FM}}(t)\right] is required to obtain an appropriate four quadrant estimate with assumed positive IF. In Fig. 12 , we have provided an illustration (animated in the web version) to assist the reader with understanding of how we arrive at (82). The computationally unstable points, \{t_{0}\} occur where \hat{\sigma}_{\mathrm{FM}}(t_{0})=0, thus we can replace a small range around these points (t_{0}\epsilon,t_{0}+\epsilon) with interpolated values. Then
\hat{\theta}(t)=\text{arg}\left[\hat{s}_{\mathrm{FM}}(t)+j\hat{\sigma}_{% \mathrm{FM}}(t)\right]  (83) 
and the IF is obtained with (41). We can optionally smooth the resulting IF estimate. Our complete IMF demodulation algorithm is listed in Algorithm 8.
We note to the reader that the proposed IMF demodulation procedure, which does not involve the HT or FT, can in some cases return the same results but in general this is not true. For example, demodulation of \cos(\omega_{0}t) with either the IMF or HT demodulation returns a SHC with HC. On the other hand, HT demodulation of the triangle waveform returns a single AMFM component with HC as given in (27), while IMF demodulation returns IF given by (32) and constant IA.
2.5 Proposed Algorithm for HSA Assuming IMFs
The various recommendations and improvements listed this section appear to have never been combined and investigated together. Therefore, we propose to incorporate the most desirable features of the complete EEMD and tone masking for addressing the modemixing problem, Rato’s improvements to the sifting algorithm, and our proposed demodulation method, into a single algorithm given in Algorithm 9.
Although the masking signal v^{(i,k)}(t) is usually a carefully chosen realvalued AM–FM signal senroy2007two (); guanlei2009time (), we choose v^{(i,k)}(t) as a white, Gaussian, lowpassfiltered noise with cutoff frequency below the amplitudeweighted IF strom1977amplitude () of the previously found IMF. This replaces the need to use sifted white noise in complete EEMD, however, we do allow for more sophisticated approaches for choosing the masking signal. Note that this introduces a feedback loop between decomposition and demodulation. This changes EMD from simply a decomposition method, to a full HSA method because the IA and IF estimates of the kth component are inherently computed and then used to design the masking signal for the (k+1)th component.
3 Examples using the HSA–IMF Algorithm
In this section, we demonstrate the HSA–IMF algorithm for signal analysis and compare to conventional STFT using both synthetic and realworld signals. In these examples, we use the proposed visualization method from Part II orthographically projected onto the timefrequency plane, for plotting the instantaneous parameters resulting from HSA. We also plot the STFT magnitude (STFTM) using the same colormap to facilitate comparisons NiccoliPMK (). The STFT is computed using a 4096point hamming window with a 1 sample frame advance. We have provided matlab functions for these examples HSAlink ().
3.1 Synthetic Signals
We provide two examples using synthetic signals with known underlying signal models. The synthetic signal examples demonstrate the proposed algorithm on a single component in a noiseless environment. In this case, mode mixing is not a problem and we initialize I=1, \alpha=0.95, and \beta_{k}=0 in Algorithm 9.
In the first example, we analyze a signal with a slowvarying AM and fastvarying FM given by
x(t)=\Re\left\{a(t)\exp\left[j\left(6000\pi t+\int\limits_{\infty}^{t}m(\tau)% d\tau\right)\right]\right\}  (84) 
where the IA is
a(t)=e^{(t0.5)^{2}/25}  (85) 
and the FM message is
m(t)=250\sin(140\pi t)+2000[\exp(4t)1].  (86) 
Fig. 13 shows the STFTM where we see classic harmonic structure resulting from the inherent assumption of SHCs despite (84) consisting of only a single component. Fig. 13 shows a single component in the Hilbert spectrum with a fast FM oscillation consistent with the underlying signal model in (84). The IA in (85), consisting of a Gaussian envelope centered at t=0.5, is visible as color variation in Fig. 13. The FM message in (86) consists of a 70 Hz oscillation superimposed on a decaying exponential; this FM message is offset by a 3000 Hz carrier. The FM message is reflected by the vertical behavior in Fig. 13 where we see the IF sweeping from 3000 Hz to 1000 Hz with a 70 Hz oscillation. Due to more appropriate assumptions of an underlying component, this example shows that the HSA–IMF parametrization can more closely match the underlying signal model than traditional methods.
In the second example, we analyze a signal with a fastvarying AM and slowvarying FM given by
x(t)=\Re\left\{a(t)\exp\left[j\left(1000\pi t+\int\limits_{\infty}^{t}m(\tau)% d\tau\right)\right]\right\}  (87) 
where the IA is
a(t)=\dfrac{1}{2}+\dfrac{1}{3}\sin(100\pi t)+\dfrac{1}{5}\sin(200\pi t)  (88) 
and the FM message is
m(t)=150\sin(2\pi t).  (89) 
Fig. 14 shows the STFTM where we again see classic harmonic structure resulting from the inherent assumption of SHCs despite (87) consisting of only a single component. Fig. 14 shows a single component in the Hilbert spectrum with a fast AM oscillation, captured by color variation in the plot line, consistent with the underlying signal model in (87).
As described in Part II, signals composed of narrowband components can have similar Fourier and Hilbert spectra, however, signals composed of wideband components can have spectra which are quite different. Fourier analysis of wideband signals may result in multiple ridges in terms of spectral frequencies, which may be further decomposed into narrowband components. However, for many realworld signals, a wideband component, such as the AM–FM component in (36), is the more appropriate choice rather than multiple, narrowband components. The appearance of structure in the Fourier spectrum may indicate the presence of a wideband component(s) in the signal as demonstrated in Example 1, Fig. 13 (fastvarying IF) and in Example 2, Fig. 14 (fastvarying IA).
Both decompositions, i.e. STFT and HSA–IMF are equally valid models for the real signal x(t)—the resulting decompositions simply correspond to the different assumptions of the underlying components, i.e. constantfrequency components and IMF components, respectively. The complex extensions assumed in STFT and implied in HSA–IMF are fundamentally different and because of this, the quadrature signal y(t) can be different for the two analysis methods. This ultimately results in a different z(t) for each analysis and consequently different instantaneous parameterizations despite both models produce the same x(t) in (2).
3.2 RealWorld Signals
We provide two examples using realworld audio signals. In Algorithm 9, we initialize I=200, \alpha=0.95, and through experimentation, \beta_{k}=4 for all k. In addition, we have applied a 1 ms movingaverage filter to smooth the IF estimate.
In the first example, we analyze a recording (f_{s}=22.050 kHz) of a single note played on a cello. Fig. 15 shows the STFTM where we again see classic harmonic structure resulting from the inherent assumption of SHCs. We note the fundamental frequency is approximately 67 Hz (15 spectral lines evenlyspaced over 1000 Hz) and two dominant spectral lines at the second harmonic (\sim133 Hz) and at the fifth harmonic (\sim333 Hz). In addition, there is a brief dominant spectral line at t=2 s corresponding to the ninth harmonic at \sim600 Hz.
Fig. 15 plots the five components returned by the HSA–IMF algorithm, where we see three dominant components. The lower two components, range in IF from \sim120140 Hz and from \sim300360 Hz corresponding to the dominant spectral lines in the Fourier spectrum. The upper component also exhibits significant energy at t=2 at \sim550750 Hz corresponding to the brief dominant spectral line, i.e. ninth harmonic in the Fourier spectrum.
In the second example, we analyze a recording (f_{s}=44.1 kHz) of the word “shoot.” Fig. 16 shows the STFTM where we see the spectral energy of the fricative “SH” over 0\leq t\leq 0.15 s, scattered over the range 0 to 8 kHz. The spectral energy for the vowel “UW” over 0.15\leq t\leq 0.25, is concentrated at a fundamental of \sim230 Hz. The spectral energy for the stop “T” over 0.37\leq t\leq 0.4 s, is very weak and spread across the band and hence not visible in the plot.
Fig. 16 plots the five components returned by the HSA–IMF algorithm. The “SH” fricative appears in three components with IF ranges of \sim60007000 Hz (zeroth component), \sim25005000 Hz (first component), and \sim10002500 Hz (second component) but is mostly captured in the second component with quickly varying AM and FM. The vowel “UW” is clearly captured in a single component near 230 Hz exhibiting some FM variation conjectured to be natural jitter. Unlike in the Fourier spectrum, the stop “T” is clearly captured in the Hilbert spectrum by the first component near t=0.4.
4 Discussion
4.1 Resolving CloselySpaced Components
EMD has been criticized for its inability to resolve closelyspaced components and there have been numerous studies and analyzes on the resolving ability DeeringMasking (); rilling2008one (); feldman2009analytical (). Rilling has investigated EMD’s (Algorithm 1) ability to resolve two tones as a function of a relative amplitude parameter and a relative frequency spacing parameter. This analysis describes regions in this parameter space where EMD returns one or two components rilling2008one (). As Rilling noted, the goal in EMD is not resolving closelyspaced components but rather resolving components that are suitably matched to an underlying signal model or compatible with assumptions made on the signal model rilling2008one ().
As an example, consider two infinitelylong tones, \cos\left(\omega_{a}t\right) and \cos\left(\omega_{b}t\right), with \omega_{b}>\omega_{a}. We can express the sum of these tones as
\displaystyle x(t)  \displaystyle=\Re\left\{\exp\left(j\omega_{a}t\right)+\exp\left(j\omega_{b}t% \right)\right\}  (90a)  
\displaystyle=\Re\left\{2\cos\left[\left(\omega_{b}\omega_{a}\right)t/2\right% ]\exp\left[j\left(\omega_{b}+\omega_{a}\right)t/2\right]\right\}.  (90b) 
If \omega_{a} and \omega_{b} are sufficiently far apart, both Fourier analysis and EMD will resolve two SHCs as in (90a). On the other hand, if \omega_{a} and \omega_{b} are not sufficiently far apart, EMD will resolve a single IMF as in (90b). As is well known, when \omega_{a} and \omega_{b} are closelyspaced, the signal exhibits a beat effect. In the human auditory system, these closelyspaced tones are not perceived as two distinct tones but rather a single AM tone oshaughnessy (). As Deering points out, EMD may correspond to the psychoacoustics of human hearing DeeringMasking (). As Rilling points out, a decomposition into SHCs may not be an appropriate solution “…if the aim is to get a representation matched to physics (and/or perception) rather than to mathematics” rilling2008one ().
A generalized example of this beat effect was given in Example 2 in Subsection 3.1. To clarify the connection to auditory perception of beating, ignore the slowvarying FM and hence the plot lines in Fig. 13 would be horizontal and not sinusoidal. The example, when analyzed with the STFT shows five closelyspaced tones as in Fig. 13. However, if these tones are closely spaced then as noted, they may be perceived as a single tone with AM variation. This is demonstrated with an analysis using HSA–IMF and shown in Fig. 13.
A similar example regarding an FM signal, was given in Example 1 in Subsection 3.1, and is essentially the same waveform used in FM synthesis pioneered by Chowning chowning1977synthesis (). In Chowning’s work, he expressed the FM signal as a superposition of SHCs weighted by Bessel functions and showed rich spectra associated with this signal. Fig. 14 illustrates this rich spectra via the presence of multiple harmonics. In the AM–FM model, such rich spectra may be encapsulated in a single component as illustrated in Fig. 14 and as has been suggested, a decomposition into SHCs may not be an appropriate solution to describe the underlying signal model.
4.2 Remarks on Computation
The HSA–IMF algorithm contains a triplenested loop that incurs significant computation depending on the signal and parameter choices. Clearly, if the signal has many underlying components, EMD will require more computation due to the extrema searches and interpolations. Unfortunately, there is no way to predict the number of components ahead of time.
The outermost loop in Algorithm 9 Step 3 iteratively removes the IMF estimated by tone masking until termination conditions are reached. Hence, there is no way to predict the number of iterations required for termination of the loop. The middle loop, ensemble averages the IMFs returned from tone masking in Algorithm 9 Step 4. This average is controlled by a fixed number of trials, I (the choice of I is discussed in the next subsection). The inner loop results from the tone masking procedure in Algorithm 9 Step 4 calling the sifting algorithm twice, in Algorithm 5 Step 3. This in turn calls Algorithm 2, in which Step 3 iteratively estimates an IMF. Within this inner loop, the stepsize introduced in (79) and used in Step 9, also controls the speed at which termination conditions are reached. Of these loops, the innermost loop requires the most computation due to the search for extrema and interpolation. Taken together, the iterative nature of the outer and inner loops, coupled with computationally complex inner loop can require significant computation depending on the signal length and sample rate. As has been pointed out, signal oversampling is required for robust estimates of the IMFs further increasing computation. Finally, IMF demodulation in Step 5 occurs in the outermost loop and does not add significant computational burden. Much of this computation can occur in parallel waskito2010parallelizing (); chen2010gpgpu (); chang2011parallel (); waskito2011evaluation (), for example, the ensemble averaging in Algorithm 9 Step 4 and the search for extrema in Algorithm 2 Steps 4 and 5.
We have implemented Algorithm 9 in matlab where the trials are computed in parallel using a parfor loop and have timed the computation for the synthetic and realworld signal examples in this paper. Our PC consists of an eightcore AMD FX at 4.01 GHz with 32 GB RAM. The results are given in Table 3 where we see for relatively short audio signals, decomposition may require relatively large \beta and I leading to long computation times.
Example  f_{s} (kHz)  Duration (s)  \beta  I  Computation Time (s) 

Synthetic 1  44.1  1  0  1  1.3 
Synthetic 2  44.1  1  0  1  0.9 
Cello  22.05  4.21  4  200  320.1 
Speech  44.1  0.43  4  200  52.4 
4.3 HSA–IMF Robustness
The goal in HSA is to obtain an AM–FM decomposition with meaningful interpretation. Thus, proper identification of the underlying components is required from Algorithm 9. Two parameters must therefore be carefully chosen: the SNR factor \beta_{k} and the number of trials I. As a reminder, \beta_{k} weights the additive masking signal used to mitigate mode mixing and I minimizes, through ensemble averaging, the influence of the additive masking signal on the resulting IMF; both of these parameters appear in Step 4.
In our experience with realworld signals where the underlying signal model is unknown, we begin by selecting \beta_{k}=0 and I=1 and visualizing the resulting IMFs by plotting the timereal plane or timefrequency plane. Mode mixing will be evident in the timereal plane when the frequency of the waveform changes abruptly. Mode mixing will also be evident in the timefrequency plane when the IMF is similar to the illustrated IMF within the frame in Fig. 11. If mode mixing is present, we increase \beta_{k} and I. This process is repeated until reasonable IMFs are obtained keeping in mind the associated computational load for large I. Although this process for decomposition is somewhat heuristic, such refinement is typically present in all timefrequency analyzes cohen1995time (), e.g. choice of window length and type for Fourier analysis and mother wavelet selection for wavelet analysis.
The stepsize parameter \alpha introduced in (79) and used in Algorithm 2 Step 9, is of secondary importance and merely scales the trend which is removed from the signal being sifted. This scaling is used to minimize the impact of possible overshoot in the trend. In our experience, selecting \alpha=0.95 gives satisfactory performance, noting that lower values lead to additional iterations and more computation.
The termination condition in the sifting algorithm may be too restrictively chosen hence preventing convergence. In our implementation we include a maximum number of iterations, typically 50, to guarantee an exit. As a final point, IMFs with excessively large values are omitted from the ensemble.
5 Future Research in HSA–IMF
In the course of this research, several avenues for further algorithm improvement are apparent. These include: alternative masking signals, alternate interpolators, and improvements to IMF demodulation. In the HSA–IMF algorithm, we use lowpass filtered noise as the masking signal. However, for certain signal analysis problems, more sophisticated masking signals have been proposed senroy2007improved (); senroy2007two (); guanlei2009time (). Ideally, the masking signals would have properties similar to the underlying components which would likely require prior knowledge of the signal model. At the heart of EMD is the iterative estimation of the IMFs which are completely determined by interpolation. As noted earlier, the cubic spline interpolator may be susceptible to overfitting which can lead to poor IMF estimates. Alternatives to cubic spline interpolation have been investigated including Bspline and Akima interpolators riemenschneider2005b (); qin2006envelope (); chen2006b (); rato2008hht (); kopsinis2008improved (); wang2010intrinsic (); bouchikhi2012multicomponent (), however, these do not appear to offer significant improvements. Nevertheless, improvements in interpolation may lead to more robust decomposition and demodulation. Finally, IMF demodulation requires estimation of the IF which uses Huang’s iterative normalization procedure. Unfortunately, the required interpolation in the normalization procedure can result in an overfitting of the cubic spline leading to incorrect instantaneous estimates. As noted earlier, changing the cubic spline interpolator in Algorithm 8 changes the IMF. Thus a change of the interpolator in the IA estimator requires the same change in the sifting algorithm.
Although IMFs are latent AM–FM components, there are other classes of AM–FM components that are not IMFs opening up alternate possibilities. The IMF is only important due to its computability via sifting. Hence, it may be possible to define other useful AM–FM components (not defined by C1 and C2) corresponding to different assumptions on the form of the component (not A1 and A2), that ultimately lead to a replacement of the sifting algorithm as the component tracker, and thus new AM–FM decomposition methods.
6 Summary
In the final part of this paper, we have reported an endtoend, solution for the estimation of instantaneous parameters of the AM–FM model assuming IMF components. By leveraging the theory developed in Part II of these papers to interpret and demodulate the results returned by EMD, we have provided a complete numerical method for HSA. We provided examples of HSA–IMF on synthetic signals and argued that the resulting decompositions were more representative of the underlying signal models as compared to conventional Fourier analysis. Examples of HSA–IMF on realworld signals were shown to allow for alternative and possibly more useful interpretation of the underlying signal model. Finally, we discussed computational aspects of the proposed algorithm.
\pdfstringdef\BKM@titleConclusions
Conclusions
In this paper, we began by reframing the classic complex extension problem as a Latent Signal Analysis (LSA) problem where the objective is to determine the complexvalued latent signal, z(t) from the realvalued observation, x(t)=\Re\{z(t)\}. We used this framework to argue against the assumption of Harmonic Correspondence (HC), and hence the use of Gabor’s quadrature method and the Hilbert Transform (HT) for complex extension. By relaxing the HC assumption, there are many choices for the quadrature and furthermore, many of these complex extensions can be useful in modeling physical phenomena.
Next, we presented a theory for the Hilbert spectrum framed as generalized LSA problem in which we seek a representation of z(t) consisting of a superposition of latent AM–FM components parameterized by a set of Instantaneous Amplitude (IA)/Instantaneous Frequency (IF) pairs. The use of latent AM–FM components admits many possible forms of the component therefore for a given x(t), there is considerable freedom in the signal model. We presented the analogue of the LSA problem in the frequency domain where we showed that a latent spectrum cannot be uniquely tied to a to the spectrum of the real observation because of the structure imposed by the real operator. We showed that without the HC assumption, there are many choices for the latent spectrum and these spectra will not have in Hermitian symmetry. We proposed a novel 3D visualization of the Hilbert spectrum which plots \omega(t) vs. s(t) vs. t and coloring with respect to a(t) and allows for the visualization of the instantaneous parameters. We have recast timefrequency analysis and Gabor’s analogies to quantum mechanics to an analysis method where uncertainty is in the quadrature signal and not in frequency. Further, by moving away from Simple Harmonic Components (SHCs) with HC, we allow for a new and powerful way to analyze nonstationary signals.
We recognized that an Intrinsic Mode Function (IMF) is a latent AM–FM component and leveraged HSA theory to interpret both the IMF and Empirical Mode Decomposition (EMD). With this interpretation we show that the definition of an IMF unambiguously forces a unique complex extension. Furthermore, we also recognize fundamental problems with the HilbertHuang Transform (HHT), i.e. that the HT is inappropriate for IMF demodulation and proposed an IMF demodulation method that is compatible the the IMF definition. Finally, we utilized the EMD algorithm with our modifications and IMF demodulation to calculate the IA/IF parameters of x(t), thus providing a numerical method for Hilbert Spectral Analysis (HSA).
Acknowledgments
The authors wish to thank Prof. Joe Lakey of the Dept. of Mathematical Sciences at New Mexico State University and Prof. Antonia PapandreouSuppappola of the School of Electrical, Computer and Energy Engineering at Arizona State University for reviewing an early draft of this manuscript.
References
 (1) L. Cohen, TimeFrequency Analysis, Prentice Hall, 1995.
 (2) R. Shankar, Fundamentals of Physics: Mechanics, Relativity and Thermodynamics, Yale University Press, 2014.
 (3) L. Kinsler, A. Frey, A. Coppens, J. Sanders, Fundamentals of Acoustics, 3rd Edition, Wiley Publishing, 1982.
 (4) J. R. Carson, T. C. Fry, Variable frequency electric circuit theory, Bell System Technical Journal 16 (1937) 513–540.
 (5) B. V. D. Pol, The fundamental principles of frequency modulation, J. Inst. Electr. Eng. 3 93 (23) (1946) 153–158.
 (6) D. Gabor, Theory of communication, part 1: the analysis of information, J. Inst. Electr. Eng. 3 93 (26) (1946) 429–441.
 (7) J. Shekel, ‘Instantaneous’ frequency, Proc. IRE 41 (4) (1953) 548–548.
 (8) B. Picinbono, On instantaneous amplitude and phase of signals, IEEE Trans. Signal Process. 45 (3) (1997) 552–560.
 (9) B. Boashash, Estimating and interpreting the instantaneous frequency of a signal–part 1: fundamentals, Proc. IEEE 80 (4) (1992) 520 –538.
 (10) J. Hurpert, Instantaneous frequency, Proc. IRE 41 (1953) 1188–1188.
 (11) A. W. Rihaczek, E. Bedrosian, Hilbert transforms and the complex representation of real signals, Proc. IEEE 54 (3) (1966) 434–435.
 (12) D. E. Vakman, On the definition of concepts of amplitude phase and instantaneous frequency, Radio Eng. and Electron. Phys. 17 (1972) 754–759.
 (13) D. E. Vakman, Do we know what are the instantaneous frequency and instantaneous amplitude of a signal, Radio Eng. and Electron. Phys. 21 (1976) 95–100.
 (14) D. E. Vakman, Measuring the frequency of an analytical signal, Radio Eng. and Electron. Phys. 24 (1979) 63–69.
 (15) P. J. Loughlin, B. Tacer, On the amplitude and frequencymodulation decomposition of signals, J. Acoust. Soc. Am. 100 (3) (1996) 1594–1601.
 (16) P. J. Loughlin, Do bounded signals have bounded amplitudes?, Multidim. Syst. Sig. Proc. 9 (4) (1998) 419–424.
 (17) D. Vakman, On the analytic signal, the teagerkaiser energy algorithm and other methods for defining amplitude and frequency, IEEE Trans. Signal Process. 44 (4) (1996) 791–797.
 (18) D. Vakman, Signals, Oscillations and Waves, Artech House, 1998.
 (19) E. C. Titchmarch, Theory of Fourier Integrals, 2nd Edition, Oxford University Press, 1948.
 (20) C. W. Therrien, The LeeWiener legacy [statistical theory of communication], IEEE Signal Process. Mag. 19 (6) (2002) 33–34.
 (21) D. E. Vakman, L. A. Vainshtein, Amplitude, phase, frequency–fundamental concepts in the theory of oscillations, Uspekhi Fizicheskikh Nauk 123 (1977) 657–682.
 (22) L. Cohen, P. Loughlin, D. Vakman, On an ambiguity in the definition of the amplitude and phase of a signal, Signal Process. 79 (3) (1999) 301–307.
 (23) B. Boashash, Time Frequency Signal Analysis and Processing, Elsevier, 2003.
 (24) J. Ville, Theorie et applications de la notion de signal analytique, Cables et Transmission 2a (1948) 61–74.
 (25) E. Bedrosian, The analytic signal representation of modulated waveforms, Proc. IRE 50 (10) (1962) 2071–2076.
 (26) X. G. Xia, L. Cohen, On analytic signals with nonnegative instantaneous frequency, in: Proc. IEEE Int. Conf. Acoust. Speech Signal Process., 1999.
 (27) J. W. Brown, R. V. Churchill, Complex Variables and Applications, McGrawHill Higher Education, 2009.
 (28) B. V. der Pol, The nonlinear theory of electric oscillations, Proc. IRE 22 (9) (1934) 1051–1086.
 (29) S. Farlow, Partial differential equations for scientists and engineers, Courier Corporation, 2012.
 (30) E. Bedrosian, A product theorem for Hilbert transforms, Proc. IEEE 51 (5) (1963) 868–869.
 (31) A. H. Nuttall, E. Bedrosian, On the quadrature approximation to the Hilbert transform of modulated signals, Proc. IEEE 54 (10) (1966) 1458–1459.
 (32) J. R. Carson, Notes on the theory of modulation, Proc. IRE 10 (1) (1922) 57–64.
 (33) M. B. Priestley, Nonlinear and Nonstationary Time Series Analysis, Academic Press, 1988.
 (34) L. M. Fink, Relations between the spectrum and instantaneous frequency of a signal, Problemy Peredachi Informatsii 2 (4) (1966) 26–38.
 (35) L. Mandel, Interpretation of instantaneous frequencies, Am. J. of Phys. 42 (10) (1974) 840–846.
 (36) M. S. Gupta, Definition of instantaneous frequency and frequency measurability, Am. J. of Phys. 43 (12) (1975) 1087–1088.
 (37) P. J. Loughlin, B. Tacer, Comments on the interpretation of instantaneous frequency, IEEE Signal Process. Lett. 4 (5) (1997) 123–125.
 (38) P. M. Oliveira, V. Barroso, On the concept of instantaneous frequency, in: Proc. IEEE Int. Conf. Acoust. Speech Signal Process., Vol. 4, IEEE, 1998, pp. 2241–2244.
 (39) P. M. Oliveira, V. Barroso, Instantaneous frequency of multicomponent signals, IEEE Signal Process. Let. 6 (4) (1999) 81–83.
 (40) W. Nho, P. J. Loughlin, When is instantaneous frequency the average frequency at each time?, IEEE Signal Process. Lett. 6 (4) (1999) 78–80.
 (41) C. K. Chui, H. K. Mhaskar, Signal decomposition and analysis via extraction of frequencies, Appl. Comput. Harmonic Anal.
 (42) N. E. Huang, Z. Shen, S. R. Long, M. C. Wu, H. H. Shih, Q. Zheng, N. C. Yen, C. C. Tung, H. H. Liu, The empirical mode decomposition and the Hilbert spectrum for nonlinear and nonstationary time series analysis, Proc. R. Soc. London Ser. A 454 (1971) (1998) 903–995.
 (43) B. Boashash, TimeFrequency Signal Analysis, Longman Publishing Group, 1992.
 (44) L. Cohen, What is a multicomponent signal?, in: Proc. IEEE Int. Conf. Acoust. Speech Signal Process., Vol. 5, 1992, pp. 113–116.
 (45) D. Wei, A. C. Bovik, On the instantaneous frequencies of multicomponent AM–FM signals, IEEE Signal Process. Lett. 5 (4) (1998) 84–86.
 (46) L. Cohen, C. Lee, Instantaneous frequency, its standard deviation and multicomponent signals, in: Proc. SPIE Int. Soc. Opt. Eng., 1988, pp. 186–208.
 (47) J. Fourier, On the propagation of heat in solid bodies.
 (48) J. Fourier, The analytical theory of heat, Firmin Didot, 1822.
 (49) A. V. Oppenheim, A. Willsky, S. H. Nawab, Signals and Systems, 2nd Edition, Prentice Hall, 1997.
 (50) J. B. Allen, L. Rabiner, A unified approach to shorttime fourier analysis and synthesis, Proc. IEEE 65 (11) (1977) 1558–1564.
 (51) L. Cohen, Timefrequency distributions–a review, Proc. IEEE 77 (7) (1989) 941–981.
 (52) G. Jones, B. Boashash, Instantaneous frequency, instantaneous bandwidth and the analysis of multicomponent signals, in: Proc. IEEE Int. Conf. Acoust. Speech Signal Process., 1990, pp. 2467–2470.
 (53) L. Cohen, Instantaneous ‘anything’, in: Proc. IEEE Int. Conf. Acoust. Speech Signal Process., 1993, pp. 105–108.
 (54) B. C. Lovell, R. C. Williamson, B. Boashash, The relationship between instantaneous frequency and timefrequency representations, IEEE Trans. Signal Process. 41 (3) (1993) 1458–1461.
 (55) S. Qian, C. Dapang, Joint timefrequency analysis, IEEE Signal Process. Mag. 16 (2) (1999) 52–67.
 (56) G. Strang, Introduction to Linear Algebra, WellesleyCambridge Press, 2009.
 (57) E. H. Armstrong, Some recent developments in the audion receiver, Proc. IRE 3 (3) (1915) 215–238.
 (58) H. B. Voelcker, Toward a unified theory of modulation part I: Phaseenvelope relationships, Proc. IEEE 54 (3) (1966) 340–353.
 (59) H. B. Voelcker, Toward a unified theory of modulationï¿½ part II: zero manipulation, Proc. IEEE 54 (5) (1966) 735–755.
 (60) J. M. Chowning, Origins of FM Synthesis, http://www.youtube.com/watch?v=w4g92vX1YF4 (2012).
 (61) R. Bedaux, Micro frequency modulation in sound synthesis, J. New Music Res. 3 (2) (1974) 89–108.
 (62) J. M. Chowning, The synthesis of complex audio spectra by means of frequency modulation, J. Audio Eng. Soc. 21 (7) (1977) 526–534.
 (63) M. Feldman, Nonlinear system vibration analysis using Hilbert transform–I. free vibration analysis method ’FREEVIB’, Mech. Syst. and Signal Process. 8 (2) (1994) 119–127.
 (64) A. Rao, R. Kumaresan, On decomposing speech into modulated components, IEEE Trans. Speech Audio Process. 8 (3) (2000) 240–254.
 (65) F. Gianfelici, G. Biagetti, P. Crippa, C. Turchetti, Multicomponent AM–FM representations: an asymptotically exact approach, IEEE Trans. Audio Speech Lang. Process. 15 (3) (2007) 823–837.
 (66) M. Feldman, Hilbert Transform Applications in Mechanical Vibration, Wiley, 2011.
 (67) R. McAulay, T. F. Quatieri, Speech analysis/synthesis based on a sinusoidal representation, IEEE Trans. Acoust., Speech, Signal Process. 34 (4) (1986) 744–754.
 (68) P. Rao, F. J. Taylor, Estimation of instantaneous frequency using the discrete Wigner distribution, Electron. Lett. 26 (4) (1990) 246–248.
 (69) Y. Pantazis, O. Rosec, Y. Stylianou, Adaptive AM–FM Signal Decomposition With Application to Speech Analysis, IEEE Trans. Audio Speech Lang. Process. 19 (2) (2011) 290–300.
 (70) B. Boashash, G. Azemi, J. O’Toole, Timefrequency processing of nonstationary signals: Advanced tfd design to aid diagnosis with highlights from medical applications, IEEE Signal Process. Mag. 30 (6) (2013) 108–119.
 (71) P. Maragos, J. F. Kaiser, T. F. Quatieri, Energy separation in signal modulations with application to speech analysis, IEEE Trans. Signal Process. 41 (10) (1993) 3024–3051.
 (72) A. C. Bovik, P. Maragos, T. F. Quatieri, AM–FM energy detection and separation in noise using multiband energy operators, IEEE Trans. Signal Process. 41 (12) (1993) 3245–3265.
 (73) L. B. Fertig, J. H. McClellan, Instantaneous frequency estimation using linear prediction with comparisons to the DESAs, IEEE Signal Process. Lett. 3 (2) (1996) 54–56.
 (74) A. Potamianos, P. Maragos, Speech analysis and synthesis using an AM–FM modulation model, Speech Commun. 28 (3) (1999) 195–209.
 (75) A. O. Boudraa, J. C. Cexus, F. Salzenstein, L. Guillon, If estimation using empirical mode decomposition and nonlinear teager energy operator, in: Int. Symp. Cont., Comm. Signal Process., 2004, pp. 45–48.
 (76) A. O. Boudraa, Instantaneous frequency estimation of fm signals by \psibenergy operator, Electron. Lett. 47 (10) (2011) 623–624.
 (77) T. F. Quatieri, T. E. Hanna, G. C. O’Leary, AM–FM separation using auditorymotivated filters, IEEE Trans. Speech Audio Process. 5 (5) (1997) 465–480.
 (78) F. Gianfelici, C. Turchetti, P. Crippa, Multicomponent AM–FM demodulation: the state of the art after the development of the iterated Hilbert transform, in: Proc. Int. Conf. Sig. Process. and Commun., 2007, pp. 1471–1474.
 (79) B. Boashash, P. O’Shea, M. J. Arnold, Algorithms for instantaneous frequency estimation: a comparative study, in: Proc. SPIE Int. Soc. Opt. Eng., 1990, pp. 126–148.
 (80) B. Boashash, Estimating and interpreting the instantaneous frequency of a signal–part II: algorithms and applications, Proc. IEEE 80 (4) (1992) 540–568.
 (81) M. Feldman, Theoretical analysis and comparison of the Hilbert transform decomposition methods, Mech. Syst. and Signal Process. 22 (3) (2008) 509–519.
 (82) T. F. Quatieri, DiscreteTime Speech Signal Processing, Prentice Hall, 2002.
 (83) A. Potamianos, P. Maragos, A comparison of the energy operator and the hilbert transform approach to signal and speech demodulation, Signal Process. 37 (1) (1994) 95–120.
 (84) E. S. Diop, A. O. Boudraa, F. Salzenstein, A joint 2D AM–FM estimation based on higher order Teager–Kaiser energy operators, Signal, Image and Video Process. 5 (1) (2011) 61–68.
 (85) B. Santhanam, P. Maragos, Multicomponent AM–FM demodulation via periodicitybased algebraic separation and energybased demodulation, IEEE Trans. Commun. 48 (3) (2000) 473–490.
 (86) J. A. Moorer, The synthesis of complex audio spectra by means of discrete summation formulas, J. Audio Eng. Soc. 24 (9) (1976) 717–727.
 (87) C. Dodge, T. A. Jerse, Computer Music: Synthesis, Composition and Performance, Macmillan Library Reference, 1997.
 (88) P. J. Loughlin, L. Cohen, The uncertainty principle: global, local, or both?, IEEE Trans. Signal Process. 52 (5) (2004) 1218–1227.
 (89) Hilbert Spectral Analysis, http://www.HilbertSpectrum.com (2015).
 (90) D. Borland, R. M. Taylor, Rainbow color map (still) considered harmful., IEEE Trans. Visual. Comput. Graphics 27 (2) (2007) 14–17.
 (91) M. Niccoli, S. Lynch, A more perceptual color palette for structure maps, in: Proc. GeoConvention, 2012.
 (92) Z. Wu, N. E. Huang, Ensemble empirical mode decomposition: a noiseassisted data analysis method, Adv. Adapt. Data Anal. 1 (01) (2009) 1–41.
 (93) R. Deering, J. F. Kaiser, The use of a masking signal to improve empirical mode decomposition, in: Proc. IEEE Int. Conf. Acoust. Speech Signal Process., 2005, pp. 485–488.
 (94) M. E. Torres, M. A. Colominas, G. Schlotthauer, P. Flandrin, A complete ensemble empirical mode decomposition with adaptive noise, in: Proc. IEEE Int. Conf. Acoust. Speech Signal Process., 2011, pp. 4144–4147.
 (95) R. Rato, M. D. Ortigueira, A. Batista, On the HHT, its problems and some solutions, Mech. Syst. and Signal Process. 22 (6) (2008) 1374–1394.
 (96) N. E. Huang, Z. Wu, S. R. Long, K. C. Arnold, X. Chen, K. Blank, On instantaneous frequency, Adv. Adapt. Data Anal. 1 (02) (2009) 177–229.
 (97) Z. Wu, N. E. Huang, A study of the characteristics of white noise using the empirical mode decomposition method, Proc. R. Soc. London Ser. A 460 (2046) (2004) 1597–1611.
 (98) P. Flandrin, G. Rilling, P. Goncalves, Empirical mode decomposition as a filter bank, IEEE Signal Process. Lett. 11 (2) (2004) 112–114.
 (99) P. Flandrin, P. Gonçalves, G. Rilling, EMD Equivalent Filter Banks, from Interpretation to Applications, in: N. Huang, S. Shen (Eds.), The HilbertHuang Transform and Its Applications, Interdisciplinary Math. Sciences, 2005, pp. 55–74.
 (100) D. Mandic, et al., Filter bank property of multivariate empirical mode decomposition, IEEE Trans. Signal Process. 59 (5) (2011) 2421–2426.
 (101) S. Meignen, V. Perrier, A new formulation for empirical mode decomposition based on constrained optimization, IEEE Signal Process. Lett. 14 (12) (2007) 932–935.
 (102) Y. Kopsinis, S. McLaughlin, Investigation and performance enhancement of the empirical mode decomposition method based on a heuristic search optimization approach, IEEE Trans. Signal Process. 56 (1) (2008) 1–13.
 (103) T. Y. Hou, Z. Shi, Datadriven time–frequency analysis, Appl. Comput. Harmonic Anal. 35 (2) (2013) 284–308.
 (104) a. D. P. M. D. Looney, A machine learning enhanced empirical mode decomposition, in: Proc. IEEE Int. Conf. Acoust. Speech Signal Process., 2008, pp. 1897–1900.
 (105) E. Deléchelle, J. Lemoine, O. Niang, Empirical mode decomposition: an analytical approach for sifting process, IEEE Signal Process. Lett. 12 (11) (2005) 764–767.
 (106) R. Sharpley, V. Vatchev, Analysis of the intrinsic mode functions, Constr. Approx. 24 (1) (2006) 17–47.
 (107) V. Vatchev, R. Sharpley, Decomposition of functions into pairs of intrinsic mode functions, Proc. R. Soc. Lond. A Math. Phys. Sci. 464 (2097) (2008) 2265–2280.
 (108) E. H. Diop, R. Alexandre, A. O. Boudraa, A pde characterization of the intrinsic mode functions, in: Proc. IEEE Int. Conf. Acoust. Speech Signal Process., 2009, pp. 3429–3432.
 (109) E. S. Diop, R. Alexandre, A. O. Boudraa, Analysis of intrinsic mode functions: a pde approach, IEEE Signal Process. Lett. 17 (4) (2010) 398–401.
 (110) S. D. E. Hadji, R. Alexandre, V. Perrier, A pde based and interpolationfree framework for modeling the sifting process in a continuous domain, Adv. Comput. Math. 38 (4) (2013) 801–835.
 (111) O. Niang, E. Deléchelle, J. Lemoine, A spectral approach for sifting process in empirical mode decomposition, IEEE Trans. Signal Process. 58 (11) (2010) 5612–5623.
 (112) C. Damerval, S. Meignen, V. Perrier, A fast algorithm for bidimensional emd, IEEE Signal Process. Lett. 12 (10) (2005) 701–704.
 (113) G. Rilling, P. Flandrin, P. Gonçalves, J. M. Lilly, Bivariate empirical mode decomposition, IEEE Signal Process. Lett. 14 (12) (2007) 936–939.
 (114) Z. Wu, N. E. Huang, X. Chen, The multidimensional ensemble empirical mode decomposition method, Adv. Adapt. Data Anal. 1 (03) (2009) 339–372.
 (115) A. Linderhed, Image empirical mode decomposition: A new tool for image processing, Adv. Adapt. Data Anal. 1 (02) (2009) 265–294.
 (116) J. C. Nunes, E. Deléchelle, Empirical mode decomposition: Applications on signal and image processing, Adv. Adapt. Data Anal. 1 (01) (2009) 125–175.
 (117) D. P. Mandic, N. U. Rehman, W. Zhaohua, N. E. Huang, Empirical mode decompositionbased timefrequency analysis of multivariate signals: The power of adaptive data analysis, IEEE Signal Process. Mag. 30 (6) (2013) 74–86.
 (118) T. Tanaka, D. P. Mandic, Complex empirical mode decomposition, IEEE Signal Process. Lett. 14 (2) (2007) 101–104.
 (119) M. Feldman, Timevarying vibration decomposition and analysis based on the Hilbert transform, J. Sound Vib. 295 (3) (2006) 518–530.
 (120) X. Chen, Z. Wu, N. E. Huang, The timedependent intrinsic correlation based on the empirical mode decomposition, Adv. Adapt. Data Anal. 2 (02) (2010) 233–265.
 (121) D. I, J. Lu, H. T. Wu, Synchrosqueezed wavelet transforms: an empirical mode decompositionlike tool, Appl. Comput. Harmonic Anal. 30 (2) (2011) 243–261.
 (122) H. T. Wu, P. Flandrin, I. Daubechies, One or two frequencies? the synchrosqueezing answers, Adv. Adapt. Data Anal. 3 (01n02) (2011) 29–39.
 (123) T. Y. Hou, Z. Shi, Adaptive data analysis via sparse timefrequency representation, Adv. Adapt. Data Anal. 3 (01n02) (2011) 1–28.
 (124) G. Rilling, P. Flandrin, On the influence of sampling on the empirical mode decomposition, in: Proc. IEEE Int. Conf. Acoust. Speech Signal Process., 2006, pp. 444–447.
 (125) R. Rato, M. Ortigueira, A Modified EMD Algorithm for Application in Biomedical Signal Processing, in: Int. Conf. Comput. Intell. Med. Healthcare, 2005.
 (126) N. E. Huang, M. L. C. Wu, S. R. Long, S. S. P. Shen, W. Qu, P. Gloersen, K. L. Fan, A confidence limit for the empirical mode decomposition and hilbert spectral analysis, Proc. R. Soc. London Ser. A 459 (2037) (2003) 2317–2345.
 (127) G. Rilling, P. Flandrin, P. Goncalves, et al., On empirical mode decomposition and its algorithms, in: IEEEEURASIP Workshop on Nonlinear Signal and Image Process., Vol. 3, 2003, pp. 8–11.
 (128) G. Rilling, P. Flandrin, One or two frequencies? The empirical mode decomposition answers, IEEE Trans. Signal Process. 56 (1) (2008) 85–95.
 (129) N. Senroy, S. Suryanarayanan, Two techniques to enhance empirical mode decomposition for power quality applications, in: Proc. IEEE Power Eng. Soc. General Meeting, 2007, pp. 1–6.
 (130) X. Guanle, W. Xiaotong, X. Xiaogang, Timevarying frequencyshifting signalassisted empirical mode decomposition method for AM–FM signals, Mech. Syst. and Signal Process. 23 (8) (2009) 2458–2469.
 (131) X. M. Li, C. C. Bao, M. S. Jia, A sinusoidal audio and speech analysis/synthesis model based on improved EMD by adding pure tone, in: IEEE Int. Workshop Mach. Learn. Signal Process., 2011, pp. 1–5.
 (132) S. D. Hawley, L. E. Atlas, H. J. Chizeck, Some properties of an empirical mode type signal decomposition algorithm, IEEE Signal Process. Lett. 17 (1) (2010) 24–27.
 (133) T. Strom, On amplitudeweighted instantaneous frequencies, IEEE Trans. Acoust., Speech, Signal Process. 25 (4) (1977) 351–353.
 (134) M. Feldman, Analytical basics of the emd: Two harmonics decomposition, Mech. Syst. and Signal Process. 23 (7) (2009) 2059–2071.
 (135) D. O’Shaughnessy, Speech communication: human and machine, AddisonWesley Pub. Co., 1987.
 (136) P. Waskito, S. Miwa, Y. Mitsukura, H. Nakajo, Parallelizing HilbertHuang transform on a GPU, in: Networking and Computing (ICNC), 2010 First International Conference on, 2010, pp. 184–190.
 (137) D. Chen, D. Li, M. Xiong, H. Bao, X. Li, GPGPUaided ensemble empiricalmode decomposition for EEG analysis during anesthesia, IEEE Trans. Inf. Technol. Biomed. 14 (6) (2010) 1417–1427.
 (138) L. W. Chang, M. T. Lo, N. Anssari, K. H. Hsu, N. E. Huang, W. M. Hwu, Parallel implementation of multidimensional ensemble empirical mode decomposition, in: Proc. IEEE Int. Conf. Acoust. Speech Signal Process., 2011, pp. 1621–1624.
 (139) P. Waskito, Y. Mitsukura, H. Nakajo, Evaluation of GPUBased Empirical Mode Decomposition for OffLine Analysis, IEICE Trans. on Inf. and Syst. 94 (12) (2011) 2328–2337.
 (140) N. Senroy, S. Suryanarayanan, P. F. Ribeiro, An improved Hilbert–Huang method for analysis of timevarying waveforms in power quality, IEEE Trans. Power Syst. 22 (4) (2007) 1843–1850.
 (141) S. Riemenschneider, B. Liu, Y. Xu, N. E. Huang, Bspline based empirical mode decomposition, Interdisciplinary Math. Sciences 5 (2005) 27–56.
 (142) S. Qin, Y. M. Zhong, A new envelope algorithm of Hilbert–Huang transform, Mech. Syst. and Signal Process. 20 (8) (2006) 1941–1952.
 (143) Q. Chen, N. Huang, S. Riemenschneider, Y. Xu, A bspline approach for empirical mode decompositions, Adv. Comput. Math. 24 (14) (2006) 171–195.
 (144) Y. Kopsinis, S. McLaughlin, Improved EMD using doublyiterative sifting and high order spline interpolation, J. Adv. Signal Process. 2008 (2008) 120.
 (145) F. Wang, X. Y. CHEN, F. L. Qiao, Z. Wu, N. E. Huang, On intrinsic mode function, Adv. Adapt. Data Anal. 2 (03) (2010) 277–293.
 (146) A. Bouchikhi, A. O. Boudraa, Multicomponent am–fm signals analysis based on emd–bsplines esa, Signal Process. 92 (9) (2012) 2214–2228.