Anas Husseis
email: anas.husseis@telecom-paristech.fr

Anissa Mokraoui
email: anissa.mokraoui@univ-paris13.fr

Basarab Matei
email: matei@lipn.univ-paris13.fr

L2ti

REVISITED HISTOGRAM EQUALIZATION AS HDR IMAGES TONE MAPPING OPERATORS

High Dynamic Range (HDR) image provides higher perceptual quality such that it appears considerably more realistic and attractive for the human observer. Since most of current screens are Low Dynamic Range (LDR) screens, lots of researches have been proposed to design tone mapping algorithms converting the HDR images into a range that is suitable to display these tone mapped images on standard LDR screens. For this purpose, this paper first investigates the pixels distribution in the HDR image through the study of various traditional Histogram Equalization (HE) algorithms which were originally developed to improve the visual quality of LDR images such as the contrast enhancement of the latter. Moreover, a modification of the Histogram Adjustment based Linear to Equalized Quantizer (HALEQ), developed for HDR images, is proposed. Simulation results show that the proposed modification preserves more details than the original version of the algorithm in most parts of the HDR image.

INTRODUCTION

High Dynamic Range (HDR) images consists of pixels, which contains greater range of colors and brightness levels than the range of Low Dynamic Range (LDR) images. Therefore displaying the HDR images on standard LDR screens remains a challenging problem. Indeed reducing the high into low dynamic range relies on a lossy transformation. This latter is known as Tone Mapping (TM) operation. Recently, significant number of Tone Mapping Operators (TMOs) converting the HDR images to LDR images have been proposed using different strategies. Among the developed methods, some of them were designed with respect to the Human Visual System (HVS) [START_REF] Tumblin | Two Methods for Display of High Contrast Images[END_REF], another were based on the frequency domain [START_REF] Tumblin | LCIS : A Boundary Hierarchy for Detail Preserving Contrast Reduction[END_REF]- [START_REF] Thai | Performance Evaluation of High Dynamic Range Image Tone Mapping Operators Based on Separable Non-linear Multiresolution Families[END_REF], segmentation operators were also designed to combine multiple functions to different segments [START_REF] Lischinski | Interactive Local Adjustment of Tonal Values[END_REF], [START_REF] Krawczyk | Lightness Perception in Tone Reproduction for High Dynamic Range Images[END_REF]. However, tone mapped images suffer from data loss because of the reduction process into low dynamic range, affecting the quality of the tone mapped images.

This paper focuses on tone mapping operators that consider both Histogram Equalization (HE) and human sensitivity to the light function. To do so, a logarithmic transform is first applied on the HDR luminance to produce a distribution that is close to the HVS [START_REF] Drago | Adaptive logarithmic mapping for displaying high contrast scenes[END_REF], [START_REF] Stockham | Image processing in the context of a visual model[END_REF]. The target then is to redistribute this logarithmic luminance on a given number of intervals corresponding to the standard LDR display such that it is equally distributed on each interval. This ensures the randomness maximization on the redistributed probability density function (pdf) respecting the human sensitivity to the light.

Several well known HE techniques have been deployed to enhance the LDR image contrast by redistributing the pdf of the pixels to obtain a flat histogram. This paper first investigates these traditional HE methods to render the HDR images at LDR displays to not only preserve the contrast and brightness but also achieve maximum randomness in the tone mapped images. Then, a modification is made to the Histogram Ajustement Linear to Equalization algorithm (HALEQ) dedicated for tonemapped HDR images [START_REF] Jiang Duan | Tone-mapping high dynamic range images by novel histogram adjustment[END_REF].

The rest of this paper is organized as follows. Next section revisits the traditional histogram equalization methods, originally developed to improve the visual quality of LDR images, to be considered as TMO for HDR images. Section 3 presents the generic Histogram Adjustment based Linear to Equalized Quantizer (HALEQ) TMO. A modification is then made to improve its performance in terms of objective and subjective quality. Section 4 analyses and discusses the simulation results. Section 5 concludes this paper.

REVISITED HISTOGRAM EQUALIZATION AS HDR IMAGE TMO

This section concerns HE methods to map HDR into LDR images while preserving the most important details in the HDR images in accordance to the range of the LDR device. These methods are derived from those initially deployed for LDR images to enhance their visual quality. These techniques, considered in this paper as TMOs, aim to produce images with an applicable dynamic range on the standard LDR display. Note that the probability density function of the HDR pixel values usually behaves in Laplacian manner and manipulating the linear values of luminance causes heavy distortion because of the insignificant relation to the visual perception [START_REF] Boschetti | High Dynamic Range image tone mapping based on local Histogram Equalization[END_REF]. Weber law can be used to conclude a straightforward relationship between the sensitivity function and the logarithmic transform of the luminance. Therefore to take into account the human sensitivity function to light, the derived HE techniques will operate on the logarithmic luminance of the HDR images.

A HE technique is a nonlinear operation that modifies, in an appropriate way, the initial distribution (i.e. pdf) of the pixels in the HDR image to produce a uniform distribution (i.e. a flat histogram) over the range of the LDR display. The mathematical explanation behind this strategy is based on the entropy of the source. Consider a source of logarithmic luminance events {x 1 , . . . , x n } with their probabilities {p(x 1 ), . . . , p(x n )} (i.e. the source pdf). The entropy, denoted H(x), is given by the average amount of information as follows [START_REF] Shannon | A mathematical theory of communication[END_REF] :

