A Real time Alternative to the Hilbert Huang Transform Based on Internal Model Principle

and n(t) is measurement noise. These are signals that are the sums of n periodic components with each component composed of mi harmonics. The periods, the harmonic amplitudes and relative phases can vary slowly in time. By identification, we mean determining the values ωi, Āij and φij φ11. Several techniques have been developed in the literature to solve this problem. The most traditional technique is the fast Fourier transform. Newer techniques include wavelet analysis. These approaches suffer from not allowing continuous estimations of the frequencies and have difficult trade-offs between time and frequency resolutions. Other approaches are based on the use of adaptive notch filters [1] and output regulation [2]. A new approach that has been widely applied is the Hilbert Huang Transform (HHT) [3]. Control engineers treat similar problems where exact tracking of reference signals or rejection of disturbances is required. Approaches that accomplish this include repetitive controllers [4] and adaptive feed-forward cancellation (AFC) [5]. The repetitive controller is based on a fundamental control theory principle called the internal model principle (IMP). This principle was presented by Francis and Wonham and states that the output error can be driven asymptotically to zero by placing a model of exogenous signals in a stable feedback loop [6]. Unfortunately small errors in this model can lead to significant degradation in the performance of internal model principle controllers. This problem of uncertainty in the signal model can be overcome with adaptive controllers [7]. In achieving asymptotically perfect rejection of disturbances it is inherent that the disturbance is completely identified. Thus, these types of controllers can be turned into signal processing algorithms by replacing the process to be controlled with tuning functions [8]. Unfortunately, to successfully implement this algorithm requires being able to tune a stable feedback control loop for the entire range of possible frequencies in the model given by equation (1). Fortunately, it has been shown that in the signal processing framework, the simplest tuning solution, i.e. selecting all of the gains to be one, is guaranteed to be stable. This algorithm has been successfully applied to the problem of the repeatable disturbances seen in disk drive head control [9]. Unfortunately, by resorting to this simple tuning approach, there is no control over the dynamics and noise rejection characteristics of the algorithm.


Introduction
In this article, we are interested in the problem of identifying signals of the following form  (2) and n(t) is measurement noise. These are signals that are the sums of n periodic components with each component composed of m i harmonics. The periods, the harmonic amplitudes and relative phases can vary slowly in time. By identification, we mean determining the values ω i , Ā ij and φ ij -φ 11 .
Several techniques have been developed in the literature to solve this problem. The most traditional technique is the fast Fourier transform. Newer techniques include wavelet analysis. These approaches suffer from not allowing continuous estimations of the frequencies and have difficult trade-offs between time and frequency resolutions. Other approaches are based on the use of adaptive notch filters [1] and output regulation [2]. A new approach that has been widely applied is the Hilbert Huang Transform (HHT) [3]. Control engineers treat similar problems where exact tracking of reference signals or rejection of disturbances is required. Approaches that accomplish this include repetitive controllers [4] and adaptive feed-forward cancellation (AFC) [5]. The repetitive controller is based on a fundamental control theory principle called the internal model principle (IMP). This principle was presented by Francis and Wonham and states that the output error can be driven asymptotically to zero by placing a model of exogenous signals in a stable feedback loop [6]. Unfortunately small errors in this model can lead to significant degradation in the performance of internal model principle controllers. This problem of uncertainty in the signal model can be overcome with adaptive controllers [7]. In achieving asymptotically perfect rejection of disturbances it is inherent that the disturbance is completely identified. Thus, these types of controllers can be turned into signal processing algorithms by replacing the process to be controlled with tuning functions [8].
Unfortunately, to successfully implement this algorithm requires being able to tune a stable feedback control loop for the entire range of possible frequencies in the model given by equation (1). Fortunately, it has been shown that in the signal processing framework, the simplest tuning solution, i.e. selecting all of the gains to be one, is guaranteed to be stable. This algorithm has been successfully applied to the problem of the repeatable disturbances seen in disk drive head control [9]. Unfortunately, by resorting to this simple tuning approach, there is no control over the dynamics and noise rejection characteristics of the algorithm.
When the frequencies are known a priori, the report [10] shows how the dynamics of the algorithm can be completely specified. Unfortunately this article requires solving a set more than 2 2 t i n m = ∑ coupled linear equations which are a function of the signal's frequencies.
Unless the sample rate is less than 1Hz this will not be feasible to do each sample. This article shows how these parameters can be explicitly solved by simply evaluating some frequency response functions at certain frequencies.
In Section II, an instantaneous Fourier decomposition (IFD) algorithm [11] that is similar in approach to the HHT is presented. In Section III an updated formula for calculating the instantaneous frequencies are given. In Section IV, the new realtime tuned algorithm is presented. In Section V, the ability of the proposed algorithm to identify the periodic signal with uncertain frequencies is demonstrated. Conclusions are drawn in Section VI.
A preliminary version of this article was presented at the 30th annual IEEE Canadian Conference on Electrical and Computer Engineering (IEEE 2017 CCECE) in Windsor [12].

