EURASIP Journal on Applied Signal Processing 2005:8, 1251–1260 © 2005 Hindawi Publishing Corporation Fast Frequency Estimation by Zero Crossings of Differential Spline Wavelet Transform

Zero crossings or extrema of a wavelet transform constitute important signatures for signal analysis with the advantage of great simplicity. In this paper, we introduce a fast frequency-estimation method based on zero-crossing counting in the transform domain of a family of differential spline wavelets. The resolution and order of the vanishing moments of the chosen wavelets have a close relation with the frequency components of a signal. Theoretical results on estimating the highest and the lowest frequency components are derived, which are particularly useful for frequency estimation of harmonic signals. The results are illustrated with the help of several numerical examples. Finally, we discuss the connection of this approach with other frequency estimation methods, with the high-order level-crossing analysis in statistics, and with the scaling theorem in computer vision.


INTRODUCTION
Zero crossings or extrema of a signal constitute very important features that have been used in signal processing applications for detection of weak signals, search for periodicity, white-noise testing, and so on, because of their great simplicity [1]. They have also found applications in many engineering fields. For example, the higher-order zero crossings were utilized for fast detection of contractions in uterine electromyography [2] and for discrimination of discontinuous breath sounds [3]. In computer vision, zero crossings were usually used for detection of edges in an image [4]. In wavelet theory, zero crossings at multiple scales were used as critical signatures from which a signal can be reconstructed [5,6].
The counting of the number of zero crossings is a "dirty" but fast approximate way to estimate frequency without having to resort to spectral analysis calculations. It provides an alternative to the traditional spectral methods and has been implemented in several existing software packages. Specifically, the frequency of a pure sinusoidal wave can be directly estimated by simply counting the number of times the waveform crosses the point of zero amplitude. Even though zero crossings have been used in the wavelet domain for signal reconstruction, it remains unclear how they are related to the spectrum of a signal. An intuitive observation is that a signal of high frequency has more zero crossings in the wavelet transform and vice versa. In addition, with the increase of scale, the number of zero crossings will be decreased because of the smoothing operation. Thus, it motivates us to explore the intrinsic connection between the wavelet signatures of a signal and its spectrum. We have found that a mathematically rigorous characterization of their connection can be easily established using a special class of differential wavelets introduced in [7]. In order to keep the paper self-contained, we will briefly review the differential spline wavelets in Section 2. In Section 3, we will introduce the notations and preliminary results for zero-crossing estimation. In Section 4, the main results will be derived pertaining to the connection of zero crossings of a differential spline wavelet transform with the frequency of a signal. The illustrations of these theoretical results will be presented with the help of several numerical examples. We then conclude the paper with a discussion on the connection of the proposed method with other frequency estimation approaches, with high-order level-crossing analysis in statistics, and with the scaling theorem in computer vision.

Definitions
The differential spline wavelet transforms are defined as a class of continuous wavelet transforms when the wavelets are chosen as the derivatives of B-spline functions [7]. The continuous wavelet transform of a signal X(t) is defined as where ψ s (x) = (1/s)ψ(x/s) ∈ L 2 (R) is the wavelet at scale s. The mother wavelet ψ(t) is usually taken as the kth-order derivative of a B-spline function of order n, that is, where β n (x) is the continuous B-spline of order n generated by the repeated n + 1 convolution of a B-spline of order 0, and β 0 (x) is the 0th-order B-spline or the rectangular pulse. The discrete sampled B-spline b n m (k) of order n and integer coarseness m ≥ 1 is obtained by directly sampling the nth-order continuous B-spline at scale m: The discrete B-spline of order n at scale m is defined by where B 0 m = (1/m)[1, 1, . . . , 1] is the sampled pulse of width m.
The relation between the continuous B-spline β n (x) and the discrete B-splines is given by the following m-scale relation [8]:

Filter bank implementation
When wavelets are defined as the derivatives of B-splines, the wavelets can be written as a linear combination of translated B-spline bases where operators g are the difference operators and β n2 is the B-spline of order n 2 (here we use n 2 to distinguish it from the B-spline of order n 1 in (11)). Specifically, according to definition (2), the kth-order derivative filters can be derived as where l j = l!/ j!(l − j)! is the binomial coefficient. The differential filters of order l are simply the lth difference of the binomials and the corresponding wavelet has the vanishing moment of l. In the frequency domain, it has the transfer function G(ω) = (2i sin(ω/2)) l .
We list a few of the most often used difference filters of lower orders.
(i) g (1) is the first-order difference operator. The corresponding finite impulse responses (FIRs) are (2) is the second-order difference operator. The corresponding FIRs are g (2) (−1) = 1, g (2) (0) = −2, g (2) (1) = 1, Since the spline spaces S n h provide close and stable approximations of L 2 (R), the discrete signal can be represented using B-spline bases. We use the translated B-splines of order n 1 as the bases to approximate the signal By using the m-refinable property (6) of the B-splines, we can derive a cascaded filter bank algorithm for fast implementation of a continuous wavelet transform at rational scales [8]: where ↑ m and ↓ m denote upsampling and downsampling operations by a factor of m. The computational complexity of the above algorithm lies in the convolutions with B n1 m2 and B n2 m1 , which can be fast implemented with only additions [8].
When the scale is dyadic, the computation in (12) can be realized using a two-scale recursive formula where W (m) j+1 X(k) denotes the mth-order differential wavelet transform. It is the mth difference of a smoothed sequence. In this case, the above process becomes a binary summation and difference process. The smoothing filter {h(k)} is a binomial filter, that is, its corresponding FIRs are In the frequency domain, it has the transfer function H(ω) = (cos(ω/2)) k .

Notations and definitions
Let {X t }, t = . . . , −1, 0, 1, . . . , be a discrete-time random signal or process. If it is stationary, its covariance function will only depend on the lag l, where µ = E(X t ) is the mean of X t . For simplicity and without loosing generality, we assume that µ = 0 throughout the paper. When l = 0, γ l becomes the variance The autocorrelation function is obtained from normalizing γ l by γ 0 : It is obvious that |ρ l | ≤ 1.
The zero crossings of a stationary Gaussian time series {X t } are associated with a clipped binary series {Z t } defined in [9], that is, If we define an indicator function at time t, then d t is either 0 or 1. If d t = 1, a zero crossing occurs at time t. The total number of zero crossings in {X t } is given by the sum where ZC#(X t ) denotes the number of zero crossings in sequence X t .

Zero-crossing counting for continuous spline wavelet transform
The density of local extrema D for a stationary Gaussian process f can be obtained based on the second-and fourthorder derivatives of its autocorrelation function R, or equivalently, based on the second-and fourth-order moments of its spectrum P(ω) [10]: where the autocorrelation functions are related to the spectrum p(ω) by Using this relation, we have derived a formula for the estimation of the average number of zero crossings for white noise, or more generally, for a Brownian motion (fBm). For an fBm signal f , its power spectrum is inversely proportional to the power of its frequency, that is, the 1/ f law. More specifically, where λ is a constant and 0 < H < 1 is the Hurst parameter, which controls the "roughness" of the fBm [11]. The realization of fBm is fractals with dimension d = 2−H. The smaller the Hurst parameter, the more rough or higher fractional dimension of a curve. We have proven that the density of zero crossing in the continuous spline wavelet transform has the following relationship with the scale and the order of spline [12].

Theorem 1. For a continuous B-spline wavelet transform, the density of zero crossing of a fractal Brownian motion at scale s is
From this theorem, we know that the density of zero crossings can be used to estimate the dimension of a fractal signal. It indicates that the roughness of a fractal curve can be determined by the number of its zero crossings. A more rough curve with smaller Hurst parameter H will have the higher density of zero crossings.
If we take the scale s to be dyadic, s = 2 j , we can derive that where c is a constant. It implies that the zero-crosssing number at the next scale is reduced to be half of the current one, and at the log scale, the zero crossing number of an fBm is linearly decreasing with dyadic scale j.