H(x) = - n i=1 p(x i )log(p(x i )).
(

) 1
According to this definition, the entropy is maximized when the events (logarithmic luminance levels) are uniformly distributed. Therefore applying HE operation to the source ensures that maximum possible amount of information is contained in the tone-mapped image. Therefore this paper exploits HE techniques to map the logarithmic luminance from HDR into LDR such that the reduced range corresponds to the LDR display, taking into account the pixels population distribution. In the following, four HE methods are investigated : starting with Global Histogram Equalization (GHE), continuing with HE development in order to maximize the brightness preservation using Brightness Preserving Bi-Histogram Equalization (BBHE) [START_REF] Kim | Contrast enhancement using brightness preserving Bi-Histogram equalization[END_REF], Dualistic Sub-Image Histogram Equalization (DSIHE) [START_REF] Wang | Image enhancement based on equal area dualistic sub-image histogram equalization method[END_REF], and Recursive Mean-Separate Histogram Equalization (RMSHE) [START_REF] Chen | Contrast enhancement using recursive mean separate histogram equalization for scalable brightness preservation[END_REF].

Before presenting these techniques, let us introduce some notations. Assume that the HDR image is of size N = n × m pixels. The HE divides the logarithmic luminance HDR into K intervals in accordance to the bit depth of the LDR screen, such that each interval contains n×m K pixels. Denote L w (i, j) the logarithmic luminance of the pixel at position (i, j) on the HDR image. Denote L wmin (respectively L wmax ) the smallest (respectively highest) logarithmic luminance of the HDR image.

Global Histogram Equalization

Global Histogram Equalization (GHE) maps the logarithmic luminance HDR into LDR of K levels such that all reproduced levels have equal amount of pixels. GHE objective is to produce a luminance that takes values in the sub-intervals

[D k , D k+1 ] such that L d (i, j) ∈ {D 0 , D 1 , . . . , D K-1 }.
The probability density function (pdf) according to the previous model is formulated as follows :

p(D k ) = N k N , (2)
where D k is the display level, N k is the pixel population at each interval [D k , D k+1 ] and N is the total amount of pixels in the image. The cumulative distribution function (cdf) of D k is given by :

C(D k ) = k j=0 p(D k ) = k j=0 N k N , k ∈ [0, K -1]. (3)
C(D k ) is a non-decreasing function because it is the empirical probability distribution function and it takes the maximum value at the last point such that C(D K-1 ) = 1.

To equalize and map the histogram, primarily the discrete cumulative normalized histogram C(D k ) has to be calculated and applied at each logarithmic luminance value according to the following formula :

L d (i, j) = f (L w (i, j)) = D 0 + (D K-1 -D 0 )C(L w (i, j)), (4)
where L w (i, j) is the logarithmic luminance value and C(L w (i, j)) is its cumulative probability which considers K intervals corresponding to the LDR screen. The display luminance L d has approximately a flat histogram not a perfectly flatted one, because histogram bins are never decreased in amplitude, though they may grow if various luminance levels map to the same value.

GHE produces L d logarithmic luminance values which are suitable with LDR displays, taking into account pixel population distribution. The main drawback of GHE is that the mean value of L d does not take into account the mean value of L w . In other words, the brightness of the HDR image is not preserved in the tone mapped image. Several techniques were then proposed to bypass this drawback as it will be discussed in the next sections.

Brightness Preserving Bi-Histogram Equalization

Brightness Preserving Bi-Histogram Equalization (BBHE) technique is an intuitive technique, originally used to preserve the brightness of the enhanced LDR images [START_REF] Kim | Contrast enhancement using brightness preserving Bi-Histogram equalization[END_REF]. This approach is exploited as a tone mapping operator. It consists of cutting the logarithmic luminance HDR histogram into two sub-intervals using a selected threshold L T . This latter is deduced as the average of the logarithmic luminance values ([L wmin , L wmax ]) of the HDR image. The GHE strategy, described in section 2.1, is then applied on each sub-interval (i.e. [L w,min , L T ] and [L T , L w,max ]). This separation into two sub-intervals grants the preservation of the mean brightness of the HDR image.

Dualistic Sub-Image Histogram Equalization

Dualistic Sub-Image Histogram Equalization (DSHE) technique has been developed to enhance the visual quality of LDR images [START_REF] Wang | Image enhancement based on equal area dualistic sub-image histogram equalization method[END_REF]. In this section, it is exploited as a TMO for HDR images. The original logarithmic luminance of the HDR image is decomposed into two sub-intervals where each of them contains the half amount of the total number of pixels of the HDR image. This proposition is based on Shannon entropy function :

h = -plog(p) -(1 -p)log(1 -p), (5)
where p corresponds to the probability distribution of the first interval. The entropy of the aggregation dualistic segmentation is maximized when p = 0.5, consequently (1 -p) represents the probability distribution of the other interval. In order to satisfy this assumption, L T is defined as the median of L w . Then each interval is equalized using the GHE strategy (see section 2.1). In this technique, significant shifts on the average logarithmic luminance is avoided especially for the large regions with same luminance level.