Abstract
This article presents a new tuning approach for an adaptive internal-model-principle based signal identification algorithm whose computational costs are low enough to allow a realtime implementation. The algorithm allows an instantaneous Fourier decomposition of non-stationary signals that have a strongly predictable component. The algorithm is implemented as a feedback loop resulting in a closed loop system with a frequency response of a bandpass filter with notches at the frequencies of the Fourier decomposition. This is achieved through real time selection of the coefficients of the transfer functions in the feedback loop. Previously these coefficients were selected by solving a large set of coupled linear equations. Rules for explicitly solving for these parameters are given that only involve evaluating frequency responses at the frequencies of the instantaneous Fourier decomposition. This allows realtime implementation on a low cost lap top with sampling rates up to 10 kHz.

Adaptive Algorithm and Comparison to HHT
The HHT proceeds from the realization that the Hilbert transform gives a mathematically precise definition of instantaneous frequency that agrees with our intuitive understanding when applied to narrowband signals. In this narrowband case, the instantaneous frequency can be approximated as the derivative of the angle of the narrowband signal and 1 − times the quadrature of that signal where the quadrature can be approximated by either a scaled version of the derivative or integral of the signal. The HHT uses an empirical method to break down signals into narrowband signals. This empirical method is numerically intensive and not compatible with a realtime implementation.
Our algorithm uses the same approximations to estimate the instantaneous frequencies as the HHT but uses an alternative, notch filter based approach that simultaneous calculates the quadrature signals and decomposes the signal into narrow band signals. The structure of the adaptive instantaneous frequency decomposition is shown in Figure 1, where G(s) is a tuning function.
Each of the transfer functions IM i,j are an internal model for a sinusoid of frequency i j ω *  . When the model frequencies and the signal frequencies match, i.e., i i ω ω =  and the closed loop system is stable, each u ij will be a single sinusoidal and meet the HHT definition of an intrinsic function. The basic algorithm is the state space based implementation of the internal models given by This is taken from [11] with minor modifications to fit the signal model that was given in equations (1) and (2). The gains K 1ij , K 2ij have been moved to the input vector from the output vector so that adjustments in their value do not directly change u, i.e., a bumpless transfer. Consequently, the responses at x 1ij (t) and x 2ij (t) in steady state are: i.e., the second state is the sinusoidal component of the original signal and the first state is its quadrature. While the states are time varying, when the signal parameters are time invariant Since the state variables x 1i1 and x 2i1 are orthogonal to each other then, as with the HHT, the derivative of the angle of 1i1 Thus, using an integral controller can be used to update the frequency estimates.
Thus a quasi-periodic signal can be decomposed into a sum of narrow band signals, {u ij }={x 2ij }, and a real time Fourier representation of the reference can be obtained. The signal u(t) is the estimate of the signal of interest and can be represented by In ref. [11], it is establish for sufficiently small K ai the algorithm is locally exponentially stable when G(s) and the K 1ij , K 2ij are chosen so that the feedback loop in Figure 1 is stable at each point in time. Designing these controller parameters is a challenging problem as it is assumed that there is limited knowledge about the {ω i } and during transients there can be a significant difference between {ω i } and { } i ω  .

Off-line tuning
As with ref. [10], we satisfy the above stability assumption by designing the closed loop system to incorporate a bandpass filter with notch filter. Let a 2 nd order desirable bandpass filter be given by We choose the controller parameters to be such that the transfer function from d to e is where ij ε are small real numbers, and i jω  are the notches frequency. The presence of the numerator of the second term is a fundamental consequence of the internal model principle. Therefore, the ability of the algorithm to improve noise rejection is achieved.
An analysis of Figure 1 gives Where, ( ) The contribution of this article is to develop a less computationally intensive algorithm for calculating the controller parameter to meet the realtime requirement.

