- Email: [email protected]
Analytical formulation of the fractal dimension of filtered stochastic signals
ARTICLE IN PRESS Signal Processing 90 (2010) 2690–2697
Contents lists available at ScienceDirect
Signal Processing journal homepage: www.elsevier.com/locate/sigpro
Analytical formulation of the fractal dimension of filtered stochastic signals Jean-Marc Girault a,, Denis Kouame´ b, Abdeldjalil Ouahabi c a b c
Universite´ Franc- ois Rabelais de Tours, UMRS ‘‘Imaging and Brain’’ U930 and CNRS ERL 3106, Tours, France Universite´ PAUL Sabatier de Toulouse 3, IRIT UMR 5505, France Universite´ Franc- ois Rabelais de Tours, CERTIHS, Tours, France
a r t i c l e i n f o
abstract
Article history: Received 10 September 2009 Received in revised form 12 February 2010 Accepted 15 March 2010 Available online 24 March 2010
The aim of this study was to investigate the effects of a linear filter on the regularity of a given stochastic process in terms of the fractal dimension. This general approach, described in a continuous time domain, is new and is characterized by its simplicity. The framework of this problem is general since it emerges when a fractal process undertakes a transformation, as is the case in denoising or measurement processes. & 2010 Elsevier B.V. All rights reserved.
Keywords: Length Fractal dimension Stochastic signals Lorentzian filters
1. Introduction and definition 1.1. Introduction Fractal processes are processes whose ubiquity spans innumerable fields. In contrast to the well-known family of ARMA processes, these intrinsically scale-invariant processes, also known as 1/f processes, are characterized by long-term correlation structure. Moreover, these fractal signals are continuous and underivable signals presenting deterministic or statistical self-affinities [27,29]. Invariants such as the fractal dimension D, the ¨ Hurst exponent H or the Holder exponent h are often evaluated to characterize such irregular signals. Stochastic signals can be continuous and not differentiable, presenting strong irregularities. Our study was performed in the context of analysis and synthesis of fractal signals where fractal and stochastic signals can be related through fractional
Corresponding author.
E-mail addresses: [email protected], [email protected] (J.-M. Girault). 0165-1684/$ - see front matter & 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.sigpro.2010.03.019
differentiation [16] or filtering [26]. Others have studied 1/f noise in electronic devices [3,4,18], roughness and texture of surfaces [12,24,35,17] and biomedical signals and images [33,22,30,20]. Issues dealing with filtering and fractal signals have previously been addressed in three ways: the denoising of fractal signals corrupted by an additive noise [23,34,8,6,7], the filtering of fractal signals by a measurement system [28,12,30,22,35,17,20] and the synthesis of fractal signals from stochastic processes [32,26,10,19,14]. In the synthesis of 1=f b noise, the framework which predominated for many years was naturally the spectral domain [3,18,26,32]. The relationships between the spectral and irregularity properties of the process appeared later in reports by Berry [1] who demonstrated the link between the slope b of the spectrum and the fractal dimension D ¼ ð5bÞ=2 for a fBm. The spectral estimator consisted of evaluating the slope b in a log– log plot (power/frequency), and then applying the previously introduced Berry equation. Although much used in the past, due to its simple formulation (for a continuous or a discrete framework), this spectral estimator is less often used today. The main reasons for its obsolescence are its poor performance and to the need for
ARTICLE IN PRESS J.-M. Girault et al. / Signal Processing 90 (2010) 2690–2697
setting the spectral bandwidth of interest when the power spectral density presents several slopes. The approaches based on maximum likelihood [23,34,7] and the estimators based on multi-scale analysis [10,34,5,15] were mainly used in denoising of fractal signals (or images). Because the likelihood function was a non-linear function according to D, it was not possible to obtain an explicit formula and only a numerical computation was possible. For the approaches based on wavelet analysis, they were mainly focused on estimation of the slope of variance b ¼ 2H þ 1 of the wavelet coefficients dj in the log–log plot. This estimator was very interesting since it could easily be extended to the study of multi-fractal processes. When the fractal signals (mono or multi-dimensional) are accessible only by measurement, as is the case for example when measuring surface roughness [12,35,17] or image texture [30,22,20], the measurement or acquisition systems act like a linear filter. Quantification of the impact of such linear filtering on the regularity of the processes studied is thus necessary. For these applications, the most commonly used fractal dimension estimators have been based on the box counting method (for example see [33]), the spectral slope (for example see [11]) and length measurement [2]. There is not, to our knowledge, a simple framework which directly shows the link between the filter characteristics and the regularity property of such filtered signals through explicit relationships. We therefore
2691
investigated the impact of linear filtering on irregularity for a given stochastic process in terms of the fractal dimension. The novelty of our approach is based on an explicit formulation, in continuous time, linking the spectral properties of a given filter and the fractal properties of a given process. Our approach shows the relationship between the fractal dimensions at the input and the output of a filter. Due to its simple formulation and the good performance of its estimator, we focused on the fractal dimension formulation through the length measurement. The only a priori knowledge required to use our approach is the probability density of the input process and the correlation coefficient of the filter. Our approach is original since the issue has never been previously addressed in denoising, synthesis or measuring processes of fractal signals. Our study is complementary to other studies based mostly on a discrete formulation. The general framework and the simplicity of our formulation are its main assets (Fig. 1). Before explaining the details of our approach, we summarize the definitions used related to the length measurement of fractal stochastic processes. 1.2. Definition The Euclidian length of a deterministic signal x(t) for a limited duration T can be defined as follows: L0 ¼ T l ,
ð1Þ
Fractal signal (H=0.5) and filtered fractal signal 1
amplitude
0.5 0 −0.5 −1 −1.5 0
200
400
0
200
400
600
800
1000
600
800
1000
1
Amplitude
0.5 0 −0.5 −1 −1.5 sample Fig. 1. Low pass filtering increases the regularity of the fractal signal. (a) Fractal signal (H= 0.5). (b) Filtered fractal signal.
ARTICLE IN PRESS J.-M. Girault et al. / Signal Processing 90 (2010) 2690–2697
where Z 1 T=2 l¼ lðtÞ dt T T=2
0.2 0
ð2Þ
The segment λ is variable and the projector τ is constant.
−0.2
is the arithmetic mean of the elementary length and where sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 dxðtÞ lðtÞ ¼ 1 þ : ð3Þ dt On the hypothesis that x(t) is homogeneous to Volt, then l is homogeneous to Volt per second. The elementary length l can be rewritten according to the Pythagorean theorem as qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 ð4Þ lðtÞ ¼ lim t2 þ Dðt, tÞ2 , t-0 t
amplitude
2692
−0.4 −0.6 λ
−0.8 −1 −1.2
τ
−1.4 30
35
40
45
50
time Fig. 2. Principle of length measurement of a stochastic signal.
where the increment Dðt, tÞ is written in a symmetric form as
Dðt, tÞ ¼ xðt þ t=2Þxðtt=2Þ:
ð5Þ
Finally, the signal length is written as Z T=2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 Lx ¼ lim t2 þ Dðt, tÞ2 dt:
The definition of Eq. (8) extends easily to stationary and ergodic stochastic signals by replacing the arithmetic mean by the mathematical expectation E[ ]. The definition then becomes
ð6Þ
Lx ðlÞ ¼ TE½l,
t-0 T=2
t
Note that we had already used such a definition for deterministic signals (see [13]). Unlike deterministic signals, which are almost always derivable and have a finite length, stochastic signals do not follow the laws of Euclidian geometry. Stochastic signals cannot be described by continuous mathematical curves and their irregularities are still present even when examined with higher and higher resolutions. This leads to the paradoxical conclusion that the lengths of such signals are not finite. To overcome the fact that (as fractal signals) stochastic signals do not have a finite length, Richardson [31] assumed that the length Lx of a fractal signal depends on the resolution at which it is measured, i.e. the gauge l of the measurement. To overcome this indetermination for l-0, Richardson proposed measuring the curve length for each value l, 8l 40. Mandelbrot demonstrated that measuring Lx ðlÞ as a function of l was one way to determine the fractal dimension: 1Dx
Lx ðlÞ ¼ Kx l
,
ð7Þ
where Dx is the fractal dimension and Kx is a constant. Length definition becomes a function of l and is written like Eq. (1): Lx ðlÞ ¼ T lðtÞ,
ð8Þ
where
lðtÞ ¼
1 T
Z
T=2
lðt, tÞ dt,
ð9Þ
T=2
and where the ffi elementary length is qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi lðt, tÞ ¼ ð1=tÞ t2 þ Dðt, tÞ2 . The corresponding length estimator of Eq. (8) illustrated in Fig. 2 is Nk qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 X ðkT s Þ2 þ Dðn,kT s Þ2 , ð10Þ L^ x ðkT s Þ ¼ kT s n ¼ 1 where Ts is the sampling period and N is the number of elementary lengths between ½T=2,T=2.
ð11Þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 where l ¼ ð1=tÞ t2 þ D . Finally, the definition becomes qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi T ð12Þ t2 þ D2 : Lx ðtÞ ¼ E
t
The calculation of the mathematical expectation of Eq. (12) does not lead to an analytical formulation, and only a numerical resolution is possible. In order to present an analytical formulation, we propose using another definition based on that proposed by Burlaga and Klein [2]: Lx ðtÞ ¼
T
t
E½jDj:
ð13Þ
This definition supposes that t 5 D, which is often true, particularly for a high sampling frequency. However, if this is not the case, an amplitude gain can be applied. A discussion of the sense of such a definition is reported in [21]. The length estimator that we finally used which corresponded to a discrete version of Eq. (13) was Lx ðkÞ ¼
Nk X
jDðn,kÞj:
n¼1
The performance of such estimator in its average version, known as the Burlaga and Klein [2] estimator LBK(k), was excellent and outperformed the spectral estimators and the maximum likelihood estimators (see [11]). However, the performance was similar to that obtained with estimators based on wavelet analysis (see [11]). Note that Flandrin [10] demonstrated the link between length offfiq the wavelet measurement LBK(k) and variance ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi coeffipffiffiffiffiffiffiffiffi cient dj(k): E½LBK ð2j1 Þ ¼ 235j=2 2=p varðdj ðkÞÞ. Though this estimator, which is based on wavelet analysis, gave an excellent performance, it did not lead to a simple formulation in a continuous framework. The strong point of this estimator based on wavelet analysis was its adaptability to study multi-fractal processes.
ARTICLE IN PRESS J.-M. Girault et al. / Signal Processing 90 (2010) 2690–2697
The paper is organized as follows. Section 2 presents the approach in detail, setting out the notion of length ratio of colored processes, then stating the link between the length ratio and the fractal dimension of colored processes. Section 3 presents several examples of filters and the simulation results are discussed. Finally we propose a conclusion and some future applications.
t
E½jx1 x2 j,
2 T sx qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Lx ðtÞ ¼ pffiffiffiffi 1rx ðtÞ:
p t
ð20Þ
Similarly the length of the output signal y(t) is written as 2 T sy qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Ly ðtÞ ¼ pffiffiffiffi 1ry ðtÞ: ð21Þ
p t
Lx ðtÞ ¼
The method that we propose is new and was investigated in a general stochastic framework with an explicit formulation in continuous time. This approach requires some a priori information such as the probability density of the input process and the correlation function of the filter. In a linear system where the impulse response is h(t), the input signal x(t) and the output signal y(t), and the length Lx ðtÞ of the input signal, can be evaluated by T
and finally
Note that for a white noise ðrx ¼ 0Þ, the length of the input signal x(t) is written as
2. Length ratio and fractal dimension of filtered stochastic signals
Lx ðtÞ ¼
2693
ð14Þ
T sx
t
g,
ð22Þ
pffiffiffiffi where the factor g is equal to 2= pffiffiffi p for a Gaussian noise and the factor g is equal to 2= 3 for a noise uniformly distributed pffiffiffiffiffiffi between ½a=2,a=2 with a standard deviation sx ¼ a= 12. The factor g (and hence the length) seems to reveal the statistic nature of the white noise. In order to establish a link between the impulse response or the correlation function of the filter ry ðtÞ and the fractal dimension Dy of the filtered signal y(t), it should be remembered that in the case of fractal signals, the length Ly of the output signal is related to the fractal dimension Dy as follows (see also Eq. (7) for the input): Ly ðtÞ ¼ Ky t1Dy ,
where x1 ¼ xðt1 þ t=2Þ, x2 ¼ xðt2 t=2Þ 8t1 ,t2 and where Z Z jx1 x2 jpx1 ,x2 ðx1 ,x2 Þ dx1 dx2 : ð15Þ E½jx1 x2 j ¼
where Ky is a constant. Consequently, the fractal dimension Dy of the signal y(t) can be defined as follows:
px1 ,x2 ðx1 ,x2 Þ is the joint probability density.
Dy ¼ 1lim
In the special case where x1 and x2 are two random Gaussian variables, the joint probability density px1 ,x2 ðx1 ,x2 Þ is written as ð1r2x Þ1=2 ð1=ð2ð1r2x ÞÞÞðx2 =s2 þ x2 =s2 2rx x1 x2 =s1 s2 Þ 1 1 2 2 e , 2ps1 s2 ð16Þ
where s21 and s22 are the variances of the two random variables x1 and x2, respectively. rx is the correlation coefficient defined by rx ¼ E½x1 x2 =E½x21 E½x22 . The probability density of the increment D ¼ x2 x1 is written as Z þ1 px1 x2 ðx1 ,x1 þ DÞ dx1 : ð17Þ pD ðDÞ ¼ 1
After calculation, the probability density of the increment is reduced to a normal law Nð0, s0 2 Þ: 2 1 02 ð18Þ pD ðDÞ ¼ pffiffiffiffiffiffi eD =2s , 2ps0 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi with s0 ¼ s21 þ s22 2rx s1 s2 . x1 and x2 belong to the ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sx and the standard same signal x(t), then s1 ¼ s2 p deviation s0 is simplified, s0 ¼ sx 2ð1rx Þ. The mathematical expectation of the increment becomes rffiffiffiffi Z þ1 2 0 E½jDj ¼ 2 DpD ðDÞ dD ¼ s, ð19Þ
0
logðLy ðtÞÞ : logðtÞ
ð24Þ
logðLx ðtÞÞ : logðtÞ
ð25Þ
Similarly:
2.1. Filtered Gaussian noise
px1 ,x2 ðx1 ,x2 Þ ¼
t-0
ð23Þ
p
Dx ¼ 1lim
t-0
By using Eqs. (20), (25) and (21), (24) the equations become Dy ¼ 2lim
1 logð1ry ðtÞÞ , logðtÞ
ð26Þ
Dx ¼ 2lim
1 logð1rx ðtÞÞ : logðtÞ
ð27Þ
t-0 2
t-0 2
For unidimensional signal, the fractal dimension of signals x(t) and y(t) was the topological dimension of the space minus a term related to the inner correlation of x(t) and y(t). In order to simply qualify the influence of linear filters, we propose evaluating a new parameter, i.e. the length ratio O between the input and the output of the linear system:
OðtÞ ¼
Ly ðtÞ : Lx ðtÞ
By using Eqs. (20) and (21), Eq. (28) becomes sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sy 1ry ðtÞ : OðtÞ ¼ sx 1rx ðtÞ
ð28Þ
ð29Þ
Furthermore, the following function LðtÞ is defined as
LðtÞ ¼ lim
t-0
logðOðtÞÞ : logðtÞ
ð30Þ
LðtÞ is an important parameter since it quantifies the filtering impact in terms of reduction irregularity. By
ARTICLE IN PRESS 2694
J.-M. Girault et al. / Signal Processing 90 (2010) 2690–2697
using Eqs. (7), (23) and (28), it is easy to show that
L ¼ Dx Dy :
2.3. Impact of the correlation coefficient ð31Þ
In this case Dy is defined by Dy ¼ Dx lim
t-0
logðOðtÞÞ : logðtÞ
ð32Þ
When the input signal is a Gaussian noise and when the system is linear, by using Eqs. (29) and (32), this becomes 1 logð1ry ðtÞÞ 1 logð1rx ðtÞÞ þ lim : Dy ¼ Dx lim t-0 2 t-0 2 logðtÞ logðtÞ ð33Þ The fractal dimension of the filtered signal is then equal to the fractal dimension of the signal before filtering minus a term related to the filter plus a term explaining the inner correlation of the input signal x(t). Before testing different types of filter, we present the filtered fractional Brownian motion.
From definitions (33) and (37), it becomes obvious that the asymptotic behavior of the correlation coefficient when t-0 is a key point in the estimation of fractal dimension. Indeed, by exciting filters by a white Gaussian noise, all correlation coefficients ry ðtÞ whose Taylor expansion satisfies the equation limr ðtÞ 1kt t-0 y
n
,
ð38Þ
where k is a known constant and where n 2 R and 0 rn r 2, involve a fractal dimension and a Hurst exponent described as follows: 8 n > < D y ¼ Dx , 2 ð39Þ > H ¼ H þ n, : y x 2 where Dx = 2 and Hx =0 for rx ¼ 0. Note that for 2 Z Dy Z 1, n satisfies 0 r n r2. Moreover, by using definition (33) and Eq. (38) for rx ¼ 0, it becomes 1 logðktn Þ n ¼ 2 , t-0 2 logðtÞ 2
2.2. Filtered fractional Brownian motion
Dy ¼ 2lim
Let us consider the fractional Brownian motion (fBm) [25] as being in the input of the filter. This self-similar process which generalizes the Brownian motion or the Wiener process is composed of stationary and independent increments. These increments follow a Gaussian law1: x1 x2 Nð0,2t2Hx Þ where Hx is the Hurst exponent. Note that an average power spectrum density of fBm can be defined in the frequency domain [9] and it turns out that it is proportional to 1=f a with a ¼ 2Hx þ1. By using Eqs. (14)–(20), we show that
and n= 4 2Dy. Without making out an exhaustive list, we propose the following correlation coefficients (or a combination) satisfying Eq. (38) for t 4 0:
2 Lx ðtÞ ¼ pffiffiffiffi T tHx 1 :
p
ð40Þ
ry ðtÞ ¼ a þ ð1aÞcosðktn=2 Þ. For a ¼ 0:5 and for n= 2 we recognize the Hanning function;
ry ðtÞ ¼ ð1ktn Þ. For n= 1 we recognize the triangle function; n
ry ðtÞ ¼ ekt . For n =2 we recognize the Gaussian function and for n = 1 we recognize the Lorentzian function.
ð34Þ
pffiffiffiffi Note that we recognize the factor 2= p which underlines the Gaussian nature of the process. By using Eqs. (7)–(25), it becomes
3. Examples
Dx ¼ 2Hx ,
The impulse response h(t) of a Lorentzian filter is defined by
ð35Þ
where logðLx ðtÞÞ : Hx ¼ 1þ lim t-0 logðtÞ
3.1. Lorentzian filters
hðtÞ ¼ oc eoc t cosðo0 tÞuðtÞ, ð36Þ
Note that Hy has the same definition as Hx. As for Eq. (33), it becomes 1 logð1ry ðtÞÞ : ð37Þ Hy ¼ Hx lim t-0 2 logðtÞ Similarly, the Hurst exponent Hy of the filter output is the Hurst exponent Hx before filtering minus a term related to the inner filter properties. 1 Note that the increments follow a Gaussian law written as follows: x1 x2 Nð0, t2Hx Þ. However, in order to make a direct link between equations in this subsection and equations in Section 2.1, we voluntary modified this formula by doubling the variance and this modification does not modify the definition of fBm.
ð41Þ
where u(t) is the Heaviside function, oc the cut-off pulsation and o0 the central pulsation of the filter. This linear system of impulse response h(t) is excited by a Gaussian white noise. For linear filters, the crosscorrelation function of the output signal is written as Ry ðtÞ ¼ Rh ðtÞ Rx ðtÞ,
ð42Þ
where is the convolution operator. Note that the process coming from this type of filter is also known as being the ‘‘first order Markov process’’. For the Lorentzian filter the correlation function is Rh ðtÞ ¼ ðoc =4Þeoc jtj cosðo0 tÞ, the correlation function of the input signal is Rx ðtÞ ¼ s2x dðtÞ, and the correlation function of the output signal reduces to Ry ðtÞ ¼ s2y eoc jtj cosðo0 tÞ, where pffiffiffiffiffiffi sy ¼ sx oc =2. The correlation coefficient defined by ry ðtÞ ¼ Ry ðtÞ=Ry ð0Þ, is reported in Table 1 for t 40.
ARTICLE IN PRESS J.-M. Girault et al. / Signal Processing 90 (2010) 2690–2697
The output signal is a white noise, i.e. an irregular signal which fills all the space.
Table 1 Correlation coefficients and length ratios 8t 4 0. Filters
h(t)
ry ðtÞ
Wide band oc eoc t cosðo0 tÞuðtÞ eoc t cosðo0 tÞ cosðo0 tÞuðtÞ
Narrow band Base band All pass
cosðo0 tÞ
oc eoc t uðtÞ
eoc t
dðtÞ
0
OðtÞ=G pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1eoc t cosðo0 tÞ o tpffiffiffi 0 sin 2 2 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1eoc t 1
When the input is a fractional Brownian motion, the Hurst exponent of the output signal is:
for a base band filter o0 -0 and for t-0, oc t 5 1: 1 logðoc tÞ : t-0 2 logðtÞ
Hy ¼ Hx þ lim
Where rx ¼ 0 8t 40 and where the gain is G ¼ sy =sx .
As already mentioned, the length of the white Gaussian noise is 8t 4 0: 2T sx Lx ðtÞ ¼ pffiffiffiffi t,
p
ð43Þ
2T sy pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Ly ðtÞ ¼ pffiffiffiffi t 1eoc t cosðo0 tÞ:
ð44Þ
The length ratio 8t 40 is written as follows: pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi OðtÞ ¼ G 1eoc t cosðo0 tÞ,
ð45Þ
p
where the gain G is written as G ¼ sy =sx . In order to test different kinds of Lorentzian filter, four filters were considered: a wide band filter, a narrow band filter, a base band filter and an all pass filter. The correlation coefficients and length ratios for these different kinds of filter are summarized in Table 1. We verified that for these different kinds of filter, the parameter O=G was related to the intrinsic filter parameters, oc and o0 . When the input is a white Gaussian noise, the fractal dimension Dy of the output is 1 logð1eoc t cosðo0 tÞÞ : ð46Þ Dy ¼ 2lim t-0 2 logðtÞ
By imposing Hx = 0.3 (high irregularity), the Hurst exponent of the output is Hy ¼ 0:3 þ0:5 ¼ 0:8. This confirms the fact that the greater the filtering the greater the regularity; for a narrow band filter oc -0 and for t-0, o0 t 5 1: 1 logðo20 t2 =2Þ : logðtÞ t-0 2
since the correlation coefficient is rx ¼ 0, 8t 4 0. The length of the output signal is 8t 4 0:
Hy ¼ Hx þ lim
By imposing Hx = 0.3 (high irregularity), the Hurst exponent of the output is Hy ¼ 0:3 þ 1 ¼ 1:3. As 0 r Hy r 1, this means that Hy =1. 3.2. Simulations As an illustration, Figs. 3 and 4 show the length ratios for the four types of Lorentzian filter. Time series of 10 000 samples are simulated with a sampling frequency of 10 000 Hz for different values of the cut-off frequency fc ¼ f1,10,100,1000,5000g Hz and for different values of the central frequency of the filter f0 ¼ f1,10g Hz. The minimum and maximum values of the gauge are tmin ¼ 104 s and tmax ¼ 0:5 s. The sampling frequency is thus an important setting parameter since it fixes the minimum resolution. Fig. 3 shows that for a cut-off frequency fc =5000 Hz and a central frequency f0 =1 Hz (similar to an all pass
Length ratio of a low pass filtered white gaussian noise (f0=1 Hz)
For a base band filter o0 -0 and for t-0, oc t 51: Dy ¼ 212 ¼ 32,
2695
ð47Þ 100
with
1/ω c
D=2
limð1eoc t Þ oc t:
t-0
For a narrow band filter oc -0 and for t-0, o0 t 5 1: ð48Þ
with limð1cosðo0 tÞÞ o20 t2 =2:
t-0
The output signal is a sinusoid, i.e. a continuous signal whose fractal dimension is the same as a line. Note that for a Gaussian correlation coefficient 2 ry ðtÞ ¼ eoc t , the fractal dimension is also Dy =1 since 2 limt-0 ð1eoc t Þ oc t2 . For an all pass filter: 1 logð1Þ ¼ 20 ¼ 2: Dy ¼ 2lim t-0 2 logðtÞ
Ω/G
Dy ¼ 21 ¼ 1,
10−1 D=1.5
10−4
10−3
10−2
10−1
gauge (τ) Fig. 3. Length ratio of colored stochastic signals. First order Lorentzian low pass filter for different cut-off pulsations and for a fixed central pulsation.
ARTICLE IN PRESS 2696
J.-M. Girault et al. / Signal Processing 90 (2010) 2690–2697
filter), the length ratio O of such a signal has a slope equal to zero, involving a fractal dimension Dy =2. The output signal is a white Gaussian noise which entirely fills the space. Note that for a white noise, we verified that Oðtmin Þ=G 1. When the cut-off frequency was reduced to fc =1 Hz and f0 =1 Hz (low pass filter), the slope for the length ratio was 1/2 for a gauge lower than the cut-off gauge t o tc ¼ 1=oc , involving a fractal dimension Dy =3/2. The output signal was a standard Brownian motion characterized by its Hurst exponent H= 1/2. Beyond the cut-off gauge, the slope was equal to zero. Fig. 4 shows that for a central frequency f0 =100 Hz and for a cut-off frequency fc = 10 Hz, oscillations where the position of its first minimum corresponds to t ¼ 1=f0 appear in O. For a gauge less than the cut-off gauge tc ¼ 1=oc , the fractal dimension is between 1.5 and 1 for a narrow band filter and between 2 and 1.5 for a wide band filter. When the cut-off frequency tends towards zero, the fractal dimension tends to 1, thus involving the presence of a sinusoidal signal. When the cut-off frequency tends towards infinity (whatever the central frequency) then
Length ratio of band pass filtered white gaussian noise (f0=100 Hz) D=2
100
D=1.5
Ω/G
1/f
0
10−1 D=1
10−4
10−3
10−2
10−1
Gauge (τ) Fig. 4. Length ratio of colored stochastic signals. First order Lorentzian low pass filter for different cut-off pulsations and for a fixed central pulsation.
Synthesized fractal signals (n−"Lorentzian filter")
Length ratio of filtered stochastic signals
50 D=1.5, n=1 45
D=2, n=0 0
10 40
30 Ω/G
justapoxed signals
35
25
20
15
−1
10
10
5
D=1.5, n=1 D=2, n=0
0 0
1000
2000 3000 sample
4000
5000
−3
10
−2
10
−1
10
gauge (τ)
Fig. 5. Length ratio of colored stochastic signals. (a) Synthesized fractal signals with an n-order ‘‘Lorentzian filter’’. (b) Log–log plot of length ratio O=G versus gauge t.
ARTICLE IN PRESS J.-M. Girault et al. / Signal Processing 90 (2010) 2690–2697
the fractal dimension tends to 2, involving the presence of a white noise. n By using an impulse response of the form hðtÞ ¼ eoc t with 0 rn r 1, time series of 5000 samples are simulated with a sampling frequency of 10 000 Hz and a cut-off frequency fc =1 Hz. The input signal is a white Gaussian noise N(0,1). The minimum and maximum values of the gauge are tmin ¼ 104 s and tmax ¼ 0:5 s. Simulation results are depicted in Fig. 5. Fig. 5a shows synthesized fractal signals: the higher the order n, the greater the regularity of the time series. Fig. 5b shows the length ratio O=G with respect to the gauge t: the higher the order n, the higher the slope and the lower the fractal dimension. For n= 0, the corresponding filter is an all pass filter, the output signal is a white Gaussian noise and the corresponding fractal dimension is Dy = 2. For n = 1, the corresponding filter is a base band filter, the output signal is a Brownian motion and the corresponding fractal dimension is Dy = 1.5 (Hy = 0.5). For 0 r n r1 the corresponding fractal dimensions of the output filter are between 2 ZDy Z 1:5. 4. Discussion and conclusion We propose here a simple explicit formulation of the filtering effect on the regularity of a given stochastic process in a general framework and in a continuous time context. The simulations studied validated the theoretical calculus but also permitted extraction of the filter characteristics such as cut-off frequency and central frequency. Though limited to mono-fractal processes, this approach is interesting since it requires only commonly used a priori knowledge such as the probability law of the process and the spectral characteristic of the filter. The need for a priori information such as the probability density and the correlation coefficient of the filter constitute both an advantage and a drawback. It is an advantage since they provide a simple formulation, but it is also a disadvantage since the study presented here was limited to this information. Another drawback is that our approach is currently limited to mono-fractal processes. The next step will be to extend our framework to the study of multi-fractal processes. References [1] M.V. Berry, Diffractals, J. Phys. A Math. Gen. 12 (1979) 781. [2] L. Burlaga, L. Klein, Fractal structure of the interplanetary magnetic field, J. Geophys. Res. 91 (1986) 347–350. [3] J.A. Barnes, D.W. Allan, A statistical model of flicker noise, Proc. IEEE 54 (1966) 176–185. [4] J.A. Barnes, S. Jarvis, Efficient numerical and analog modeling of flicker noise processes, National Bureau of Standard Technical Note 604, 1971. [5] B.S. Chen, G.W. Lin, Multiscale Wiener filter for the restoration of fractal signals: wavelet filter bank approach, IEEE Trans. Signal Process. 42 (1994) 2972–2982. [6] B.S. Chen, W.S. Hou, Deconvolution filter design for fractal signal transmission systems: a multiscale Kalman filter bank approach, IEEE Trans. Signal Process. 45 (1997) 1359–1364.
2697
[7] B.S. Chen, B.K. Lee, S.C. Peng, Maximum likelihood parameter estimation of F-ARIMA processes using the genetic algorithm in the frequency domain, IEEE Trans. Signal Process. 50 (2002) 2208–2220. [8] M. Deriche, A.H. Tewfik, Signal modeling with filtered discrete fractional noise processes, IEEE Trans. Signal Process. 41 (1993) 2839–2849. [9] P. Flandrin, On the spectrum of fractional Brownian motions, IEEE Trans. Inform. Theory 35 (1989) 197–199. [10] P. Flandrin, Wavelet analysis and synthesis of fractional Brownian motion, IEEE Trans. Inform. Theory 38 (1992) 910–917. [11] N. Gache, P. Flandrin, D. Garreau, Fractal dimension estimators for fractional Brownian motions, in: Proceedings of the ICASSP, 1991, pp. 3557–3560. [12] J. Gagnepain, C. Roques-Carmes, Fractal approach to two-dimensional and three-dimensional surface roughness, in: International Conference on Metrology and Properties of Engineering Surfaces, vol. 109, 1986, pp. 19–126. [13] J.M. Girault, D. Kouame´, Length and frequency of band-limited signals, IEEE Signal Process. Lett. 9 (2002) 371–374. [14] M. Guglielmi, 1=f a signal synthesis with precision control, Signal Processing 86 (2006) 2548–2553. [15] G.A. Hirchoren, C.E. D’Attellis, On the optimal number of scales in estimation of fractal signals using wavelets and filter banks, Signal Processing 63 (1997) 55–63. [16] J.R.M. Hosking, Fractional differencing, Biometrika 68 (1981) 165–176. [17] J. James, Fracture surface examination of dental ceramics using fractal analysis, Dental Mater. 21 (2005) 586–589. [18] M.S. Keshner, 1/f noise, Proc. IEEE 70 (1982) 211–219. [19] N.J. Kasdin, Discrete simulation of colored noise and stochastic process and 1/f power law noise generation, Proc. IEEE 83 (1995) 802–826. [20] W.-L. Lee, K.-S. Hsieh, A robust algorithm for the fractal dimension of images and its applications to the classification of natural images and ultrasonic liver images, Signal Processing 9 (6) (2010) 1894–1904. [21] E. Le Tavernier, La me´thode de Higuchi pour la dimension fractal, Signal Processing 65 (1998) 115–128. [22] T. Loussot, R. Harbar, G. Jacquet, C.L. Benhamou, E. Lespesailles, A. Julien, An oriented fractal analysis for the characterization of texture: application to bone radiographs, Eusipco Signal Processing 1 (1996) 371–374. [23] T. Lundahl, W.J. Ohley, S.M. Kay, R. Siffert, Fractional Brownian motion: a maximum likelihood estimator and its application to image texture, IEEE Trans. Med. Imaging 5 (1986) 151–161. [24] A. Majumdar, B. Bhushan, Fractal model of elastic–plastic contact between rough surfaces, J. Tribol. 113 (1991) 1–11. [25] B.B. Mandelbrot, J.W. Van Ness, Fractional Brownian motions, fractional noises and applications, SIAM Rev. 10 (1968) 422–437. [26] B.B. Mandelbrot, A fast fractional Gaussian noise generator, Water Resource Res. 7 (1971) 543–553. [27] B. Mandelbrot, The Fractal Geometry of Nature, Freeman, San Francisco, 1982. [28] B.B. Mandelbrot, D.E. Passoja, A.J. Paulay, Fractal character of fracture surfaces of metals, Nature 308 (1984) 721–722. [29] B. Mandelbrot, Self-affine fractals and fractal dimension, Phys. Scr. 32 (1986) 257–260. [30] P. Qingrui, Y. Zhuangzhi, Mammogram analysis based on wavelets and fractals, in: Proceedings of IEEE Engineering in Medicine and Biology Society, 1996, pp. 1159–1162. [31] L. Richardson, The problem of contiguity: an appendix of statistic of deadly quarrels, General Systems Yearbook 61 (1961) 139–187. [32] R. Saletti, A comparison between two methods to generate 1/f noise, Proc. IEEE 74 (1986) 1595–1600. [33] N. Sarkar, B.B. Chaudhuri, An efficient differential box-counting approach to compute fractal dimension of image, IEEE Trans. Syst. Man Cybernet. 24 (1994) 115–120. [34] G.W. Wornell, A.V. Oppenheim, Estimation of fractal signals from noisy measurements using wavelets, IEEE Trans. Signal Process. 40 (1992) 611–622. [35] J.J. Wu, Characterization of fractal surfaces, Wear 239 (2000) 36–47.
Contents lists available at ScienceDirect
Signal Processing journal homepage: www.elsevier.com/locate/sigpro
Analytical formulation of the fractal dimension of filtered stochastic signals Jean-Marc Girault a,, Denis Kouame´ b, Abdeldjalil Ouahabi c a b c
Universite´ Franc- ois Rabelais de Tours, UMRS ‘‘Imaging and Brain’’ U930 and CNRS ERL 3106, Tours, France Universite´ PAUL Sabatier de Toulouse 3, IRIT UMR 5505, France Universite´ Franc- ois Rabelais de Tours, CERTIHS, Tours, France
a r t i c l e i n f o
abstract
Article history: Received 10 September 2009 Received in revised form 12 February 2010 Accepted 15 March 2010 Available online 24 March 2010
The aim of this study was to investigate the effects of a linear filter on the regularity of a given stochastic process in terms of the fractal dimension. This general approach, described in a continuous time domain, is new and is characterized by its simplicity. The framework of this problem is general since it emerges when a fractal process undertakes a transformation, as is the case in denoising or measurement processes. & 2010 Elsevier B.V. All rights reserved.
Keywords: Length Fractal dimension Stochastic signals Lorentzian filters
1. Introduction and definition 1.1. Introduction Fractal processes are processes whose ubiquity spans innumerable fields. In contrast to the well-known family of ARMA processes, these intrinsically scale-invariant processes, also known as 1/f processes, are characterized by long-term correlation structure. Moreover, these fractal signals are continuous and underivable signals presenting deterministic or statistical self-affinities [27,29]. Invariants such as the fractal dimension D, the ¨ Hurst exponent H or the Holder exponent h are often evaluated to characterize such irregular signals. Stochastic signals can be continuous and not differentiable, presenting strong irregularities. Our study was performed in the context of analysis and synthesis of fractal signals where fractal and stochastic signals can be related through fractional
Corresponding author.
E-mail addresses: [email protected], [email protected] (J.-M. Girault). 0165-1684/$ - see front matter & 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.sigpro.2010.03.019
differentiation [16] or filtering [26]. Others have studied 1/f noise in electronic devices [3,4,18], roughness and texture of surfaces [12,24,35,17] and biomedical signals and images [33,22,30,20]. Issues dealing with filtering and fractal signals have previously been addressed in three ways: the denoising of fractal signals corrupted by an additive noise [23,34,8,6,7], the filtering of fractal signals by a measurement system [28,12,30,22,35,17,20] and the synthesis of fractal signals from stochastic processes [32,26,10,19,14]. In the synthesis of 1=f b noise, the framework which predominated for many years was naturally the spectral domain [3,18,26,32]. The relationships between the spectral and irregularity properties of the process appeared later in reports by Berry [1] who demonstrated the link between the slope b of the spectrum and the fractal dimension D ¼ ð5bÞ=2 for a fBm. The spectral estimator consisted of evaluating the slope b in a log– log plot (power/frequency), and then applying the previously introduced Berry equation. Although much used in the past, due to its simple formulation (for a continuous or a discrete framework), this spectral estimator is less often used today. The main reasons for its obsolescence are its poor performance and to the need for
ARTICLE IN PRESS J.-M. Girault et al. / Signal Processing 90 (2010) 2690–2697
setting the spectral bandwidth of interest when the power spectral density presents several slopes. The approaches based on maximum likelihood [23,34,7] and the estimators based on multi-scale analysis [10,34,5,15] were mainly used in denoising of fractal signals (or images). Because the likelihood function was a non-linear function according to D, it was not possible to obtain an explicit formula and only a numerical computation was possible. For the approaches based on wavelet analysis, they were mainly focused on estimation of the slope of variance b ¼ 2H þ 1 of the wavelet coefficients dj in the log–log plot. This estimator was very interesting since it could easily be extended to the study of multi-fractal processes. When the fractal signals (mono or multi-dimensional) are accessible only by measurement, as is the case for example when measuring surface roughness [12,35,17] or image texture [30,22,20], the measurement or acquisition systems act like a linear filter. Quantification of the impact of such linear filtering on the regularity of the processes studied is thus necessary. For these applications, the most commonly used fractal dimension estimators have been based on the box counting method (for example see [33]), the spectral slope (for example see [11]) and length measurement [2]. There is not, to our knowledge, a simple framework which directly shows the link between the filter characteristics and the regularity property of such filtered signals through explicit relationships. We therefore
2691
investigated the impact of linear filtering on irregularity for a given stochastic process in terms of the fractal dimension. The novelty of our approach is based on an explicit formulation, in continuous time, linking the spectral properties of a given filter and the fractal properties of a given process. Our approach shows the relationship between the fractal dimensions at the input and the output of a filter. Due to its simple formulation and the good performance of its estimator, we focused on the fractal dimension formulation through the length measurement. The only a priori knowledge required to use our approach is the probability density of the input process and the correlation coefficient of the filter. Our approach is original since the issue has never been previously addressed in denoising, synthesis or measuring processes of fractal signals. Our study is complementary to other studies based mostly on a discrete formulation. The general framework and the simplicity of our formulation are its main assets (Fig. 1). Before explaining the details of our approach, we summarize the definitions used related to the length measurement of fractal stochastic processes. 1.2. Definition The Euclidian length of a deterministic signal x(t) for a limited duration T can be defined as follows: L0 ¼ T l ,
ð1Þ
Fractal signal (H=0.5) and filtered fractal signal 1
amplitude
0.5 0 −0.5 −1 −1.5 0
200
400
0
200
400
600
800
1000
600
800
1000
1
Amplitude
0.5 0 −0.5 −1 −1.5 sample Fig. 1. Low pass filtering increases the regularity of the fractal signal. (a) Fractal signal (H= 0.5). (b) Filtered fractal signal.
ARTICLE IN PRESS J.-M. Girault et al. / Signal Processing 90 (2010) 2690–2697
where Z 1 T=2 l¼ lðtÞ dt T T=2
0.2 0
ð2Þ
The segment λ is variable and the projector τ is constant.
−0.2
is the arithmetic mean of the elementary length and where sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 dxðtÞ lðtÞ ¼ 1 þ : ð3Þ dt On the hypothesis that x(t) is homogeneous to Volt, then l is homogeneous to Volt per second. The elementary length l can be rewritten according to the Pythagorean theorem as qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 ð4Þ lðtÞ ¼ lim t2 þ Dðt, tÞ2 , t-0 t
amplitude
2692
−0.4 −0.6 λ
−0.8 −1 −1.2
τ
−1.4 30
35
40
45
50
time Fig. 2. Principle of length measurement of a stochastic signal.
where the increment Dðt, tÞ is written in a symmetric form as
Dðt, tÞ ¼ xðt þ t=2Þxðtt=2Þ:
ð5Þ
Finally, the signal length is written as Z T=2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 Lx ¼ lim t2 þ Dðt, tÞ2 dt:
The definition of Eq. (8) extends easily to stationary and ergodic stochastic signals by replacing the arithmetic mean by the mathematical expectation E[ ]. The definition then becomes
ð6Þ
Lx ðlÞ ¼ TE½l,
t-0 T=2
t
Note that we had already used such a definition for deterministic signals (see [13]). Unlike deterministic signals, which are almost always derivable and have a finite length, stochastic signals do not follow the laws of Euclidian geometry. Stochastic signals cannot be described by continuous mathematical curves and their irregularities are still present even when examined with higher and higher resolutions. This leads to the paradoxical conclusion that the lengths of such signals are not finite. To overcome the fact that (as fractal signals) stochastic signals do not have a finite length, Richardson [31] assumed that the length Lx of a fractal signal depends on the resolution at which it is measured, i.e. the gauge l of the measurement. To overcome this indetermination for l-0, Richardson proposed measuring the curve length for each value l, 8l 40. Mandelbrot demonstrated that measuring Lx ðlÞ as a function of l was one way to determine the fractal dimension: 1Dx
Lx ðlÞ ¼ Kx l
,
ð7Þ
where Dx is the fractal dimension and Kx is a constant. Length definition becomes a function of l and is written like Eq. (1): Lx ðlÞ ¼ T lðtÞ,
ð8Þ
where
lðtÞ ¼
1 T
Z
T=2
lðt, tÞ dt,
ð9Þ
T=2
and where the ffi elementary length is qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi lðt, tÞ ¼ ð1=tÞ t2 þ Dðt, tÞ2 . The corresponding length estimator of Eq. (8) illustrated in Fig. 2 is Nk qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 X ðkT s Þ2 þ Dðn,kT s Þ2 , ð10Þ L^ x ðkT s Þ ¼ kT s n ¼ 1 where Ts is the sampling period and N is the number of elementary lengths between ½T=2,T=2.
ð11Þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 where l ¼ ð1=tÞ t2 þ D . Finally, the definition becomes qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi T ð12Þ t2 þ D2 : Lx ðtÞ ¼ E
t
The calculation of the mathematical expectation of Eq. (12) does not lead to an analytical formulation, and only a numerical resolution is possible. In order to present an analytical formulation, we propose using another definition based on that proposed by Burlaga and Klein [2]: Lx ðtÞ ¼
T
t
E½jDj:
ð13Þ
This definition supposes that t 5 D, which is often true, particularly for a high sampling frequency. However, if this is not the case, an amplitude gain can be applied. A discussion of the sense of such a definition is reported in [21]. The length estimator that we finally used which corresponded to a discrete version of Eq. (13) was Lx ðkÞ ¼
Nk X
jDðn,kÞj:
n¼1
The performance of such estimator in its average version, known as the Burlaga and Klein [2] estimator LBK(k), was excellent and outperformed the spectral estimators and the maximum likelihood estimators (see [11]). However, the performance was similar to that obtained with estimators based on wavelet analysis (see [11]). Note that Flandrin [10] demonstrated the link between length offfiq the wavelet measurement LBK(k) and variance ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi coeffipffiffiffiffiffiffiffiffi cient dj(k): E½LBK ð2j1 Þ ¼ 235j=2 2=p varðdj ðkÞÞ. Though this estimator, which is based on wavelet analysis, gave an excellent performance, it did not lead to a simple formulation in a continuous framework. The strong point of this estimator based on wavelet analysis was its adaptability to study multi-fractal processes.
ARTICLE IN PRESS J.-M. Girault et al. / Signal Processing 90 (2010) 2690–2697
The paper is organized as follows. Section 2 presents the approach in detail, setting out the notion of length ratio of colored processes, then stating the link between the length ratio and the fractal dimension of colored processes. Section 3 presents several examples of filters and the simulation results are discussed. Finally we propose a conclusion and some future applications.
t
E½jx1 x2 j,
2 T sx qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Lx ðtÞ ¼ pffiffiffiffi 1rx ðtÞ:
p t
ð20Þ
Similarly the length of the output signal y(t) is written as 2 T sy qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Ly ðtÞ ¼ pffiffiffiffi 1ry ðtÞ: ð21Þ
p t
Lx ðtÞ ¼
The method that we propose is new and was investigated in a general stochastic framework with an explicit formulation in continuous time. This approach requires some a priori information such as the probability density of the input process and the correlation function of the filter. In a linear system where the impulse response is h(t), the input signal x(t) and the output signal y(t), and the length Lx ðtÞ of the input signal, can be evaluated by T
and finally
Note that for a white noise ðrx ¼ 0Þ, the length of the input signal x(t) is written as
2. Length ratio and fractal dimension of filtered stochastic signals
Lx ðtÞ ¼
2693
ð14Þ
T sx
t
g,
ð22Þ
pffiffiffiffi where the factor g is equal to 2= pffiffiffi p for a Gaussian noise and the factor g is equal to 2= 3 for a noise uniformly distributed pffiffiffiffiffiffi between ½a=2,a=2 with a standard deviation sx ¼ a= 12. The factor g (and hence the length) seems to reveal the statistic nature of the white noise. In order to establish a link between the impulse response or the correlation function of the filter ry ðtÞ and the fractal dimension Dy of the filtered signal y(t), it should be remembered that in the case of fractal signals, the length Ly of the output signal is related to the fractal dimension Dy as follows (see also Eq. (7) for the input): Ly ðtÞ ¼ Ky t1Dy ,
where x1 ¼ xðt1 þ t=2Þ, x2 ¼ xðt2 t=2Þ 8t1 ,t2 and where Z Z jx1 x2 jpx1 ,x2 ðx1 ,x2 Þ dx1 dx2 : ð15Þ E½jx1 x2 j ¼
where Ky is a constant. Consequently, the fractal dimension Dy of the signal y(t) can be defined as follows:
px1 ,x2 ðx1 ,x2 Þ is the joint probability density.
Dy ¼ 1lim
In the special case where x1 and x2 are two random Gaussian variables, the joint probability density px1 ,x2 ðx1 ,x2 Þ is written as ð1r2x Þ1=2 ð1=ð2ð1r2x ÞÞÞðx2 =s2 þ x2 =s2 2rx x1 x2 =s1 s2 Þ 1 1 2 2 e , 2ps1 s2 ð16Þ
where s21 and s22 are the variances of the two random variables x1 and x2, respectively. rx is the correlation coefficient defined by rx ¼ E½x1 x2 =E½x21 E½x22 . The probability density of the increment D ¼ x2 x1 is written as Z þ1 px1 x2 ðx1 ,x1 þ DÞ dx1 : ð17Þ pD ðDÞ ¼ 1
After calculation, the probability density of the increment is reduced to a normal law Nð0, s0 2 Þ: 2 1 02 ð18Þ pD ðDÞ ¼ pffiffiffiffiffiffi eD =2s , 2ps0 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi with s0 ¼ s21 þ s22 2rx s1 s2 . x1 and x2 belong to the ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sx and the standard same signal x(t), then s1 ¼ s2 p deviation s0 is simplified, s0 ¼ sx 2ð1rx Þ. The mathematical expectation of the increment becomes rffiffiffiffi Z þ1 2 0 E½jDj ¼ 2 DpD ðDÞ dD ¼ s, ð19Þ
0
logðLy ðtÞÞ : logðtÞ
ð24Þ
logðLx ðtÞÞ : logðtÞ
ð25Þ
Similarly:
2.1. Filtered Gaussian noise
px1 ,x2 ðx1 ,x2 Þ ¼
t-0
ð23Þ
p
Dx ¼ 1lim
t-0
By using Eqs. (20), (25) and (21), (24) the equations become Dy ¼ 2lim
1 logð1ry ðtÞÞ , logðtÞ
ð26Þ
Dx ¼ 2lim
1 logð1rx ðtÞÞ : logðtÞ
ð27Þ
t-0 2
t-0 2
For unidimensional signal, the fractal dimension of signals x(t) and y(t) was the topological dimension of the space minus a term related to the inner correlation of x(t) and y(t). In order to simply qualify the influence of linear filters, we propose evaluating a new parameter, i.e. the length ratio O between the input and the output of the linear system:
OðtÞ ¼
Ly ðtÞ : Lx ðtÞ
By using Eqs. (20) and (21), Eq. (28) becomes sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sy 1ry ðtÞ : OðtÞ ¼ sx 1rx ðtÞ
ð28Þ
ð29Þ
Furthermore, the following function LðtÞ is defined as
LðtÞ ¼ lim
t-0
logðOðtÞÞ : logðtÞ
ð30Þ
LðtÞ is an important parameter since it quantifies the filtering impact in terms of reduction irregularity. By
ARTICLE IN PRESS 2694
J.-M. Girault et al. / Signal Processing 90 (2010) 2690–2697
using Eqs. (7), (23) and (28), it is easy to show that
L ¼ Dx Dy :
2.3. Impact of the correlation coefficient ð31Þ
In this case Dy is defined by Dy ¼ Dx lim
t-0
logðOðtÞÞ : logðtÞ
ð32Þ
When the input signal is a Gaussian noise and when the system is linear, by using Eqs. (29) and (32), this becomes 1 logð1ry ðtÞÞ 1 logð1rx ðtÞÞ þ lim : Dy ¼ Dx lim t-0 2 t-0 2 logðtÞ logðtÞ ð33Þ The fractal dimension of the filtered signal is then equal to the fractal dimension of the signal before filtering minus a term related to the filter plus a term explaining the inner correlation of the input signal x(t). Before testing different types of filter, we present the filtered fractional Brownian motion.
From definitions (33) and (37), it becomes obvious that the asymptotic behavior of the correlation coefficient when t-0 is a key point in the estimation of fractal dimension. Indeed, by exciting filters by a white Gaussian noise, all correlation coefficients ry ðtÞ whose Taylor expansion satisfies the equation limr ðtÞ 1kt t-0 y
n
,
ð38Þ
where k is a known constant and where n 2 R and 0 rn r 2, involve a fractal dimension and a Hurst exponent described as follows: 8 n > < D y ¼ Dx , 2 ð39Þ > H ¼ H þ n, : y x 2 where Dx = 2 and Hx =0 for rx ¼ 0. Note that for 2 Z Dy Z 1, n satisfies 0 r n r2. Moreover, by using definition (33) and Eq. (38) for rx ¼ 0, it becomes 1 logðktn Þ n ¼ 2 , t-0 2 logðtÞ 2
2.2. Filtered fractional Brownian motion
Dy ¼ 2lim
Let us consider the fractional Brownian motion (fBm) [25] as being in the input of the filter. This self-similar process which generalizes the Brownian motion or the Wiener process is composed of stationary and independent increments. These increments follow a Gaussian law1: x1 x2 Nð0,2t2Hx Þ where Hx is the Hurst exponent. Note that an average power spectrum density of fBm can be defined in the frequency domain [9] and it turns out that it is proportional to 1=f a with a ¼ 2Hx þ1. By using Eqs. (14)–(20), we show that
and n= 4 2Dy. Without making out an exhaustive list, we propose the following correlation coefficients (or a combination) satisfying Eq. (38) for t 4 0:
2 Lx ðtÞ ¼ pffiffiffiffi T tHx 1 :
p
ð40Þ
ry ðtÞ ¼ a þ ð1aÞcosðktn=2 Þ. For a ¼ 0:5 and for n= 2 we recognize the Hanning function;
ry ðtÞ ¼ ð1ktn Þ. For n= 1 we recognize the triangle function; n
ry ðtÞ ¼ ekt . For n =2 we recognize the Gaussian function and for n = 1 we recognize the Lorentzian function.
ð34Þ
pffiffiffiffi Note that we recognize the factor 2= p which underlines the Gaussian nature of the process. By using Eqs. (7)–(25), it becomes
3. Examples
Dx ¼ 2Hx ,
The impulse response h(t) of a Lorentzian filter is defined by
ð35Þ
where logðLx ðtÞÞ : Hx ¼ 1þ lim t-0 logðtÞ
3.1. Lorentzian filters
hðtÞ ¼ oc eoc t cosðo0 tÞuðtÞ, ð36Þ
Note that Hy has the same definition as Hx. As for Eq. (33), it becomes 1 logð1ry ðtÞÞ : ð37Þ Hy ¼ Hx lim t-0 2 logðtÞ Similarly, the Hurst exponent Hy of the filter output is the Hurst exponent Hx before filtering minus a term related to the inner filter properties. 1 Note that the increments follow a Gaussian law written as follows: x1 x2 Nð0, t2Hx Þ. However, in order to make a direct link between equations in this subsection and equations in Section 2.1, we voluntary modified this formula by doubling the variance and this modification does not modify the definition of fBm.
ð41Þ
where u(t) is the Heaviside function, oc the cut-off pulsation and o0 the central pulsation of the filter. This linear system of impulse response h(t) is excited by a Gaussian white noise. For linear filters, the crosscorrelation function of the output signal is written as Ry ðtÞ ¼ Rh ðtÞ Rx ðtÞ,
ð42Þ
where is the convolution operator. Note that the process coming from this type of filter is also known as being the ‘‘first order Markov process’’. For the Lorentzian filter the correlation function is Rh ðtÞ ¼ ðoc =4Þeoc jtj cosðo0 tÞ, the correlation function of the input signal is Rx ðtÞ ¼ s2x dðtÞ, and the correlation function of the output signal reduces to Ry ðtÞ ¼ s2y eoc jtj cosðo0 tÞ, where pffiffiffiffiffiffi sy ¼ sx oc =2. The correlation coefficient defined by ry ðtÞ ¼ Ry ðtÞ=Ry ð0Þ, is reported in Table 1 for t 40.
ARTICLE IN PRESS J.-M. Girault et al. / Signal Processing 90 (2010) 2690–2697
The output signal is a white noise, i.e. an irregular signal which fills all the space.
Table 1 Correlation coefficients and length ratios 8t 4 0. Filters
h(t)
ry ðtÞ
Wide band oc eoc t cosðo0 tÞuðtÞ eoc t cosðo0 tÞ cosðo0 tÞuðtÞ
Narrow band Base band All pass
cosðo0 tÞ
oc eoc t uðtÞ
eoc t
dðtÞ
0
OðtÞ=G pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1eoc t cosðo0 tÞ o tpffiffiffi 0 sin 2 2 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1eoc t 1
When the input is a fractional Brownian motion, the Hurst exponent of the output signal is:
for a base band filter o0 -0 and for t-0, oc t 5 1: 1 logðoc tÞ : t-0 2 logðtÞ
Hy ¼ Hx þ lim
Where rx ¼ 0 8t 40 and where the gain is G ¼ sy =sx .
As already mentioned, the length of the white Gaussian noise is 8t 4 0: 2T sx Lx ðtÞ ¼ pffiffiffiffi t,
p
ð43Þ
2T sy pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Ly ðtÞ ¼ pffiffiffiffi t 1eoc t cosðo0 tÞ:
ð44Þ
The length ratio 8t 40 is written as follows: pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi OðtÞ ¼ G 1eoc t cosðo0 tÞ,
ð45Þ
p
where the gain G is written as G ¼ sy =sx . In order to test different kinds of Lorentzian filter, four filters were considered: a wide band filter, a narrow band filter, a base band filter and an all pass filter. The correlation coefficients and length ratios for these different kinds of filter are summarized in Table 1. We verified that for these different kinds of filter, the parameter O=G was related to the intrinsic filter parameters, oc and o0 . When the input is a white Gaussian noise, the fractal dimension Dy of the output is 1 logð1eoc t cosðo0 tÞÞ : ð46Þ Dy ¼ 2lim t-0 2 logðtÞ
By imposing Hx = 0.3 (high irregularity), the Hurst exponent of the output is Hy ¼ 0:3 þ0:5 ¼ 0:8. This confirms the fact that the greater the filtering the greater the regularity; for a narrow band filter oc -0 and for t-0, o0 t 5 1: 1 logðo20 t2 =2Þ : logðtÞ t-0 2
since the correlation coefficient is rx ¼ 0, 8t 4 0. The length of the output signal is 8t 4 0:
Hy ¼ Hx þ lim
By imposing Hx = 0.3 (high irregularity), the Hurst exponent of the output is Hy ¼ 0:3 þ 1 ¼ 1:3. As 0 r Hy r 1, this means that Hy =1. 3.2. Simulations As an illustration, Figs. 3 and 4 show the length ratios for the four types of Lorentzian filter. Time series of 10 000 samples are simulated with a sampling frequency of 10 000 Hz for different values of the cut-off frequency fc ¼ f1,10,100,1000,5000g Hz and for different values of the central frequency of the filter f0 ¼ f1,10g Hz. The minimum and maximum values of the gauge are tmin ¼ 104 s and tmax ¼ 0:5 s. The sampling frequency is thus an important setting parameter since it fixes the minimum resolution. Fig. 3 shows that for a cut-off frequency fc =5000 Hz and a central frequency f0 =1 Hz (similar to an all pass
Length ratio of a low pass filtered white gaussian noise (f0=1 Hz)
For a base band filter o0 -0 and for t-0, oc t 51: Dy ¼ 212 ¼ 32,
2695
ð47Þ 100
with
1/ω c
D=2
limð1eoc t Þ oc t:
t-0
For a narrow band filter oc -0 and for t-0, o0 t 5 1: ð48Þ
with limð1cosðo0 tÞÞ o20 t2 =2:
t-0
The output signal is a sinusoid, i.e. a continuous signal whose fractal dimension is the same as a line. Note that for a Gaussian correlation coefficient 2 ry ðtÞ ¼ eoc t , the fractal dimension is also Dy =1 since 2 limt-0 ð1eoc t Þ oc t2 . For an all pass filter: 1 logð1Þ ¼ 20 ¼ 2: Dy ¼ 2lim t-0 2 logðtÞ
Ω/G
Dy ¼ 21 ¼ 1,
10−1 D=1.5
10−4
10−3
10−2
10−1
gauge (τ) Fig. 3. Length ratio of colored stochastic signals. First order Lorentzian low pass filter for different cut-off pulsations and for a fixed central pulsation.
ARTICLE IN PRESS 2696
J.-M. Girault et al. / Signal Processing 90 (2010) 2690–2697
filter), the length ratio O of such a signal has a slope equal to zero, involving a fractal dimension Dy =2. The output signal is a white Gaussian noise which entirely fills the space. Note that for a white noise, we verified that Oðtmin Þ=G 1. When the cut-off frequency was reduced to fc =1 Hz and f0 =1 Hz (low pass filter), the slope for the length ratio was 1/2 for a gauge lower than the cut-off gauge t o tc ¼ 1=oc , involving a fractal dimension Dy =3/2. The output signal was a standard Brownian motion characterized by its Hurst exponent H= 1/2. Beyond the cut-off gauge, the slope was equal to zero. Fig. 4 shows that for a central frequency f0 =100 Hz and for a cut-off frequency fc = 10 Hz, oscillations where the position of its first minimum corresponds to t ¼ 1=f0 appear in O. For a gauge less than the cut-off gauge tc ¼ 1=oc , the fractal dimension is between 1.5 and 1 for a narrow band filter and between 2 and 1.5 for a wide band filter. When the cut-off frequency tends towards zero, the fractal dimension tends to 1, thus involving the presence of a sinusoidal signal. When the cut-off frequency tends towards infinity (whatever the central frequency) then
Length ratio of band pass filtered white gaussian noise (f0=100 Hz) D=2
100
D=1.5
Ω/G
1/f
0
10−1 D=1
10−4
10−3
10−2
10−1
Gauge (τ) Fig. 4. Length ratio of colored stochastic signals. First order Lorentzian low pass filter for different cut-off pulsations and for a fixed central pulsation.
Synthesized fractal signals (n−"Lorentzian filter")
Length ratio of filtered stochastic signals
50 D=1.5, n=1 45
D=2, n=0 0
10 40
30 Ω/G
justapoxed signals
35
25
20
15
−1
10
10
5
D=1.5, n=1 D=2, n=0
0 0
1000
2000 3000 sample
4000
5000
−3
10
−2
10
−1
10
gauge (τ)
Fig. 5. Length ratio of colored stochastic signals. (a) Synthesized fractal signals with an n-order ‘‘Lorentzian filter’’. (b) Log–log plot of length ratio O=G versus gauge t.
ARTICLE IN PRESS J.-M. Girault et al. / Signal Processing 90 (2010) 2690–2697
the fractal dimension tends to 2, involving the presence of a white noise. n By using an impulse response of the form hðtÞ ¼ eoc t with 0 rn r 1, time series of 5000 samples are simulated with a sampling frequency of 10 000 Hz and a cut-off frequency fc =1 Hz. The input signal is a white Gaussian noise N(0,1). The minimum and maximum values of the gauge are tmin ¼ 104 s and tmax ¼ 0:5 s. Simulation results are depicted in Fig. 5. Fig. 5a shows synthesized fractal signals: the higher the order n, the greater the regularity of the time series. Fig. 5b shows the length ratio O=G with respect to the gauge t: the higher the order n, the higher the slope and the lower the fractal dimension. For n= 0, the corresponding filter is an all pass filter, the output signal is a white Gaussian noise and the corresponding fractal dimension is Dy = 2. For n = 1, the corresponding filter is a base band filter, the output signal is a Brownian motion and the corresponding fractal dimension is Dy = 1.5 (Hy = 0.5). For 0 r n r1 the corresponding fractal dimensions of the output filter are between 2 ZDy Z 1:5. 4. Discussion and conclusion We propose here a simple explicit formulation of the filtering effect on the regularity of a given stochastic process in a general framework and in a continuous time context. The simulations studied validated the theoretical calculus but also permitted extraction of the filter characteristics such as cut-off frequency and central frequency. Though limited to mono-fractal processes, this approach is interesting since it requires only commonly used a priori knowledge such as the probability law of the process and the spectral characteristic of the filter. The need for a priori information such as the probability density and the correlation coefficient of the filter constitute both an advantage and a drawback. It is an advantage since they provide a simple formulation, but it is also a disadvantage since the study presented here was limited to this information. Another drawback is that our approach is currently limited to mono-fractal processes. The next step will be to extend our framework to the study of multi-fractal processes. References [1] M.V. Berry, Diffractals, J. Phys. A Math. Gen. 12 (1979) 781. [2] L. Burlaga, L. Klein, Fractal structure of the interplanetary magnetic field, J. Geophys. Res. 91 (1986) 347–350. [3] J.A. Barnes, D.W. Allan, A statistical model of flicker noise, Proc. IEEE 54 (1966) 176–185. [4] J.A. Barnes, S. Jarvis, Efficient numerical and analog modeling of flicker noise processes, National Bureau of Standard Technical Note 604, 1971. [5] B.S. Chen, G.W. Lin, Multiscale Wiener filter for the restoration of fractal signals: wavelet filter bank approach, IEEE Trans. Signal Process. 42 (1994) 2972–2982. [6] B.S. Chen, W.S. Hou, Deconvolution filter design for fractal signal transmission systems: a multiscale Kalman filter bank approach, IEEE Trans. Signal Process. 45 (1997) 1359–1364.
2697
[7] B.S. Chen, B.K. Lee, S.C. Peng, Maximum likelihood parameter estimation of F-ARIMA processes using the genetic algorithm in the frequency domain, IEEE Trans. Signal Process. 50 (2002) 2208–2220. [8] M. Deriche, A.H. Tewfik, Signal modeling with filtered discrete fractional noise processes, IEEE Trans. Signal Process. 41 (1993) 2839–2849. [9] P. Flandrin, On the spectrum of fractional Brownian motions, IEEE Trans. Inform. Theory 35 (1989) 197–199. [10] P. Flandrin, Wavelet analysis and synthesis of fractional Brownian motion, IEEE Trans. Inform. Theory 38 (1992) 910–917. [11] N. Gache, P. Flandrin, D. Garreau, Fractal dimension estimators for fractional Brownian motions, in: Proceedings of the ICASSP, 1991, pp. 3557–3560. [12] J. Gagnepain, C. Roques-Carmes, Fractal approach to two-dimensional and three-dimensional surface roughness, in: International Conference on Metrology and Properties of Engineering Surfaces, vol. 109, 1986, pp. 19–126. [13] J.M. Girault, D. Kouame´, Length and frequency of band-limited signals, IEEE Signal Process. Lett. 9 (2002) 371–374. [14] M. Guglielmi, 1=f a signal synthesis with precision control, Signal Processing 86 (2006) 2548–2553. [15] G.A. Hirchoren, C.E. D’Attellis, On the optimal number of scales in estimation of fractal signals using wavelets and filter banks, Signal Processing 63 (1997) 55–63. [16] J.R.M. Hosking, Fractional differencing, Biometrika 68 (1981) 165–176. [17] J. James, Fracture surface examination of dental ceramics using fractal analysis, Dental Mater. 21 (2005) 586–589. [18] M.S. Keshner, 1/f noise, Proc. IEEE 70 (1982) 211–219. [19] N.J. Kasdin, Discrete simulation of colored noise and stochastic process and 1/f power law noise generation, Proc. IEEE 83 (1995) 802–826. [20] W.-L. Lee, K.-S. Hsieh, A robust algorithm for the fractal dimension of images and its applications to the classification of natural images and ultrasonic liver images, Signal Processing 9 (6) (2010) 1894–1904. [21] E. Le Tavernier, La me´thode de Higuchi pour la dimension fractal, Signal Processing 65 (1998) 115–128. [22] T. Loussot, R. Harbar, G. Jacquet, C.L. Benhamou, E. Lespesailles, A. Julien, An oriented fractal analysis for the characterization of texture: application to bone radiographs, Eusipco Signal Processing 1 (1996) 371–374. [23] T. Lundahl, W.J. Ohley, S.M. Kay, R. Siffert, Fractional Brownian motion: a maximum likelihood estimator and its application to image texture, IEEE Trans. Med. Imaging 5 (1986) 151–161. [24] A. Majumdar, B. Bhushan, Fractal model of elastic–plastic contact between rough surfaces, J. Tribol. 113 (1991) 1–11. [25] B.B. Mandelbrot, J.W. Van Ness, Fractional Brownian motions, fractional noises and applications, SIAM Rev. 10 (1968) 422–437. [26] B.B. Mandelbrot, A fast fractional Gaussian noise generator, Water Resource Res. 7 (1971) 543–553. [27] B. Mandelbrot, The Fractal Geometry of Nature, Freeman, San Francisco, 1982. [28] B.B. Mandelbrot, D.E. Passoja, A.J. Paulay, Fractal character of fracture surfaces of metals, Nature 308 (1984) 721–722. [29] B. Mandelbrot, Self-affine fractals and fractal dimension, Phys. Scr. 32 (1986) 257–260. [30] P. Qingrui, Y. Zhuangzhi, Mammogram analysis based on wavelets and fractals, in: Proceedings of IEEE Engineering in Medicine and Biology Society, 1996, pp. 1159–1162. [31] L. Richardson, The problem of contiguity: an appendix of statistic of deadly quarrels, General Systems Yearbook 61 (1961) 139–187. [32] R. Saletti, A comparison between two methods to generate 1/f noise, Proc. IEEE 74 (1986) 1595–1600. [33] N. Sarkar, B.B. Chaudhuri, An efficient differential box-counting approach to compute fractal dimension of image, IEEE Trans. Syst. Man Cybernet. 24 (1994) 115–120. [34] G.W. Wornell, A.V. Oppenheim, Estimation of fractal signals from noisy measurements using wavelets, IEEE Trans. Signal Process. 40 (1992) 611–622. [35] J.J. Wu, Characterization of fractal surfaces, Wear 239 (2000) 36–47.