Recursive Mean-Separate Histogram Equalization

Recursive Mean-Separate Histogram Equalization (RMSHE) technique has been developed for better brightness preserving of the LDR images [START_REF] Chen | Contrast enhancement using recursive mean separate histogram equalization for scalable brightness preservation[END_REF]. It is exploited in this paper as a TMO for HDR images. Instead of applying a single cut, as in BBHE, a recursive process takes place to cut in an iterative way the sub-intervals according to different cutting points {L T1 , L T2 , . . . , L T N I-1 } (with r the number of iterations that produces N I = 2 r subintervals) deduced as in section 2.2. The GHE technique is applied at each sub-interval as discussed before. Note that choosing r = 0 is equivalent to GHE, while r = 1 corresponds to BBHE. So this method provides a scalable brightness preservation higher than that in BBHE and DSIHE while choosing r 2.

MODIFIED HISTOGRAM ADJUSTMENT BASED LINEAR TO EQUALIZED QUANTIZER

Before presenting the modification made to the Histogram Adjustment based Linear to Equalized Quantizer (HALEQ), developed in [START_REF] Jiang Duan | Tone-mapping high dynamic range images by novel histogram adjustment[END_REF], this section starts with the description of the generic algorithm.

HALEQ is a hybrid TMO that considers both HE and linear mapping (uniform quantization) while reducing the high dynamic range. The logarithmic transform is applied to the world luminance as a first step through the following expression [START_REF] Drago | Adaptive logarithmic mapping for displaying high contrast scenes[END_REF], [START_REF] Stockham | Image processing in the context of a visual model[END_REF] :

D(I) = (D max -D min ) × log(I+τ )-log(Imin+τ ) log(Imax+τ )-log(Imin+τ ) + D min , (6)
where I is the world luminance ; D min and D max are the minimum and maximum level values for the LDR display device (typically they are [0,255]) ; I min and I max are extracted from the HDR image to represent the lowest and highest luminance values in the scene ; and τ is a constant which controls the overall brightness of the tone mapped image.

The tone mapping operation is performed in a recursive algorithm such that at the first iteration (i.e. r = 1) the dynamic range is divided twice : (i) it is divided into two equal sub-intervals using the cutting point l n (linear mapping) ; and (ii) it is divided into two sub-intervals, using the cutting point e n , in order to get the same number of pixels in each interval (HE). An adjustment process is then performed by choosing a new cutting point l en between these two cutting points according to the following relationship :

l en = l n + β(e n -l n ), ( 7
)
where β is a fixed constant taking its value between 0 and 1. This new cutting point l en is controlled by the parameter β which is fixed once by the authors at the beginning of the process. The described process is then iterated. Indeed at each iteration the dynamic range of the sub-intervals are cut once again to produce 2 r sub-intervals where r is the number of iterations corresponding to the capability of the LDR screen.

In order to improve, in particular, the visual quality of the tone mapped HDR image, this section proposes to adapt the value of β to each sub-intervals rather than using the same constant value for all sub-intervals (see equation ( 7)). Equation ( 7) becomes :

l en = l n + β n (e n -l n ), (8)
where β n is a parameter depending on the sub-interval

[I n , e n ].
Moreover, we propose to compute the cutting point l en as the average luminance on the sub-interval [I n , e n ]. From equation [START_REF] Wang | Image enhancement based on equal area dualistic sub-image histogram equalization method[END_REF], the β n is then deduced as follows :

β n = D(I en,ln ) -l n e n -l n , (9)
where D(I en,ln ) is the average logarithmic luminance in the sub-interval [e n , l n ].

This strategy provides an automatic, simple and efficient way to assign the value β n at each cut such that the brightness preservation is maximized in the tone mapped images.

SIMULATION RESULTS

This section evaluates the performance of the revisited and modified histogram equalization methods that have been presented in the previous sections. Simulations have been conducted under Matlab environnement using the HDR Toolbox with its HDR test images [START_REF] Banterle | Advanced High Dynamic Range Imaging : Theory and Practice[END_REF]. HE methods are tested on 10 indoor and outdoor HDR images. An objective performance evaluation based on TMQI (Tone-Mapped image Quality Index) metric, that considers the Structural SIMilarity and the Structural Naturalness between the tone mapped image and the HDR image [START_REF] Yeganeh | Objective Quality Assessment of Tone-Mapped Images[END_REF], has been performed. This metric is upper-bounded by 1.