On-line frequency identification
Now the crucial question is how to choose G(s) and K 1jk , K 2jk and implement the algorithm without needing to solve a set or 2n t +4 linear equations. It can be seen that all of the terms in the denominator except the term containing Y kl will be zero if This generates 2 complex and complementary conjugate equations with 2 unknowns, i.e., the real part of either equation gives K 1jk and the imaginary gives K 2jk . The 4 a i parameters can be explicitly solved by equating the coefficients of the degree 0, 1, 2n t +2, 2n t +3 terms of the denominator. Note the second term of the denominator of equation (12) contribute nothing to these four terms. These coefficients can be calculated by utilizing the relationships between the coefficients of a polynomial and the roots of a polynomial. We have that  (12) will be zero when we substitute in = 1 k s l ω ± − and it will not be possible to calculate two pairs of internal model gains. Further, while it is theoretical possible to solve when the frequencies are extremely close, we get solutions that lead to unstable results because of numerical stability issues. To solve this problem, while calculating the controller gains, we drop the approximately redundant internal model when the frequencies become close, i.e., within 0.1%. After calculating the controller gains, the two redundant models are each assigned half of the gain. That is when jω i =lω k , we drop Internal model IM l,k from the design stage. Let 1ij K and 2ij K be the calculated controllers gains. Then 1ij 1kl 1ij K = K = 0.5K and 2ij 2kl 2ij K = K = 0.5K . It should be noted that the threshold for HHT to distinguish between close frequencies is 10%.

Simulation Results
In this particular section, the effectiveness of our real time implementation of our proposed adaptive algorithm is verified via simulation. The model configuration parameters that are used with the matlab/simulink (R2016) environment are as follows: Solver ode5 (Dormand-prince) selection with fundamental sample time is 0.0025 s. Therefore, the sampling rate in our case is selected to be 400Hz, then the Nyquist frequency is 200 Hz. The code generation with C language and tool chain (Microsoft visual C++ 2012 V11.1 n-make 164-bit windows). All random numbers were zero mean.
Our signal to be identified was produced by summing the outputs of two copies of the model shown in Figure 2. The feedback loop containing the pure delay is called a repetitive controller and is capable of producing any periodic disturbance with period T. The value T was an integrated band limited white disturbance. The frequency cutoff of this noise was 20 rad/s and the variance was 0.5. The initial conditions for both fundamental frequencies are 4.2 and 5 Hz. The disturbance input to the repetitive controller causes the amplitudes and relative phases to vary slowly with time as well. This random signal was band limited to 50 Hz and had variance 0.1. Additional measurement noise was added to the sums of these two signals. This noise was band limited to 50 Hz and had a variance of 0.1. The low pass filter had a cutoff frequency of 100 rad/s concentrating the energy in the harmonics to below the 4th and third harmonic, respectively though signal was present in all harmonics up to the Nyquist frequency.
The frequency adaption gains were chosen as K a =1.95 or with frequency 7.5% to 10% of the fundamental frequencies ( Table 1). The closed loop transfer function was chosen to be a second order Chebyshev band-pass filter with 1 dB band-pass ripple, and low and high band-pass frequencies are 1 and 50 Hz, respectively. So the bandpass filter transfer function is given by  Table 2.
Under these conditions, a 50 s Matlab simulation could be performed in under 5 s. The identified frequencies are shown in Figure  3. We can see good identification and tracking of the fundamental frequencies. Figures 4 and 5 show a close up of the actual outputs of the signal generators and the identified signals (very good tracking of amplitude and relative phases). The first component has significant DC which we have not attempted to identify. In particular, there is no       way to distinguish and hence identify the DC content of the two true signals. Again we show good matches and thus we are able to identify these periodic signals in real time.
To get an overview of the signal frequency content and the accuracy of the identified models, the FFT transforms of the signal to be identified and the error signal are shown in Figure 6. It can be seen that most of harmonics have been identified although there is a huge DC component in both the signal and (e) in the proposed algorithm, which is as anticipated as we did not attempt to identify it. Figure 7 displays the quasi-periodic signal d(t), the identified signal y(t) and their difference e(t). After a brief transient we see that e becomes quite small. Note at 42 s the 5th harmonic of the 1st signal and the 4th harmonic of the second signal both had frequencies of 21.61 Hz. When a threshold of 0.01% was chosen for eliminating the redundant internal model, the algorithm went unstable. At the threshold of 0.1% there was a brief loss (<0.1 s) of performance in the signal estimation. At a threshold of 0.5% there was no noticeable loss in quality of the signal estimation.
It can be shown that the choice for controller parameters of G(s)=1, K 1ij =0 and K 2ij =1 always results in a stable feedback loop for any possible values of i ω  . Unfortunately, this leaves the dynamics of the closed loop system uncontrollable and uncertain which may require more conservative selections of the adaption gain and increase amplification of measurement noise. Figure 8 demonstrates the low performance of the simply tuned algorithm compared with the proposed algorithm. This simple approach resulted in much longer initial transient response (not shown). The steady state error was about 10 times larger in magnitude.