Zero-crossing count for discrete filtering
Zero-crossing analysis can be traced as early as in [13]. The analysis of the statistical properties of zero crossings in a more general discrete-time setting was pioneered by the work in [14,15].
Lemma 1. For a zero-mean stationary Gaussian random signal {X t }, t = 1, . . . , N, there exists the following relation: Equivalently, the average zero-crossing number E(ZC#(X t ))/(N − 1) is given by cos −1 ρ 1 /π, the known cosine formula. This lemma establishes the connection between the expected number of zero crossings of a random signal and its autocorrelation function ρ 1 . According to the Wiener-Khintchine theorem [16], it follows that where F is called the spectral distribution function of the random signal. This formula thus gives the spectral distribution of the zero crossings. Its continuous-time analog is given by the formula in (22). If {X t } is filtered by a linear filter £ with transfer function H, the output is still a Gaussian random signal with zero mean and its spectrum is given by |H(ω)| 2 dF(ω). The zerocrossing number of the filtered signal {F(X t )} will be related by [10] as follows: . (29)

Connection of the wavelet scale with frequency components
The connection of the scale or resolution of its differential spline wavelet transform with the frequency components of a signal is given by the following theorem.
Theorem 2. For a zero-mean stationary Gaussian random sig- let ω * be the lowest frequency in the spectrum. Then the number of zero crossings of the discrete spline scale-space filtering in (13) decreases with the scale and in the limit case as j → ∞, it converges to the lowest frequency This result is a more general extension to the theorem in [1,15], which holds for any smoothing filtering. To prove this theorem, we first prove the lemma described below.
Lemma 2. For a zero-mean stationary Gaussian random signal {X t }, t = 1, . . . , N, and letting Y t = λ 1 X t + λ 2 X t−1 for any real positive λ 1 and λ 2 , then Proof. Letρ 1 be the correlation function of the summed process {Y t }. By the definition and the stationariness of the sequence, we havē and so We know that for any real a, b, the following inequality holds [9]: According to the definitions in (16), (17), and (18), expansion of the left-hand side of this inequality with normalization leads to Therefore, we have It follows that Because the relation in (27) and the cosine function is monotonously decreasing within [0, π], we obtain the inequity in (32).
We now come to prove Theorem 2. From formula (13), we see that S j+1 is obtained by convolution with a binomial filter followed by a downsampling of 2. In the simplest case, the filter is (1/2, 1/2). Thus according to Lemma 2, the number of zero crossings of S j+1 is decreasing with the scale parameter j. Also, it is easy to verify that ZC#(S j+1 ) is bounded by 0. Therefore, according to the functional theory, it must converge to its lower bound, which is ω * . From the process of this proof, it can be seen that, although we prove the theorem when the scale is dyadic, it still holds if the scale is chosen to be rational.
The asymptotic or large sample properties of a certain type of estimator are usually desirable in statistical inference. One such property is the simple consistency of a sequence of estimators. A sequence of estimators {T n } is said to be the simple consistent estimators of parameter θ if for every ε > 0, lim n→∞ P[|T n − θ| < ε] = 1 for every θ ∈ Ω. Another desired property of estimators is the asymptotic unbiasedness. A sequence of estimators {T n } is said to be asymptotically unbiased for parameter θ if lim n→∞ E[T n ] = θ for every θ ∈ Ω. With these established concepts, we give the properties of the estimators for the zero crossings in the following corollary. Proof. The proof of the corollary follows from Theorem 2 and the convergence theorems in probability [17].

Connection of the wavelet vanishing moments with frequency components
From (14), we see that an mth-order differential wavelet transform is obtained by taking the difference of its (m−1)th order wavelet transform. Thus we have the following relation: that is, the wavelet transform with vanishing moments m is a difference operation on the one with vanishing moments m−1. In such a case, it can be proven that the autocorrelation function of a difference sequence Z t = X t − X t−1 has the form Therefore, Following the same process, we can obtain similar results regarding the effect on the estimation of spectrum using derivative filters of different orders.
Theorem 3. For a zero-mean stationary Gaussian random signal {X t }, t = 1, . . . , N, let ω * be the highest frequency in the spectrum. Then the number of zero crossings of the wavelet filtering in (14) increases with the order of derivative or vanishing moments and in the limit case as j → ∞, it converges to the highest frequency This theorem indicates that when one increases the derivative order of a spline wavelet filtering or vanishing moments of a spline wavelet transform, the corresponding zerocrossing number approximates a certain high frequency.
Similar to Corollary 1, we can prove the following properties of this estimator. Because of formula (41), we can infer that if the number of zero crossings remains unchanged, the signal is a sinusoidal. In this case ρ(1) = ρ (1) , which implies that the following difference equation holds according to (41): The solution of this equation is a sinusoidal.

Corollary 3. For a Gaussian zero-mean stationary process
In this extreme case, we can see that a sinusoidal signal can be completely reconstructed by its zero-crossing counts.
The estimation of frequency components of a harmonic signal has been encountered in many applications [18]. Theorem 2 provides a way to estimate the lowest frequency components of a signal by increasing the order of splines and the scales, while Theorem 3 can be used to estimate the highest frequency components by increasing the order of derivatives or vanishing moments. The combination of these two theorems suggests an algorithm to estimate the frequency components of a harmonic signal.
An algorithm for frequency estimation of a harmonic signal (a) Estimate the highest frequency components by using the sufficient high order of the vanishing moments of differential wavelets; (b) estimate the lowest frequency components by using high-scale smoothing; (c) perform the filtering operation in the frequency domain to filter these two lowest and highest components; (d) repeat the above procedures (a-c) to estimate the next lowest and highest frequency components, until all the in-between components are obtained.

NUMERICAL EXAMPLES
We give the following numerical examples to illustrate the above theorems. In practice, because the signal is usually contaminated by noise, it is reasonable to assume that the signal is random. However, for the illustration purpose, we use the deterministic harmonic signals in the following examples. It can be seen that although the above theorems are established for random signals, we can show that they also hold for deterministic signals.

Example 1
To illustrate the above theorem, suppose that we have a sinusoidal signal composed of three different frequency components as shown in Figure 1: If we perform the cubic spline filtering at scale 0 (without smoothing) and then take the difference up to order 9 or the wavelets with vanishing moments from 1 to 9 in (2). Table 1 shows the zero-crossing counts and the frequency estimates by using πE(ZC#(W ( j) ))/(N − 1) as in (43). We see that with the increase of differential order, the highest frequency component can be estimated, as proven in Theorem 3. The middle frequency 1.5 can be estimated by the first-and second-order derivatives.
If we want to estimate the lowest frequency component, according to Theorem 2, we should perform smoothing or take the wavelet transform at higher scales. Table 2 provides the results for the wavelet transform at scale 1-9.
The theorems in Section 4 also suggest that a combination of smoothing and derivatives can be used to estimate certain frequencies between the lowest and highest frequency components. As an example, we look at the frequency estimates at smoothing scales 5 and 10, and the derivatives from 1 to 9 in Table 3. It can be seen that at lower derivatives, we obtain the estimates of lower frequencies. With the increase of the order of derivatives, the frequency 1.5 can be estimated.
This signal is plotted in Figure 2. Table 4 shows the frequency estimates for the cases when the vanishing moments of wavelets go from 1 to 9. Table 5 shows the frequency estimates for the cases when the smoothing scales vary from 1 to 9. With the smoothing, the noise is weakened, and so the lowest frequency can be computed.
If we perform the smoothing at scale 5 and then take the wavelets up to the order of 9 , we are able to estimate the frequency components between the lowest and highest ones. Table 6 demonstrates such an example.    If we take the vanishing moments of the wavelet to be infinite, according to Theorem 3, it will approximate to the highest frequency of the signal. In such an example, it will be π because the spectrum of white noise is from 0 to π. Table 7 verifies that the highest frequency is π when the vanishing moments of wavelets go up to 118.

Example 3
We test Corollary 3 by a simple example. Suppose that we have a sinusoidal signal with frequency 0.45, as plotted in Figure 3. We consider the cubic-spline wavelet transform with up to 10 vanishing moments or with the 1st-10th-order difference of cubic-spline filtering. For each wavelet transform, we present the zero-crossing number counts and the corresponding frequency estimates using formula (43). As shown in Table 8, the zero-crossing number remains unchanged except for numerical errors and the frequency is approximately estimated.
If we perform the spline wavelet filtering at different scales and followed by the same operations described above, the results are shown in Table 9. Due to smoothing, the low frequency components are smoothed out and the estimation becomes inaccurate for higher derivatives.

DISCUSSIONS
In this paper, we have investigated how the number counting of zero crossings of a wavelet transform can be used to estimate the frequency components of signals in a fast and simplistic fashion. It is well known that the FFT is the most commonly used tool to estimate a signal's spectrum, with a computational complexity of O(N log N), where N is the length of the signal. Under the approach proposed in this paper, the computation of differential wavelet transforms and ( (19), (20), and (21)) can in fact all be implemented with additions only. The complexity associated with such an approach is simply linear, and therefore much lower than that of the FFT approach. The estimation of signal frequency has also been investigated in [19], where a continuous Gaborbased wavelet transform was used. However this latter approach appears to be computationally much more demanding because of the need to compute the continuous wavelet transform.
The study of the properties of zero crossings roots from the pioneering work of Rice [13] and Logan [20]. Its application to signal processing was pioneered by Kedem [1]. The use of the zero crossings of a wavelet transform or scale-space filtering has been reported in many papers but the connection with the signal's frequency have been seldom addressed.     Figure 3: A sinusoidal signal.
In this work, we show that the zero crossings in the differential spline wavelet domain can be used to effectively approximate the frequencies of a signal. Higher frequency can be estimated by wavelets with higher order of vanishing moments, while lower frequency can be estimated by increasing the smoothing scale. Proper use of the smoothing scale and differential order can facilitate the estimation of certain frequency components. The relationship between the wavelet transform signatures of a signal and its frequency estimation based on Fourier transform is thus established. Formula (27) also offers a way to compute the autocorrelation, which is often used in signal spectrum analysis [2,3,21].
In the scale-space theory for computer vision, the zero crossings in multiscale or scale-space filtering are usually referred to as the fingerprints in the scale-space theory [22,23]. The famous scaling theorem states that the Gaussian kernel is the only filter with which the number of zero crossings of a signal does not increase with the increase of the smoothing scale [24]. This property is also known as the monotone or the embedding property that has been proven for the continuous signal case. In practice, however, the kernels used for scale-space filtering are not Gaussian. For example, some of the most commonly used kernels are spline approximations [8,25], or more general scaling filters in the wavelet theory [26]. It has been shown that for discrete signals, the above embedding property still holds for certain discrete smoothing kernels [27]. Whether the embedding property still holds for more general discrete smoothing filters, in particular the spline filters, has not been proven previously. The results presented in this paper attempt to provide a rigorous proof of this property for a random signal in the discrete setting.

ACKNOWLEDGMENTS
The authors would like to express their sincere thanks to the two anonymous reviewers for their detailed comments. In particular, the correction of a proof error and the suggestion of additional references are greatly appreciated. The first author wishes to thank Professor Akram Aldroubi of the Mathematics Department of Vanderbilt University for conversations and comments on this work. This work is supported in part by the NIH under the SBIR Grants 5R44 HD33658-03 and 1R43RR16817-01.