Table 1 provides the TMQI of the tone mapped HDR images according to the revisited HE methods using 8 bits of color precision per R,G,B color channel. The TMQI metric values are close to 1 (ranging between 0.88 and 0.97). Subjective visual quality of these tone mapped images are also analyzed. Fig. 1, Fig. 2, Fig. 3 and Fig. 4 show the "Clock Building" and "Bottles" tone mapped HDR images. One can observe the excellent contrast and the preserved details on the mapped images. Indeed the fine details on the sky for "Clock Building" image and on the wall for the "Bottles" image are very well rendered. Over and above, the brightness preservation is maximized on the RMSHE as we assumed r = 2 to obtain 4 sub-intervals providing pleasing images with high TMQI values. 1. TMQI of the revisited HE methods. On the other hand, the optimized HALEQ TMO is also tested. Table 2 compares the TMQI metric according to the values of β considered as a constant or an optimized value. The visual quality of the tone mapped images are also compared when using a constant β value (e.g. β = 0.5) during all the process and an optimized value (see Fig. 5 and Fig. 6). One can observe that fine details are better rendered on the optimized β. This can be obviously noticed from the sky on "Clock Building" image and the wall on "Bottles" image. This is explained by the fact that β is well adapted to the different sub-intervals. Fig. 7 (respectively Fig. 8) shows the different values of β that have been used by the optimized HALEQ TMO for the "Clock Building" (respectively "Bottles") HDR image. The sensitivity curve given by Fig. 12, for "Clock Building", shows that the curve for β = 0.75 is close to the curve with optimal β that is why the TQMI metric are almost equal.

Fig. 9 illustrates the pdf of the luminance "Clock Building" HDR image. The normalized histograms of the equalized HDR image according respectively to the generic and modified HALEQ are provided by Fig. 10 (with β = 0.5 and a standard deviation σ = 0.0040) and Fig. 11 (with a standard deviation σ = 0.0033). The modified HALEQ histogram is flatter than that obtained by the generic HALEQ.

Fig. 12 provides the histogram of the equalized "Clock Building" HDR image with β = 0.75 and a standard deviation σ = 0.0031. The latter is equivalent to that of the modified HALEQ algorithm. That is why the TMQI metrics are similar (0.96 see Table 2). the contrast of LDR images. These HE methods have been applied on HDR images as TMO to achieve a limited quality. Moreover, a modification has been made to the generic HALEQ tone mapping operator such that the parameter β is automatically adjusted according to the pdf of the HDR image. Simultation results show that this modification preserves the brightness and details of the HDR image by optimizing the cutting points between the linear mapping and HE.

CONCLUSION

This paper revisited several histogram equalization techniques that have been essentially proposed to enhance

Fig. 1 .

 1 Fig. 1. GHE (left image) ; BBHE (right image) using "Clock Building" HDR image.

Fig. 2 .

 2 Fig. 2. DSIHE (left image) ; RMSHE (right image) using "Clock Building" HDR image.

Fig. 3 .

 3 Fig. 3. GHE (left image) ; BBHE (right image) using "Bottles" HDR image.

Fig. 4 .

 4 Fig. 4. DSIH (left image) ; RMSHE (right image) using "Bottles" HDR image.

Fig. 5 .

 5 Fig. 5. HALEQ with β=0.5 (left image) ; HALEQ with optimized β (right image) using "Bottles" HDR image.

Fig. 6 .

 6 Fig. 6. HALEQ with β=0.5 (left image) ; HALEQ with optimized β (right image) using "Clock Building" HDR image.

Fig. 7 .

 7 Fig. 7. Optimized β for "Clock Building" HDR image. HALEQ TMO β = 0.25 β = 0.5Optimized β β = 0.75 Bottles 0.870 0.915 0.916 0.920 Clock Building 0.870 0.951 0.961 0.965 Office 0.867 0.968 0.963 0.970 Apartment Backdrop 0.839 0.929 0.940 0.941 Atrium Night 0.942 0.964 0.933 0.944 Belgium 0.905 0.961 0.958 0.962 Doll 0.821 0.882 0.887 0.888 Memorial Church 0.846 0.935 0.915 0.918 Oxford Church 0.892 0.985 0.980 0.984 Small Office 0.898 0.955 0.961 0.965

Fig. 8 .

 8 Fig. 8. Optimized β for "Bottles" HDR image.

Fig. 9 .

 9 Fig. 9. "Clock Building" HDR image histogram.

Fig. 10 .

 10 Fig. 10. "Clock Building" tone mapped image histogram using the modified HALEQ (optimized β).

Fig. 11 .

 11 Fig. 11. "Clock Building" tone mapped image histogram using HALEQ with β = 0.5.

Fig. 12 .

 12 Fig. 12. "Clock Building" tone mapped image histogram using HALEQ with β = 0.75.

Fig. 13 .

 13 Fig. 13. "Clock Building" LDR versus HDR logarithmic luminance using HALEQ and modified HALEQ.

Table 2 .

 2 TMQI of HALEQ and modified HALEQ methods. dynamic range adapted to LDR screens. These techniques are not only simple and easy to implement but also provide good performance both in terms of objective and subjective

		.870	0.915	0.916	0.920
	Clock Building	0.870	0.951	0.961	0.965
	Office	0.867	0.968	0.963	0.970
	Apartment Backdrop	0.839	0.929	0.940	0.941
	Atrium Night	0.942	0.964	0.933	0.944
	Belgium	0.905	0.961	0.958	0.962
	Doll	0.821	0.882	0.887	0.888
	Memorial Church	0.846	0.935	0.915	0.918
	Oxford Church	0.892	0.985	0.980	0.984
	Small Office	0.898	0.955	0.961	0.965


Xin Wang
email: xin.wang@oit.edu

Edwin E Yaz
email: edwin.yaz@marquette.edu

Isabelle Dufour
email: dufour@ixl.fr

