A recursive scheme for frequency estimation using the modulating ...
A recursive scheme for frequency estimation using the modulating functions method Giuseppe Fedele, Loredana Coluccio ∗ Dipartimento di Elettronica Informatica e Sistemistica, Universit` a degli Studi della Calabria, 87036, Rende (Cs), Italy
Abstract In this paper, an application of modulating functions method for estimation of the frequency of noisy sinusoids, is proposed. The unknown frequency is updated by introducing a recursive algorithm which is independent by the choice of the modulating functions type. The proposed recursive estimation formula is able to take into account possible abrupt changes or sweep in the frequency of the sinusoidal signal. The goodness of the proposed method is verified through numerical simulations. Key words: Frequency estimation, Modulating functions, Recursive estimation
1
Introduction
The process of estimating the frequencies of multi-sine wave signals, from a finite number of noisy discrete-time measurements, is an important task from both the theoretical and practical point of view. Such a problem has been the focus of research for quite some time and still is an active research area to date [1]–[19], since it is used in a wide range of applications in many fields such as control theory [19], relaying protection, intelligent instrumentation of power systems [1], [7], [9], signal processing [12,16–18], digital communications, distribution automation, biomedical engineering [13], radar applications, radio frequency, instrumentation and measurement, to name just a few. There is a vast amount of literature regarding the estimation procedures as ∗ Corresponding author. Tel.: +39-0984-494720; fax: +39-0984-494713. E-mail address: [email protected] (G. Fedele)
Preprint submitted to Applied Mathematics and Computation
well as the theoretical behavior of the different estimators; in [14] (and the references therein) a list of several algorithms is reported: adaptive notch filter, time frequency representation (TFR) based method, phase locked loop (PLL) based method, eigensubspace tracking estimation, extended Kalman filter frequency estimation, internal model based method (for an extensive list of references see [15]). The proposed approach uses trigonometric and spline type modulating functions [20]–[28] gaining advantages from their low-pass filtering property and gives a straightforward explicit recursion formula for the frequency estimation. The identification procedure, based on modulating functions, was pioneered by Shinbrot [20] in the 1957. Essentially, the use of modulating functions allows to transform a differential expression, involving input-output signals on a specified time interval, into a sequence of algebraic equations. Moreover, the modulating functions method annihilates the effects of initial conditions and allows the direct use of noisy data signals [27]. These features make the modulating functions method desirable for use in several real processes. In more recent years, many authors have focused on the choice of modulating functions type including Walsh functions [22], Hermite functions [24,25], Fourier modulating functions [23,27,28] and spline-type functions [21]. However, the technique of the various modulating functions method remains to be implemented non-recursively [27,28]. This paper is focused on recursive frequency estimation of a biased sinusoid corrupted by noise from a finite number of discrete-time measurements. The approach is based on a simple recursive formula which is independent by the modulating functions choice. To underline such a feature, two types of modulating functions will be used: trigonometric and spline functions. The paper is organized as follows. In Section 2, the background of the trigonometric and spline type modulating functions are discussed. A recursive estimation algorithm is outlined in Section 3. Experimental results are presented in Section 4. Section 5 is devoted to conclusions.
2
Modulating functions method
The modulating functions method can be used to estimate parameters directly from any differential equation possessing the following structure: n X
m di y(t) X di u(t) ai b = , i dti dti i=0 i=0
n ≥ m,
(1)
where y(t) and u(t) are the output and input signals respectively and {ai , bi } are the unknown system parameters. Without loss of generality let us assume a0 = 1. A function φK (t) ∈ C K (K-times differentiable), defined over a finite 2
time interval [0, T ], which satisfies the following terminal conditions
di φK (t) di φK (t) = = 0, dt t=0 dt t=T
∀i = 0, 1, ..., K − 1,
(2)
is called modulating function [26]. A function f (t) ∈ L1 over [0, T ] is modulated by taking the inner product with a modulating function φK (t) hf (t), φK (t)i =
Z
0
T
f (t)φK (t)dt.
(3)
The terminal constraints of Eq. (2) essentially make the boundary conditions of the function f (t) irrelevant after modulations. Moreover, they make possible the transfer of the differentiation operation from the function f (t) on to the modulating function φK (t), thereby eliminating the need to approximate time derivatives from noisy measurement data h
di f (t) di φK (t) i , φ (t)i = (−1) hf (t), i, K dti dti
i = 0, 1, ..., K − 1.
(4)
The modulating function procedure starts by multiplying (1) with the modulating function and integrating over the interval [0, T ]: n X
ai
i=0
Z
T
φK (ξ)
0
Z T m X di u(ξ) di y(ξ) φ (ξ) b dξ = dξ. K i dξ i dξ i 0 i=0
(5)
Integrating by parts and by using property (4), Eq. (5) becomes n X i=0
(−1)i ai
Z
0
T
Z T i m X di φK (ξ) d φK (ξ) i y(ξ)dξ = u(ξ)dξ. (−1) b i i dξ dξ i 0 i=0
(6)
To determine all parameters {ai , bi }, at least the same number of linearly independent algebraic equations like Eq. (6) must be generated. The class of spline-type modulating functions, namely φsK , implemented in this paper is shown in Fig. 1. Such a class is characterized by two parameters: the order K of the highest derivative of φsK (t), and the characteristic time T¯; for a function φsK (t), defined into interval [0, T ], T¯ = T /K. The ith derivative of a spline function, with order K and characteristic time T¯, is expressed as K P j K gji (t − j T¯), i = 0, 1, ..., K − 1, (−1) j
di φsK (t) = j=0 K i P dt (−1)j Kj δ(t − j T¯), i = K, j=0
3
(7)
where
1
gji (t − j T¯) = (K−i−1)! 0,
(t − j T¯)K−i−1 , t ∈ [j T¯, T ],
(8)
otherwise,
and δ(t) is the Dirac delta function. 10
8
6
4 φs 0 K
2
φs 1 K s2 K s3 φK s4 φ K s5 φK
φ 0
−2
−4
−6
−8
−10
0
1
2
3
4
5
6
7
8
9
10
Fig. 1. Spline functions with K = 6 and T = 10s.
The trigonometric modulating functions, as derived by Pearson and Lee (1985) [23], are linear combinations of cosine and sine functions of the first K + 1 frequency modes: φtK (t) = CF (ω, t), (9) where ω = 2π/T , F = [1, cos(ωt), − sin(ωt), ..., cos(Kωt), − sin(Kωt)]T ,
(10)
and C, in case of a single modulating function, is a row vector of dimension 2K + 1 chosen so that the conditions given in (2) are satisfied. According to [27], an explicit formulation for C, herein rearranged, is C= where
1 C¯
j+1 ¯ C(j) =γ ,K + 1 , 2
4
,
(11)
j = 1, 2, ..., K,
(12)
with
−
γ(j, k) =
k Q
l=1, l6=j
l2 , l2 −j 2
−1,
k 6= 1,
(13)
otherwise.
The choice of trigonometric functions allows a convenient matrix operation for performing differentiation, namely di φtK (t) = (−1)i CDi F, dti
(14)
D = ω diag (0, R, 2R, ..., KR) ,
(15)
where
R=
0 −1 . 1 0
(16)
Trigonometric modulating functions for K = 6 and T = 10s are shown in Fig. 2. 40
30
20 φt 0
10
K
φt 1 K t2
φK
0
t3
φK
φt 4 K
−10
φtK5
−20
−30
−40
0
1
2
3
4
5
6
7
8
9
10
Fig. 2. Trigonometric functions with K = 6 and T = 10s.
3
Identification procedure
Given the sinusoidal signal v(t) = A0 + A sin(ωt + φ), 5
(17)
it is effortless to verify that d3 v(t) dv(t) + θ =0 dt3 dt
(18)
where θ = ω 2 . The last equation can be used to get an estimation θˆ of the parameter θ. The modulation with the modulating functions permits to convert the Eq. (18) into an algebraic equation with the same mathematical form as follows Z T 3 Z T d φK (t) dφK (t) v(t)dt + θ v(t)dt = 0. (19) 3 dt dt 0 0 Since in the previous equation, the 3rd order derivative of φK (t) is involved, then the order K must be greater or equal to 4. The foregoing procedure has the primary advantage of involving the known derivatives of the modulating functions instead of the derivatives of the usually noisy output data v(t). Furthermore, since the modulating function is an integral method, it has low-pass filter properties and suppresses high-frequency components of the signal. The cut-off frequency is proportional to 1/T and out-of-range data are automatically filtered. Let consider a moving data window W(m) , m = 0, 1, ..., of n continuously arriving data of the signal v(t) sampled with a period Ts such that T = (n − 1)Ts : W(m) = [v(mTs ), v((m + 1)Ts ), ..., v((n + m − 1)Ts )] ,
m = 0, 1, ...
(20)
and define the following quantities v1(m) =
n−1 P
v((m + j)Ts )
j=0 n−1 P
v3(m) = −
j=0
dφK (t) dt t=jTs
v((m + j)Ts )
,
d3 φK (t) dt3 t=jTs
(21) .
Note that v1(m) Ts and −v3(m) Ts represent an approximation of the integrals involved in Eq. (19). Therefore, the problem consists in the determination of ˆ (where i is the generic step) which minimizes the following estimation θ(i) index i i2 h 1 X (22) µi+1−m v3(m) − v1(m) θˆi , i = 0, 1, ... J(i) = 2 m=0 where the parameter µ is a positive constant in the interval (0, 1] and represents a forgetting factor to exponentially discard the “old” data in the recursive scheme. This means that in the estimation process, the old measures have a weight smaller than the most recent measurements. The presence of a forgetting factor will make the recursive estimation scheme able to get a fast tracking in case of abrupt changes of the frequency of the signal v(t). The ˆ value of θ(i), which minimizes the criterion (22), is obtained by seeking the 6
ˆ value that cancels ∂J(i)/∂ θ(i); therefore it is equal to:
θˆi =
i P
m=0
µi+1−m v1(m) v3(m)
i P
.
m=0
(23)
µi+1−m v12(m)
A recursive algorithm of Eq. (23) is based on the following update law: µ ˆ θˆi+1 = αi θi + βi+1 , αi+1
i = 0, 1, ...
(24)
where αi =
i P
m=0
µi+1−m v12(m) ,
βi = v1(i) v3(i) ,
(25) i = 0, 1, ...
Note that αi+1 can be recursively calculated by using
αi+1 = µ αi + v12(i+1) .
(26)
Remark 1 A recursive implementation of the identification procedure in case of trigonometric modulating functions is already developed in [27] where an online parameter estimation of a continuous-time system is provided by shifting a fixed window of time series data one step forward at each sampling instant. The resulting recursion algorithm uses identities of discrete Fourier transforms along with some inherent commutative properties. Clearly the approach proposed in [27] cannot be applied to other types of modulating functions. Unlike the cited recursive method, the proposed one gives a simple recursion formula, which can be easily adapted to different modulating functions type. In the next section simulation results, by using the same approach for spline and trigonometric functions, will be presented.
4
Simulation results
In this section simulations are conducted in order to highlight the properties of the proposed method. In each experiment the recursive method described in the previous section is applied in case of spline and trigonometric type modulating functions by using the same parameters (order K of modulation, dimension of the sliding-window, forgetting factor, sampling time). The value of the parameter T for trigonometric modulating functions is chosen as the time dimension of the sliding window. 7
Example 1. The signal to be estimated is v(t) = 5 + 40 sin (2π50t + π/3) , corrupted with a gaussian white noise with signal-to-noise ratio 402 SN R = 2 = 5. 2σ The signal was sampled with a sampling time Ts = 5e − 5, the length of the sliding window was chosen as n = 600, the order of the modulating functions K = 5, and the forgetting factor µ = 1. Fig. 3 shows the on-line estimation. 51 actual estimation spline estimation trig
50.8
50.6
50.4
50.2
50
49.8
49.6
49.4
49.2
49
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0.5
Fig. 3. Example 1. On-line frequency estimations.
Example 2. This experiment is aimed to show the capability of the algorithm to take into account possible abrupt changes on the frequency of the signal. The signal to be estimated is v(t) = 220 sin (¯ ω (t)t + π/3) with
2π50,
ω ¯ (t) =
2π51,
0s ≤ t ≤ 12 , 1 2
< t ≤ 1.
The reference signal was distorted by adding the third harmonic of amplitude equal to 10% of the fundamental. The sampling time used is Ts = 1e − 4, n = 400, K = 5 and µ = 0.995. In Fig. 4 the estimated frequencies are shown. 8
52 actual estimation spline estimation trig 51.5
51
50.5
50
49.5
49
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Fig. 4. Example 2. Frequency tracking in the case of frequency step.
Example 3. In this example the signal v(t) is the same as defined in the previous example with 2π(2t + 50),
ω ¯ (t) =
2π(−2t + 52),
0s ≤ t ≤ 12 , 1 2
< t ≤ 1.
The simulation parameters remain the same except for µ = 0.95. In Fig. 5 are reported the frequency sweep and its estimations.
5
Conclusions
A recursive frequency estimation scheme using modulating function method has been proposed. Such a scheme is general since it is independent by the modulating functions type. It seems to be easy to implement and, according to the numerical simulations, it possesses good tracking capabilities in case of abrupt changes or sweep in the frequency of the signal. Moreover the method gives accurate estimation even in case of large noise on the signal. 9
52 actual estimation spline estimation trig 51.5
51
50.5
50
49.5
49
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Fig. 5. Example 3. Frequency tracking in the case of frequency sweep.
Acknowledgments The authors would like to thanks the useful advices and the constructive criticism of the anonymous referees.
References [1] J. Wu, J. Long, and J. Wang, “High-Accuracy, Wide-Range Frequency Estimation Methods for Power System Signals Under Nonsinusoidal Conditions”, IEEE Transaction on Power Delivery, 20, 1, 366-374, 2005. [2] H.C. So, K.W. Chan, Y.T. Chan, and K.C. Ho, “Linear prediction approach for efficient frequency estimation of multiple real sinusoids: algorithms and analyses”, IEEE Transactions on Signal Processing, 53, 7, 2290-2305, 2005. [3] K.W. Chan, and H.C. So, “Accurate Frequency Estimation for Real Harmonic Sinusoids”, IEEE Signal Processing Letters, 11, 7, 609-612, 2004. [4] H.C. So, and K.W. Chan, “A Generalized Weighted Linear Predictor Frequency Estimation Approach for a Complex Sinusoid”, IEEE Transactions on Signal Processing, 54, 4, 1304-1315, 2006. [5] S.S.Abeysekera, “Efficient Frequency Estimation Using the Pulse-Pair Method at Various Lags”, IEEE Transactions on Communication, 54, 9, 1542-1546, 2006.
10
[6] J.D. Klein, “Fast Algorithms for Single Frequency Estimation”, IEEE Transactions on Signal Processing, 54, 5, 1762-1770, 2006. [7] M.H. El-Shafey, and M.M. Mansour, “Application of a new frequency estimation technique to power systems”, IEEE Transactions on Power Delivery, 21, 3, 1045-1053, 2006. [8] J.R. Trapero, H. Sira-Ramirez, and V.F. Batlle, “An algebraic frequency estimator for a biased and noisy sinusoidal signal”, Signal Processing, 87, 6, 1188-1201, 2007. [9] A.K. Pradhan, A. Routray, and A. Basak, “Power System Frequency Estimation Using Least Mean Square Technique”, IEEE Transactions on Power Delivery, 20, 3, 1812-1816, 2005. [10] A. Leven, N. Kaneda, U.V. Koc, and Y.K. Chen, “Frequency Estimation in Intradyne Reception”, IEEE Photonics Technology Letters, 19, 6, 366-368, 2007. [11] C. Hackman, C. Levine, J. Parker, T.E. Piester, D. Becker, and J. Becker, “A Straightforward Frequency-Estimation Technique for GPS Carrier-Phase Time Transfer”, IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, 53, 9, 1570-1583, 2006. [12] A.P. Klapuri, “Multiple Fundamental Frequency Estimation Based on Harmonicity and Spectral Smoothness”, IEEE Transactions on Speech and Audio Processing, 11, 6, 804-816, 2003. ¨ [13] N. Ostlund, J. Yu, and J.S. Karlsson, ‘Improved maximum frequency estimation with application to instantaneous mean frequency estimation of surface electromyography’, IEEE Transactions on Biomedical Engineering, 51, 9, 15411546, 2004. [14] Q. Zhang, and L.J. Brown “Noise Analysis of an Algorithm for Uncertain Frequency Identification”, IEEE Transactions on Automatic Control, 51, 1, 103110, 2006. [15] P. Stoica, “List of references on spectral line analysis”, Signal Processing, 31, 3, 329-340, 1993. [16] A. Bonaventura, L. Coluccio, and G. Fedele, “Frequency estimation of multisinusoidal signal by multiple integrals”, 7th IEEE Int. Symp. on Sig. Proc. and Inf. Techn., Cairo, Egypt, 2007. [17] L. Coluccio, A. Eisinberg, and G. Fedele, “A Prony-like method for non-uniform sampling”, Signal Processing, 87, 2484-2490, 2007. [18] L. Coluccio, A. Eisinberg, and G. Fedele, “A property of the elementary symmetric functions on the frequencies of sinusoidal signals”, Signal Processing, 89, 765-777, 2009. [19] G. Fedele, A. Ferrise, and D. Frascino, “Multi-sinusoidal signal estimation by an adaptive sogi-filters bank”, 15th IFAC Symp. on Sys. Ident., 2009.
11
[20] M. Shinbrot, “On the analysis of linear and nonlinear systems”, Trans. ASME, 79, 547-552, 1957. [21] G.P. Rao, and L. Sivakumar, “Identification of deterministic time-lag systems”, IEEE Trans., 21, 527-529, 1976. [22] G.P. Rao, and K.R. Palanisamy, “Improved algorithm for parameter identification in continuous systems via Walsh functions”, IEE Proc. D, Contr. Theory Appl., 130, 1, 9-16, 1983. [23] A.E. Pearson, and F.C. Lee, “On the Identification of Polynomial Input-Output Differential Systems”, IEEE Trans. Autom. Contr., 30, 8, 778-782, 1985. [24] J.R. Jordan, and N.E. Paterson, “A modulating-function method for on-line fault detection”, J. Phys. E: Sci. Instrum., 19, 681-685, 1986. [25] J.R. Jordan, S.A. Jalali-Naini, and R.D.L. Mackie, “System identification with Hermite modulating functions”, IEEE Proceedings, 137, 2, 87-92, 1990. [26] H.A. Preising, and D.W.T. Rippin, “Theory and application of the modulating function method-I. Review and theory of the method and theory of the splinetype modulating functions”, Comput. Chem. Eng., 17, pp. 1–16, 1993. [27] T.B. Co, and S. Ungarala “Batch Scheme Recursive Parameter Estimation of Continuous-time System Using the Modulating Functions Method”, Automatica, 33, 6, pp. 1185–1191, 1997. [28] S. Ungarala, and T.B. Co, “Tyme-varying system identification using modulating functions and spline models with application to bio-processes”, Comput. Chem. Eng., 24, pp. 2739–2753, 2000.
12
Abstract In this paper, an application of modulating functions method for estimation of the frequency of noisy sinusoids, is proposed. The unknown frequency is updated by introducing a recursive algorithm which is independent by the choice of the modulating functions type. The proposed recursive estimation formula is able to take into account possible abrupt changes or sweep in the frequency of the sinusoidal signal. The goodness of the proposed method is verified through numerical simulations. Key words: Frequency estimation, Modulating functions, Recursive estimation
1
Introduction
The process of estimating the frequencies of multi-sine wave signals, from a finite number of noisy discrete-time measurements, is an important task from both the theoretical and practical point of view. Such a problem has been the focus of research for quite some time and still is an active research area to date [1]–[19], since it is used in a wide range of applications in many fields such as control theory [19], relaying protection, intelligent instrumentation of power systems [1], [7], [9], signal processing [12,16–18], digital communications, distribution automation, biomedical engineering [13], radar applications, radio frequency, instrumentation and measurement, to name just a few. There is a vast amount of literature regarding the estimation procedures as ∗ Corresponding author. Tel.: +39-0984-494720; fax: +39-0984-494713. E-mail address: [email protected] (G. Fedele)
Preprint submitted to Applied Mathematics and Computation
well as the theoretical behavior of the different estimators; in [14] (and the references therein) a list of several algorithms is reported: adaptive notch filter, time frequency representation (TFR) based method, phase locked loop (PLL) based method, eigensubspace tracking estimation, extended Kalman filter frequency estimation, internal model based method (for an extensive list of references see [15]). The proposed approach uses trigonometric and spline type modulating functions [20]–[28] gaining advantages from their low-pass filtering property and gives a straightforward explicit recursion formula for the frequency estimation. The identification procedure, based on modulating functions, was pioneered by Shinbrot [20] in the 1957. Essentially, the use of modulating functions allows to transform a differential expression, involving input-output signals on a specified time interval, into a sequence of algebraic equations. Moreover, the modulating functions method annihilates the effects of initial conditions and allows the direct use of noisy data signals [27]. These features make the modulating functions method desirable for use in several real processes. In more recent years, many authors have focused on the choice of modulating functions type including Walsh functions [22], Hermite functions [24,25], Fourier modulating functions [23,27,28] and spline-type functions [21]. However, the technique of the various modulating functions method remains to be implemented non-recursively [27,28]. This paper is focused on recursive frequency estimation of a biased sinusoid corrupted by noise from a finite number of discrete-time measurements. The approach is based on a simple recursive formula which is independent by the modulating functions choice. To underline such a feature, two types of modulating functions will be used: trigonometric and spline functions. The paper is organized as follows. In Section 2, the background of the trigonometric and spline type modulating functions are discussed. A recursive estimation algorithm is outlined in Section 3. Experimental results are presented in Section 4. Section 5 is devoted to conclusions.
2
Modulating functions method
The modulating functions method can be used to estimate parameters directly from any differential equation possessing the following structure: n X
m di y(t) X di u(t) ai b = , i dti dti i=0 i=0
n ≥ m,
(1)
where y(t) and u(t) are the output and input signals respectively and {ai , bi } are the unknown system parameters. Without loss of generality let us assume a0 = 1. A function φK (t) ∈ C K (K-times differentiable), defined over a finite 2
time interval [0, T ], which satisfies the following terminal conditions
di φK (t) di φK (t) = = 0, dt t=0 dt t=T
∀i = 0, 1, ..., K − 1,
(2)
is called modulating function [26]. A function f (t) ∈ L1 over [0, T ] is modulated by taking the inner product with a modulating function φK (t) hf (t), φK (t)i =
Z
0
T
f (t)φK (t)dt.
(3)
The terminal constraints of Eq. (2) essentially make the boundary conditions of the function f (t) irrelevant after modulations. Moreover, they make possible the transfer of the differentiation operation from the function f (t) on to the modulating function φK (t), thereby eliminating the need to approximate time derivatives from noisy measurement data h
di f (t) di φK (t) i , φ (t)i = (−1) hf (t), i, K dti dti
i = 0, 1, ..., K − 1.
(4)
The modulating function procedure starts by multiplying (1) with the modulating function and integrating over the interval [0, T ]: n X
ai
i=0
Z
T
φK (ξ)
0
Z T m X di u(ξ) di y(ξ) φ (ξ) b dξ = dξ. K i dξ i dξ i 0 i=0
(5)
Integrating by parts and by using property (4), Eq. (5) becomes n X i=0
(−1)i ai
Z
0
T
Z T i m X di φK (ξ) d φK (ξ) i y(ξ)dξ = u(ξ)dξ. (−1) b i i dξ dξ i 0 i=0
(6)
To determine all parameters {ai , bi }, at least the same number of linearly independent algebraic equations like Eq. (6) must be generated. The class of spline-type modulating functions, namely φsK , implemented in this paper is shown in Fig. 1. Such a class is characterized by two parameters: the order K of the highest derivative of φsK (t), and the characteristic time T¯; for a function φsK (t), defined into interval [0, T ], T¯ = T /K. The ith derivative of a spline function, with order K and characteristic time T¯, is expressed as K P j K gji (t − j T¯), i = 0, 1, ..., K − 1, (−1) j
di φsK (t) = j=0 K i P dt (−1)j Kj δ(t − j T¯), i = K, j=0
3
(7)
where
1
gji (t − j T¯) = (K−i−1)! 0,
(t − j T¯)K−i−1 , t ∈ [j T¯, T ],
(8)
otherwise,
and δ(t) is the Dirac delta function. 10
8
6
4 φs 0 K
2
φs 1 K s2 K s3 φK s4 φ K s5 φK
φ 0
−2
−4
−6
−8
−10
0
1
2
3
4
5
6
7
8
9
10
Fig. 1. Spline functions with K = 6 and T = 10s.
The trigonometric modulating functions, as derived by Pearson and Lee (1985) [23], are linear combinations of cosine and sine functions of the first K + 1 frequency modes: φtK (t) = CF (ω, t), (9) where ω = 2π/T , F = [1, cos(ωt), − sin(ωt), ..., cos(Kωt), − sin(Kωt)]T ,
(10)
and C, in case of a single modulating function, is a row vector of dimension 2K + 1 chosen so that the conditions given in (2) are satisfied. According to [27], an explicit formulation for C, herein rearranged, is C= where
1 C¯
j+1 ¯ C(j) =γ ,K + 1 , 2
4
,
(11)
j = 1, 2, ..., K,
(12)
with
−
γ(j, k) =
k Q
l=1, l6=j
l2 , l2 −j 2
−1,
k 6= 1,
(13)
otherwise.
The choice of trigonometric functions allows a convenient matrix operation for performing differentiation, namely di φtK (t) = (−1)i CDi F, dti
(14)
D = ω diag (0, R, 2R, ..., KR) ,
(15)
where
R=
0 −1 . 1 0
(16)
Trigonometric modulating functions for K = 6 and T = 10s are shown in Fig. 2. 40
30
20 φt 0
10
K
φt 1 K t2
φK
0
t3
φK
φt 4 K
−10
φtK5
−20
−30
−40
0
1
2
3
4
5
6
7
8
9
10
Fig. 2. Trigonometric functions with K = 6 and T = 10s.
3
Identification procedure
Given the sinusoidal signal v(t) = A0 + A sin(ωt + φ), 5
(17)
it is effortless to verify that d3 v(t) dv(t) + θ =0 dt3 dt
(18)
where θ = ω 2 . The last equation can be used to get an estimation θˆ of the parameter θ. The modulation with the modulating functions permits to convert the Eq. (18) into an algebraic equation with the same mathematical form as follows Z T 3 Z T d φK (t) dφK (t) v(t)dt + θ v(t)dt = 0. (19) 3 dt dt 0 0 Since in the previous equation, the 3rd order derivative of φK (t) is involved, then the order K must be greater or equal to 4. The foregoing procedure has the primary advantage of involving the known derivatives of the modulating functions instead of the derivatives of the usually noisy output data v(t). Furthermore, since the modulating function is an integral method, it has low-pass filter properties and suppresses high-frequency components of the signal. The cut-off frequency is proportional to 1/T and out-of-range data are automatically filtered. Let consider a moving data window W(m) , m = 0, 1, ..., of n continuously arriving data of the signal v(t) sampled with a period Ts such that T = (n − 1)Ts : W(m) = [v(mTs ), v((m + 1)Ts ), ..., v((n + m − 1)Ts )] ,
m = 0, 1, ...
(20)
and define the following quantities v1(m) =
n−1 P
v((m + j)Ts )
j=0 n−1 P
v3(m) = −
j=0
dφK (t) dt t=jTs
v((m + j)Ts )
,
d3 φK (t) dt3 t=jTs
(21) .
Note that v1(m) Ts and −v3(m) Ts represent an approximation of the integrals involved in Eq. (19). Therefore, the problem consists in the determination of ˆ (where i is the generic step) which minimizes the following estimation θ(i) index i i2 h 1 X (22) µi+1−m v3(m) − v1(m) θˆi , i = 0, 1, ... J(i) = 2 m=0 where the parameter µ is a positive constant in the interval (0, 1] and represents a forgetting factor to exponentially discard the “old” data in the recursive scheme. This means that in the estimation process, the old measures have a weight smaller than the most recent measurements. The presence of a forgetting factor will make the recursive estimation scheme able to get a fast tracking in case of abrupt changes of the frequency of the signal v(t). The ˆ value of θ(i), which minimizes the criterion (22), is obtained by seeking the 6
ˆ value that cancels ∂J(i)/∂ θ(i); therefore it is equal to:
θˆi =
i P
m=0
µi+1−m v1(m) v3(m)
i P
.
m=0
(23)
µi+1−m v12(m)
A recursive algorithm of Eq. (23) is based on the following update law: µ ˆ θˆi+1 = αi θi + βi+1 , αi+1
i = 0, 1, ...
(24)
where αi =
i P
m=0
µi+1−m v12(m) ,
βi = v1(i) v3(i) ,
(25) i = 0, 1, ...
Note that αi+1 can be recursively calculated by using
αi+1 = µ αi + v12(i+1) .
(26)
Remark 1 A recursive implementation of the identification procedure in case of trigonometric modulating functions is already developed in [27] where an online parameter estimation of a continuous-time system is provided by shifting a fixed window of time series data one step forward at each sampling instant. The resulting recursion algorithm uses identities of discrete Fourier transforms along with some inherent commutative properties. Clearly the approach proposed in [27] cannot be applied to other types of modulating functions. Unlike the cited recursive method, the proposed one gives a simple recursion formula, which can be easily adapted to different modulating functions type. In the next section simulation results, by using the same approach for spline and trigonometric functions, will be presented.
4
Simulation results
In this section simulations are conducted in order to highlight the properties of the proposed method. In each experiment the recursive method described in the previous section is applied in case of spline and trigonometric type modulating functions by using the same parameters (order K of modulation, dimension of the sliding-window, forgetting factor, sampling time). The value of the parameter T for trigonometric modulating functions is chosen as the time dimension of the sliding window. 7
Example 1. The signal to be estimated is v(t) = 5 + 40 sin (2π50t + π/3) , corrupted with a gaussian white noise with signal-to-noise ratio 402 SN R = 2 = 5. 2σ The signal was sampled with a sampling time Ts = 5e − 5, the length of the sliding window was chosen as n = 600, the order of the modulating functions K = 5, and the forgetting factor µ = 1. Fig. 3 shows the on-line estimation. 51 actual estimation spline estimation trig
50.8
50.6
50.4
50.2
50
49.8
49.6
49.4
49.2
49
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0.5
Fig. 3. Example 1. On-line frequency estimations.
Example 2. This experiment is aimed to show the capability of the algorithm to take into account possible abrupt changes on the frequency of the signal. The signal to be estimated is v(t) = 220 sin (¯ ω (t)t + π/3) with
2π50,
ω ¯ (t) =
2π51,
0s ≤ t ≤ 12 , 1 2
< t ≤ 1.
The reference signal was distorted by adding the third harmonic of amplitude equal to 10% of the fundamental. The sampling time used is Ts = 1e − 4, n = 400, K = 5 and µ = 0.995. In Fig. 4 the estimated frequencies are shown. 8
52 actual estimation spline estimation trig 51.5
51
50.5
50
49.5
49
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Fig. 4. Example 2. Frequency tracking in the case of frequency step.
Example 3. In this example the signal v(t) is the same as defined in the previous example with 2π(2t + 50),
ω ¯ (t) =
2π(−2t + 52),
0s ≤ t ≤ 12 , 1 2
< t ≤ 1.
The simulation parameters remain the same except for µ = 0.95. In Fig. 5 are reported the frequency sweep and its estimations.
5
Conclusions
A recursive frequency estimation scheme using modulating function method has been proposed. Such a scheme is general since it is independent by the modulating functions type. It seems to be easy to implement and, according to the numerical simulations, it possesses good tracking capabilities in case of abrupt changes or sweep in the frequency of the signal. Moreover the method gives accurate estimation even in case of large noise on the signal. 9
52 actual estimation spline estimation trig 51.5
51
50.5
50
49.5
49
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Fig. 5. Example 3. Frequency tracking in the case of frequency sweep.
Acknowledgments The authors would like to thanks the useful advices and the constructive criticism of the anonymous referees.
References [1] J. Wu, J. Long, and J. Wang, “High-Accuracy, Wide-Range Frequency Estimation Methods for Power System Signals Under Nonsinusoidal Conditions”, IEEE Transaction on Power Delivery, 20, 1, 366-374, 2005. [2] H.C. So, K.W. Chan, Y.T. Chan, and K.C. Ho, “Linear prediction approach for efficient frequency estimation of multiple real sinusoids: algorithms and analyses”, IEEE Transactions on Signal Processing, 53, 7, 2290-2305, 2005. [3] K.W. Chan, and H.C. So, “Accurate Frequency Estimation for Real Harmonic Sinusoids”, IEEE Signal Processing Letters, 11, 7, 609-612, 2004. [4] H.C. So, and K.W. Chan, “A Generalized Weighted Linear Predictor Frequency Estimation Approach for a Complex Sinusoid”, IEEE Transactions on Signal Processing, 54, 4, 1304-1315, 2006. [5] S.S.Abeysekera, “Efficient Frequency Estimation Using the Pulse-Pair Method at Various Lags”, IEEE Transactions on Communication, 54, 9, 1542-1546, 2006.
10
[6] J.D. Klein, “Fast Algorithms for Single Frequency Estimation”, IEEE Transactions on Signal Processing, 54, 5, 1762-1770, 2006. [7] M.H. El-Shafey, and M.M. Mansour, “Application of a new frequency estimation technique to power systems”, IEEE Transactions on Power Delivery, 21, 3, 1045-1053, 2006. [8] J.R. Trapero, H. Sira-Ramirez, and V.F. Batlle, “An algebraic frequency estimator for a biased and noisy sinusoidal signal”, Signal Processing, 87, 6, 1188-1201, 2007. [9] A.K. Pradhan, A. Routray, and A. Basak, “Power System Frequency Estimation Using Least Mean Square Technique”, IEEE Transactions on Power Delivery, 20, 3, 1812-1816, 2005. [10] A. Leven, N. Kaneda, U.V. Koc, and Y.K. Chen, “Frequency Estimation in Intradyne Reception”, IEEE Photonics Technology Letters, 19, 6, 366-368, 2007. [11] C. Hackman, C. Levine, J. Parker, T.E. Piester, D. Becker, and J. Becker, “A Straightforward Frequency-Estimation Technique for GPS Carrier-Phase Time Transfer”, IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, 53, 9, 1570-1583, 2006. [12] A.P. Klapuri, “Multiple Fundamental Frequency Estimation Based on Harmonicity and Spectral Smoothness”, IEEE Transactions on Speech and Audio Processing, 11, 6, 804-816, 2003. ¨ [13] N. Ostlund, J. Yu, and J.S. Karlsson, ‘Improved maximum frequency estimation with application to instantaneous mean frequency estimation of surface electromyography’, IEEE Transactions on Biomedical Engineering, 51, 9, 15411546, 2004. [14] Q. Zhang, and L.J. Brown “Noise Analysis of an Algorithm for Uncertain Frequency Identification”, IEEE Transactions on Automatic Control, 51, 1, 103110, 2006. [15] P. Stoica, “List of references on spectral line analysis”, Signal Processing, 31, 3, 329-340, 1993. [16] A. Bonaventura, L. Coluccio, and G. Fedele, “Frequency estimation of multisinusoidal signal by multiple integrals”, 7th IEEE Int. Symp. on Sig. Proc. and Inf. Techn., Cairo, Egypt, 2007. [17] L. Coluccio, A. Eisinberg, and G. Fedele, “A Prony-like method for non-uniform sampling”, Signal Processing, 87, 2484-2490, 2007. [18] L. Coluccio, A. Eisinberg, and G. Fedele, “A property of the elementary symmetric functions on the frequencies of sinusoidal signals”, Signal Processing, 89, 765-777, 2009. [19] G. Fedele, A. Ferrise, and D. Frascino, “Multi-sinusoidal signal estimation by an adaptive sogi-filters bank”, 15th IFAC Symp. on Sys. Ident., 2009.
11
[20] M. Shinbrot, “On the analysis of linear and nonlinear systems”, Trans. ASME, 79, 547-552, 1957. [21] G.P. Rao, and L. Sivakumar, “Identification of deterministic time-lag systems”, IEEE Trans., 21, 527-529, 1976. [22] G.P. Rao, and K.R. Palanisamy, “Improved algorithm for parameter identification in continuous systems via Walsh functions”, IEE Proc. D, Contr. Theory Appl., 130, 1, 9-16, 1983. [23] A.E. Pearson, and F.C. Lee, “On the Identification of Polynomial Input-Output Differential Systems”, IEEE Trans. Autom. Contr., 30, 8, 778-782, 1985. [24] J.R. Jordan, and N.E. Paterson, “A modulating-function method for on-line fault detection”, J. Phys. E: Sci. Instrum., 19, 681-685, 1986. [25] J.R. Jordan, S.A. Jalali-Naini, and R.D.L. Mackie, “System identification with Hermite modulating functions”, IEEE Proceedings, 137, 2, 87-92, 1990. [26] H.A. Preising, and D.W.T. Rippin, “Theory and application of the modulating function method-I. Review and theory of the method and theory of the splinetype modulating functions”, Comput. Chem. Eng., 17, pp. 1–16, 1993. [27] T.B. Co, and S. Ungarala “Batch Scheme Recursive Parameter Estimation of Continuous-time System Using the Modulating Functions Method”, Automatica, 33, 6, pp. 1185–1191, 1997. [28] S. Ungarala, and T.B. Co, “Tyme-varying system identification using modulating functions and spline models with application to bio-processes”, Comput. Chem. Eng., 24, pp. 2739–2753, 2000.
12