Susan C Schneider
email: susan.schneider@marquette.edu

MICRO-CANTILEVER BIOCHEMICAL SENSING BASED ON ROBUST FINITE-HORIZON NON-LINEAR ESTIMATION

A novel non-linear estimation based micro-cantilever biochemical rheological sensing technique is proposed in this paper. Contrary to classical rheological measurements using microcantilevers, which is either restricted to resonant phenomena, or based on the measurement of fluid properties over a range of vibration frequencies, the proposed measurement technique can be applied to arbitrary available vibration frequency. By applying non-linear estimation technique, the fluid characteristic parameters including density and viscosity, can be estimated and used to quantify fluid properties. The preliminary simulation studies demonstrate that vibrating micro-cantilevers with nonlinear estimation can be used as effective and robust fluid characterization micro-system in liquid environments. The proposed technique is encouraging for the development of a useful microrheometer on a silicon chip for fluid detection application.

INTRODUCTION

Micro-cantilevers have recently attracted much attention as promising and highly sensitive approach for minute quantities detection and have wide range of applicability in chemical, biological and medical fields. Moreover, recent technologies developed in the last few years for the fabrication technologies in micro-electromechanical systems (MEMS) and nano-electromechanical systems (NEMS) have made microcantilevers available miniaturized products. Recent literature in micro-cantilevers demonstrates that vibrating micro-cantilevers can be used to quantify fluid properties. Classical rheological measurement using micro-cantilevers is either restricted to resonant phenomena, or based on the measurement of fluid properties over a wide range of frequencies [START_REF] Belmiloud | Rheological behaviour probed by vibrating micro-cantilevers[END_REF][START_REF] Cox | Generalized model of resonant polymer-coated micro-cantilevers in viscous liquid media[END_REF][START_REF] Youssry | A straightforward determination of fluid viscosity and density using microcantilevers: from experimental data to analytical expressions[END_REF][START_REF] Dufour | Resonant microcantilever type chemical sensors: analytical modelling in view of optimization[END_REF][START_REF] Castille | Longitudinal vibration mode of piezoelectric thick-film cantilever-based sensors in liquid media[END_REF][START_REF] Belmiloud | Determination of fluid viscoelastic properties using resonant micro-cantilever[END_REF][START_REF] Heinrich | An analytical model for transient deformation of viscoelastically coated beams: applications to static-mode microcantilever chemical sensors[END_REF][START_REF] Dufour | Effect of coating viscoelasticity on quality factor and limit of detection of micro-cantilever chemical sensors[END_REF][START_REF] Dufour | Theoretical analysis of strong-axis bending mode vibrations for resonant micro-cantilever (bio)chemical sensors in gas or liquid phase[END_REF]. Contrary to the conventional methods, the proposed detection technique can be applied to any excitation vibration frequency. Based on time domain measurement at the arbitrary excitation frequency, non-linear filtering estimation techniques are used to determine the characteristic parameters of the tested liquid.

In the following section, we first derive the mechanical modelling of micro-cantilever system model. The Unscented Kalman Filter (UKF) is developed for rheological fluid properties determination. In the last section, simulation results demonstrate that the proposed non-linear estimation technique can effectively and robustly determine the fluid properties.

The complex response method is used to solve for the partial differential equation. If w = W cos(ωt) is a real-valued quantity in the time domain, let us denote w = We jωt to represent the corresponding complex notation for w, therefore, we have the property of w = Re{w} .

GENERALIZED MODEL OF RESONANT MICRO-CANTILEVERS IN VISCOUS LIQUID MEDIA

Before introducing the novel rheological fluid characterization technique with non-linear estimation, the mechanical model is derived in this section. Comparing with classical viscometers or rheometers, micro-cantilever presents some unique property such as larger frequency range measurements and compactness. These unique sensing advantages are due to its mechanical prop-erties [START_REF] Belmiloud | Rheological behaviour probed by vibrating micro-cantilevers[END_REF][START_REF] Cox | Generalized model of resonant polymer-coated micro-cantilevers in viscous liquid media[END_REF][START_REF] Rosenhead | Laminar Boundary Layers[END_REF]. The micro-cantilever rheological sensing system in liquid environment is illustrated as shown in Fig. 1. The mechanical system modelling for the micro-cantilever with the concentrated activating force applied at the free-end is derived as follows:

Generalized Equation of Motion

For transverse bending modes, a micro-cantilever sensor is immersed in fluid. The fixed end of the cantilever is denoted as the origin of the coordinate. The device is activated by the concentrated force applied at the free-end and the movement is affected by distributed force from the liquid.

For mechanics of materials, the beam sustains a bending moment M(x,t), which is related to the beam deflection w(x,t) by

M(x,t) = EI(x) ∂ 2 w(x,t) ∂ x 2 (1)
Assuming the deflection to be small enough such that the shear deflection is much smaller than w(x,t), i.e. so that the sides of the element dx do not bend, the summation of forces in y direction yields

(V (x,t) + ∂V (x,t) ∂ x dx) -V (x,t) + f (x,t)dx = ρA(x)dx ∂ 2 w(x,t) ∂t 2
(2) where V (x,t) shear force at the left end (V (x,t) + ∂V (x,t) ∂ x dx) shear force at the right end f (x,t) fluid force per unit length ρA(x)dx ∂ 2 w(x,t) ∂t 2 the inertial force of the element Notice that the assumption of small shear deformation used in the force balance of equation ( 2) is true if L/h ≥ 10, L/b ≥ 10, i.e. for long slender beam, or Euler-Bournoulli beam. h and b are the height and width of the beam, respectively.

Next, moments acting on the element dx about z axis through (Figure 2) point Q are summed, which yields

(M(x,t) + ∂ M(x,t) ∂ x dx) -M(x,t) + (V (x,t) + ∂V (x,t) ∂ x dx)dx + ( f (x,t)dx) dx 2 = 0 (3)
where f (x,t)dx is the sum of distributive force. The left hand side of the equation equals zero, since it also assumed that the rotary inertia of the element dx is negligible. Simplifying the expression yields

( ∂ M(x,t) ∂ x +V (x,t))dx + ( ∂V (x,t) ∂ x + f (x,t) 2 )dx 2 = 0 (4)
As dx → 0, dx 2 → 0 much faster. Therefore, we have

V (x,t) = - ∂ M(x,t) ∂ x (5)
This states that the shear force is proportional to the spatial change in the bending moment. Substitution of this expression for the shear force into Eqn.(2) yields

- ∂ 2 M(x,t) ∂ x 2 dx + f (x,t)dx = ρA(x)dx ∂ 2 w(x,t) ∂t 2 (6)
Further substitution of the Eqn.(1) into Eqn.( 6), and dividing by dx yields:

ρA(x) ∂ 2 w(x,t) ∂t 2 + ∂ 2 ∂ x 2 (EI(x) ∂ 2 w(x,t) ∂ x 2 ) = f (x,t) (7)
Consider the effect of the concentrated force acting at the point x = L, then Eqn.( 7) can be modified to

ρA(x) ∂ 2 w(x,t) ∂t 2 + ∂ 2 ∂ x 2 (EI(x) ∂ 2 w(x,t) ∂ x 2 ) = f (x,t) + Fδ (x -L) (8) Denote m L = ρA(x), f = f f luid , assuming EI is not a func- tion of x, then m L ∂ 2 w ∂t 2 + EI ∂ 4 w(x,t) ∂ x 4 = f f luid + Fδ (x -L) (9)
Under the incompressible and Newtonian assumption on the fluid, viscosity and density are intrinsic constants. When a micro-cantilever vibrates in a viscous fluid, the fluid offers resistance to the motion of the beam.

The fluid loading can be interpreted as the sum of two forces: an inertial force, which is proportional to the cantilever acceleration, and a viscous force, which is proportional to the cantilever velocity. Therefore, a viscous fluid exerts on the beam surface the drag force that can be expressed as

f f luid = -g 1 ( f ) ∂ w(t) ∂t -g 2 ( f ) ∂ 2 w(t) ∂t 2 (10)
From Eqn.( 9),( 10), we have

EI ∂ 4 w(x,t) ∂ x 4 + g 1 ( f ) ∂ w(t) ∂t + (g 2 ( f ) + m L ) ∂ 2 w ∂t 2 = Fδ (x -L) (11)
Hence, we have derived the system equation. Now suppose using sinusoidal excitation and adopt complex notations with time dependence e jωt , we have

EI ∂ 4 w ∂ x 4 + g 1 jωw -(g 2 + m L )ω 2 w = 0 (12)
Equivalently, we have

EI( ∂ 4 w ∂ x 4 -Λ 4 w) = 0 (13)
where

Λ 4 = m L ω 2 EI [1 + g 2 m L -j g 1 m L ω ] (14)
Let us consider the boundary conditions w| x=0 = 0 no deflection at the fixed end

∂ w ∂ x | x=0 = 0 no slope at the fixed end EI( ∂ 2 w ∂ x 2 )| x=L = 0 no bending moment at the free-end ∂ ∂ x [EI ∂ 2 w ∂ x 2 ]| x=L = F concentrated
shear force at the free-end The last two conditions are because: if a beam in transverse vibration is free at one end, the deflection and slope ∂ w ∂ x at that end are not restricted, the bending moment must be zero. Since we have force F acting at the free-end, then shear force is not zero but F. Therefore, we have the following boundary conditions: w| x=0 = 0 no deflection at the fixed end

∂ w ∂ x | x=0 = 0 no slope at the fixed end ∂ 2 w ∂ x 2 | x=L = 0
no bending moment at the free-end

∂ 3 w ∂ x 3 | x=L = F EI
concentrated shear force at the free-end

In this way, we include the effect of F by solving the boundary conditions. Solving the differential Eqn. ( 13)

EI( ∂ 4 w ∂ x 4 -Λ 4 w) = 0 ( 15
)
The characteristic equation is

λ 4 -Λ 4 = 0 (16)
Therefore, we have

λ = Λe j 2π 4 n , f or n = 0, 1, 2, 3 (17)
Therefore, we have the general solution

w = C 1 cos(Λx) +C 2 sin(Λx) +C 3 cosh(Λx) +C 4 sinh(Λx) = C 1 cos(Λx) +C 2 sin(Λx) +C 3 cos( jΛx) +C 4 (-jsin( jΛx)) (18)
Applying boundary conditions to solve for the coefficients, we have

w(x = L) = - cos (

MICRO-CANTILEVER SENSING BASED ON NON-LINEAR ESTIMATION

First, the g 1 , g 2 parameters are estimated with the Unscented Kalman Filter. Then, the important micro-cantilever sensing parameters: density and viscosity can be obtained with Maali's equation [10], which can be used to identify the fluid.

UKF Estimation for g 1 , g 2 Parameters

In this section, we set up the Unscented Kalman Filter to estimate the parameters g 1 , g 2 by the following system and measurement equations [START_REF] Wang | A New Nonlinear-Filter-Based Modulation/Demodulation Technique for Chaotic Communication[END_REF] Denote intermediate parameters

p 1 = g 2 m L p 2 = g 1 m L ω (20)
From Eqn.( 14), we have

Λ = 4 m L ω 2 EI [1 + g 2 m L -j g 1 m L ω ] = 4 m L ω 2 EI [1 + p 1 -jp 2 ] (21)
For a given cantilever geometry, for a given viscous fluid at a fixed frequency g 1 and g 2 are constant, therefore, the parameters p 1 and p 2 are constants. Therefore, we have process equation

p 1,k+1 p 2,k+1 = 1 0 0 1 × p 1,k p 2,k (22)
The measurement equations are given as follows: From Eqn.( 14), we have

Λ k = 4 m L ω 2 EI [1 + p 1,k -jp 2,k ] (23)
Denote the following notation

Ξ k = cos(Λ k L)e -Λ k L + sin(Λ k L)e -Λ k L -cos(Λ k L)e Λ k L + sin(Λ k L)e Λ k L (24)
Applying boundary conditions ( 16) to solve coefficients, we have

w k (x = L) = -Re{ Ξ k EIΛ 3 k [2 + cos(Λ k L)(e Λ k L + e -Λ k L )] Fe jωt } + v k (25)
where v k is the measurement noise, assume to be additive white Gaussian noise v k ∼ N (0,V ). Based on the system and measurement equations above, we set up the following Unscented Kalman Filter for parameters estimation.

System model: The new state vector is

χ k = p 1,k p 2,k (26)
with initial condition:

χ 0 = p 1,0 p 2,0 (27)
the total number of state variables n equals 2, κ equals 1, satisfying the condition n + κ = 3.

Propagating from time step k to (k+1): In the time update part, the 2-dimentional (n = 2) χ k with mean χ+ k and covariance P + k is approximated by (2n+1) weighted samples (sigma points) selected by (2n+1) weighted samples (sigma points) selected by

χ (0) k = χ+ k χ (i) k = χ+ k + χ(i) , i = 1, 2, ..., 2n χ(i) = ( (n + κ)P + k ), i = 1, 2, ..., n χ(i) = -( (n + κ)P + k ), i = n + 1, n + 2, ..., 2n (28)
( (n + κ)P + k )is the i th row or column of the matrix square root of (n + κ)P + k , and the (2n+1) weighting coefficients are given by

W (0)) = κ n + κ W (i)) = 1 2(n + κ) , i = 1, 2, ..., 2n (29)
Each sigma point is initiated through the process model. Since the two parameters are constants, we have

χ (i) k+1 = χ (i) k = χ (i) k (1) χ (i) k (2) = p (i) 1,k p (i) 2,k (30)
The predicted mean is computed by

χ- k+1 = 2n ∑ i=0 W (i) χ (i) k+1 = 1 n + κ {κ χ (0) k+1 + 1 2 2n ∑ i=1 χ (i) k+1 } (31)
The predicted covariance, the priori error covariance, is computed by

P - k+1 = 2n ∑ i=0 W (i) (χ (i) k+1 -χ- k+1 )(χ (i) k+1 -χ- k+1 ) T = 1 n + κ × κ((χ (0) k+1 -χ- k+1 )(χ (0) k+1 -χ- k+1 ) T )+ 1 2 ∑ 2n i=1 (χ (i) k+1 -χ- k+1 )(χ (i) k+1 -χ- k+1 ) T (32)
We implement the measurement update equations as follows. Initialize each predicted observation point using the given observation model. The state variables are the sigma points from the time update part shown above.

Λ (i) k+1 = 4 m L ω 2 EI [ 1 + p (i) 1,k+1 -jp (i) 2,k+1 ]
= 4 m L ω 2 EI [1 + χ (i) k+1 (1) -jχ (i) k+1 (2)] (33)
Therefore, we have

Ξ k+1 = cos(Λ (i) k+1 L)e -Λ (i) k+1 L + sin(Λ (i) k+1 L)e -Λ (i) k+1 L -cos(Λ (i) k+1 L)e Λ (i) k+1 L + sin(Λ (i) k+1 L)e Λ (i) k+1 L (34)
Also, we have

Y (i) k+1 = -Re{ Ξ k+1 EI(Λ (i) k+1 ) 3 [2 + cos(Λ (i) k+1 L)(e Λ (i) k+1 L + e -Λ (i) k+1 L )] Fe jωt } (35)
The predicted observation is calculated by:

ŷk+1 = 2n ∑ i=0 W (i) Y (i) k+1 = 1 n + κ {κY (0) k+1 + 1 2 2n ∑ i=1 Y (i) k+1 } (36)
Since the measurement noise is AWGN with covariance V, the covariance of the predicted measurement is:

P y = V + 2n ∑ i=0 W (i) (Y (i) k+1 -y k+1 )(Y (i) k+1 -y k+1 ) T = V + 1 n + κ × κ((Y (0) k+1 -y k+1 )(Y (0) k+1 -y k+1 ) T )+ 1 2 ∑ 2n i=1 (Y (i) k+1 -y k+1 )(Y (i) k+1 -y k+1 ) T (37)
The cross-correlation is determined by:

P xy = 2n ∑ i=0 W (i) (χ (i) k+1 -χ- k+1 )(Y (i) k+1 -y k+1 ) T = 1 n + κ × κ((χ (0) k+1 -χ- k+1 )(Y (0) k+1 -y k+1 ) T )+ 1 2 ∑ 2n i=1 (χ (i) k+1 -χ- k+1 )(Y (i) k+1 -y k+1 ) T (38)
The measurement update of the state can be performed by normal Kalman Filter equations

K k+1 = P xy P -1 y χ+ k+1 = χ- k+1 + K k+1 (y k+1 -ŷk+1 ) P + k+1 = P - k+1 -K k+1 P y K T k+1 (39)
It should be noticed that we have modified the original UKF program in Eqn.(38),(39), by replacing the original ŷ(k + 1) with y(k + 1). Simulation results show that by make this change, the performance of UKF significantly improved.

Viscosity and Density Identification

Based on the estimated results of the parameters g 1 , g 2 , the viscosity and density can be identified with Maali's equation [START_REF] Youssry | A straightforward determination of fluid viscosity and density using microcantilevers: from experimental data to analytical expressions[END_REF]:

ρ = 2π f g 2 (2a 1 b 2 -a 2 b 1 ) + a 2 ∆ π 2 b 2 f a 1 (a 1 b 2 -a 2 b 1 ) η = a 1 (2π f g 2 b 1 + ∆) 2 π(a 1 b 2 -a 2 b 1 )(2π f g 2 (2a 1 b 2 -a 2 b 1 ) + a 2 ∆) (40)
where The Maali's parameters a 1 = 1.0553, a 2 = 3.7997, b 1 = 3.8018, b 2 = 2.7364, and b is the micro-cantilever width. Therefore, based on the estimated parameters g 1 , g 2 , both of the density ρ and viscosity η can be identified.

∆ = g 1 a 2 -(2π f g 2 b 1 -a 2 g 1 ) 2 + 8π f g 2 g 1 a 1 b 2 (41)

SIMULATION RESULTS AND ANALYSIS

The following physical parameters values used in our simulation are summarized in Tab.1.

Without considering the effect of measurement noise, the simulation results are shown below:

It should be noted the percentage error results can be even further improved if more experimental measurement data is taken into account. Simulation results are given in Fig. 3 and Fig. 4, show that the free-end deflection estimation matches very well with the measured values. Based on 5 seconds of the measurement of the deflection at the free-end, the percentage error for g 1 is less than about 8 percent, while the percentage error for g 2 is less than about 3 percent. Therefore, the new method works quite accurately for identifying the liquid properties.

CONCLUSIONS

In this paper, a new mechanical model for micro-cantilever rheological sensor in liquid environment has been derived and a novel measurement technique based on non-linear estimation technique has been proposed which is significantly different from the classical detection techniques, since it is based on time domain measurement. Simulation results show that this novel tech- nique provide accurate estimation performance of the characteristic parameters of the liquid. Therefore, this novel technique can be considered for the next generation biochemical microcantilever sensor applications.

FIGURE 1 .

 1 FIGURE 1. MICRO-CANTILEVER SENSING SYSTEM.

FIGURE 2 .

 2 FIGURE 2. FORCE DIAGRAM OF AN INFINITESIMAL ELE-MENT OF THE BEAM.

  ΛL)e -ΛL + sin(ΛL)e -ΛLcos(ΛL)e ΛL + sin(ΛL)e ΛL EIΛ 3 [2 + cos(ΛL)(e ΛL + e -ΛL )] F

FIGURE 3 .

 3 FIGURE 3. DEFLECTION MEASUREMENT AND ESTIMATION AT THE FREE-END.

FIGURE 4 .

 4 FIGURE 4. PERCENTAGE ERROR OF THE CHARACTERISTIC PARAMETERS.

TABLE 1 .

 1 PHYSICAL PARAMETERS

	Items	Values
	Excitation frequency	50 Hz
	Sampling frequency	1000 Hz
	Size of micro-cantilever	100µm × 20µm × 2µm
	Density of the beam	2330kg/m 3
	Young's elastic modulus	169 GPa
	Quality factor	1.681
	Resonant frequency	2.425kHz
	Amplitude of the excitation force 1µN

TABLE 2 .

 2 SIMULATION VALUES

	Items	g 1	g 2	ρ	η
	Real Values	1	1 × 10 -4	1.931 × 10 3 0.225
	Estimated Values 1.029 1.071 × 10 -4 2.153 × 10 3 0.232