Benoît Chèze
email: benoit.cheze@ifpen.fr.

Charles Collet

Anthony Paris

Victor Court
email: victor.court@ifpen.fr

Estimating discrete choice experiments: theoretical fundamentals

Keywords: Revealed preference theory, Stated Preference / Stated Choice methods, Discrete Choice Modelling, Discrete Choice Experiment

La collection "Les Cahiers de l'Économie" a pour objectif de présenter les travaux réalisés à IFP Energies nouvelles et IFP School qui traitent d'économie, de finance ou de gestion de la transition énergétique. La forme et le fond peuvent encore être provisoires, notamment pour susciter des échanges de points de vue sur les sujets abordés. Les opinions exprimées dans cette collection appartiennent à leurs auteurs et ne reflètent pas nécessairement le point de vue d'IFP Energies nouvelles ou d'IFP School. Ni ces institutions ni les auteurs n'acceptent une quelconque responsabilité pour les pertes ou dommages éventuellement subis suite à l'utilisation ou à la confiance accordée au contenu de ces publications.

Pour toute information sur le contenu, contacter directement l'auteur.

Introduction

Discrete Choice Experiments (DCEs) have gained popularity among environmental economists in the last ten years, as shown by its considerable recent use in many fields. DCE is a Stated Preference (SP) / Stated Choice (SC) method in which preferences are elicited through repeated fictional choices made by respondents [START_REF] Hoyos | The state of the art of environmental valuation with discrete choice experiments[END_REF][START_REF] Louviere | Stated Choice Methods, Analysis and Applications[END_REF]. In comparison to Revealed Preference (RP) approaches, it presents the advantage of (i) allowing an ex ante evaluation of public policy scenarios and (ii) capturing the non-use value. 1Compared to other SP methods, such as the contingent valuation approaches, DCEs take into account several characteristics, called attributes in the literature, of the considered issue. It is particularly useful in examining the trade-offs that occur in decision-making. As DCE framework aims at considering several attributes of an issue, it delivers more detailed information than other SP methods. Especially, it makes it possible to estimate the marginal rates of substitution (MRS) between different attributes. When one of these attributes corresponds to a monetary contribution, these MRS can be interpreted as the willingness To pay (WTP), or To Accept (WTA), for changes in the attributes' levels. It thus procures a WTP/WTA for each of these attributes, rather than a global WTP/WTA, as a contingent valuation would do.

This multi-dimensional framework is especially useful for the valuation of hypothetical public policies as it not only allows to values respondent' preferences globally but also every attribute of the problem. It thus can help in setting priorities for public action by identifying the most important factors for a given population. It exhibits the trade-offs at stake for people involved. It can also identify the specific reason a public program is either successful or unsuccessful. DCE approaches rely on the economic theories of consumer choice and non-market valuation. In a DCE survey, respondents must choose from several options defined by their attributes (i.e., fundamental characteristics of the respondents' situation). Respondents receive a series of choice sets, each comprising several options. Often, three options are presented: nothing changes (i.e., the status quo) and two alternative options. The use of an opt-out option (status quo) is known to improve realism in choices [START_REF] Adamowicz | Future directions of stated choice methods for environment valuation[END_REF][START_REF] Kontoleon | Assessing the impacts of alternative 'opt-out'formats in choice experiment studies: consumer preferences for genetically modified content and production information in food[END_REF]. Respondents then choose their favorite option. Each option has different levels of the attributes. One of these attributes usually represents the monetary contribution of the respondents. Other attributes can include environmental or social implications of the issue under consideration. See [START_REF] Louviere | Stated Choice Methods, Analysis and Applications[END_REF] for a detailed description of the method.

Typically, a DCE survey begins with some information about the issue under investigation (e.g. advantages and disadvantages). Various successive choices (say ten for instance) are then proposed between two scenarios and -in most cases -a status quo option presented in the form of a set of (ten) choice cards. Fictional choices are thus made by respondents among several options. Respondents finish survey by responding to some socio-demographic, economic and/or psychological (attitudinal) questions allowing researchers to analyze the impact of these characteristics on their structure of preferences. This quantitative analysis is conducted using specific econometric estimators relying on Discrete Choice Modelling (DCM) theory.

These estimators are presented in this working paper. Section 2 presents theoretical foundations of the DCE approach, both in preference and WTP spaces, and specify the indirect utility functions usually estimated in the litterature whether or not alternatives are labeled. Section 3 presents estimators derived from econometric models used to analyze stated choices proposed in DCE surveys to reveal respondents's preferences. Section 4 concludes.

Theoretical foundations of the Discrete Choice Experiment approach

The DCE modeling framework relies on Lancaster's characteristics of value theory [START_REF] Lancaster | A new approach to consumer theory[END_REF] as well as random utility theory [START_REF] Mcfadden | Conditional logit analysis of qualitative choice behaviour[END_REF]. This Section presents both the use of the Random Utility Model (RUM) in DCEs and how indirect utility functions are usually specified.

Random Utility Model in Discrete Choice Experiment

Since a good may be described by a set of characteristics [START_REF] Lancaster | A new approach to consumer theory[END_REF], its value becomes the sum of the values of all its characteristics. According to the DCE approach, an alternative i ∈ 1; I can then be described by a set of K observable characteristics, called attributes, as denoted by X i = (x i,1 , ..., x i,k , ..., x i,K ) . An individual n ∈ 1; N is described by A economic and attitudinal characteristics, called socioeconomic variables, denoted Z n = (z n,1 , ..., z n,a , ..., z n,A ) . The (indirect) utility V n,i is thus given by:

V n,i = V (X i , Z n ) ∀ n ∈ 1; N ; i ∈ 1; I (1)
McFadden (1974) considered that individuals make choices according to a deterministic component based on both their measured characteristics Z and attribute alternatives X, along with some degree of randomness . The random utility U n,i of an alternative i ∈ 1; I for respondent n ∈ 1; N is therefore composed of a deterministic part, i.e. the (indirect) utility V n,i = V (X i , Z n ), and the stochastic element, n,i , thereby capturing the unsystematic and unobserved random element of respondent n's choice [START_REF] Louviere | Stated Choice Methods, Analysis and Applications[END_REF].

U n,i = V (X i , Z n ) + n,i (2)
Assuming the rationality of individuals, respondents choose the alternative i from a finite set of I alternatives, if its utility, U n,i , is greater than the utility derived from any other alternatives j, U n,j :

U n,i > U n,j ⇒ V n,i + n,i > V n,j + n,j ∀ j = i ; (i, j) ∈ 1; I 2 (3)
The probability that a respondent chooses alternative i is the same as the probability that the utility of alternative i is greater than the utility of any other alternative in the choice set [START_REF] Adamowicz | Stated preference approaches for measuring passive use values: Choice experiments and contingent valuation[END_REF]. Following [START_REF] Train | Discrete Choice Methods with Simulation[END_REF], this probability is

P n,i = P U n,i > U n,j ∀ j = i ; (i, j) ∈ 1; I 2 (4) ⇔ P n,i = P V n,i + n,i > V n,j + n,j ∀ j = i ; (i, j) ∈ 1; I 2 (5) ⇔ P n,i = P n,j -n,i < V n,i -V n,j ∀ j = i ; (i, j) ∈ 1; I 2 (6)

Specification of the indirect utility functions

According to Equation (2), the random utility U n,i is composed of a deterministic component, V n,i = V (X i , Z n ), and a stochastic element, n,i . Before estimating an econometric model, the deterministic part of the utility function V n,i = V (X i , Z n ) must be specified. The linear specification is often chosen in the literature for the sake of simplicity. We have thus introduced the column vector of parameters β n = (β 1,n,i , ...β K,n,i ), which are the coefficients quantifying the (linear) influence of the K attributes on utility that individual n associates with the I available alternatives.

An Alternative Specific Constant (ASC i ) term is usually introduced to capture the effect of unobserved influences (i.e. omitted variables) on the utility function for the i-th alternative. It is a dummy variable assuming a value 1 for the i-th alternative, and 0 otherwise. In its most general form, the model is specified such that the utility of an alternative is expressed as follows:

U n,i = (η i + φ n,i Z n )ASC i + β n,i X i + (X i ) Γ i Z n + n,i (7)
where φ n,i = (φ n,i,1 , ..., φ n,i,A ) are the coefficients representing the direct influence of the A socioeconomic and psychological variables on the utility of the N individuals for alternative I. The matrix Γ i of size (K × A), is composed of coefficients γ i,k,a ; it serves to capture the effect of individuals' characteristics z n,a on attribute x i,k .

Non labeled choice experiment

When alternatives are unlabeled, the choice options given are thus symmetrical, differing only in the attributes assigned X i . Even if a reference choice bias were to exists, it could be represented by a simple Alternative Specific Constant. Therefore, the previous coefficients β n,i , φ n,i and Γ i do not depend on alternative i; they can be simplified into β n , φ n and Γ. Let's observe that for each individual, φ n Z n remains constant across all alternatives i and therefor does not influence the model; hence, this item can be removed, yielding:

U n,i = β n X i + (X i ) ΓZ n + n,i (8)
This underlines that for a non-labeled choice model, one of the main solution in order to include the influence of a socio-economic variable z n,a is to introduce it through an interaction term with the alternative attributes x k,i .

WTP-space approach

In their seminal paper, [START_REF] Train | Discrete choice models in preference space and willingness-to-pay space[END_REF] found that the WTP distributions they derived from preference models had an "unreasonably" large variance in comparison with WTP distributions derived from utility models in the WTP-space. This result has been confirmed in subsequent studies [START_REF] Mabit | Representation of taste heterogeneity in willingness to pay indicators using parameterization in willingness to pay space[END_REF][START_REF] Scarpa | Utility in willingness to pay space: a tool to address confounding random scale effects in destination choice to the alps[END_REF][START_REF] Rose | A comparison of the impacts of aspects of prospect theory on wtp/wtaestimated in preference and wtp/wta space[END_REF]. According to the WTPspace approach, marginal WTP estimates are directly obtained, unlike with the preference space approach, where a ratio is to be computed between the non-cost attribute and the cost attribute.

Welfare measures can be determined in the form of marginal WTP/WTA, i.e. for a unit change of a given attribute, by computing the MRS between the considered attribute and the cost attribute [START_REF] Louviere | Stated Choice Methods, Analysis and Applications[END_REF]. Since utilities are modeled as linear functions of the attributes, the MRS between two attributes is the ratio between the corresponding coefficients.2 Thus for continuous attributes3 , the WTP for an attribute k ∈ 1; K -1 equals the ratio of the attribute's coefficient to the cost coefficient, β cost,n :4

W T P k,n = dx cost,i dx k,i = dU n /dx k,i dU n /dx cost,i = ∂V n /∂x k,i ∂V n /∂x cost,i = β k,n β cost,n (9)
where x cost,i is the cost attribute, i.e. the K-th attributes, x K,i , in X i .

Equation ( 8) corresponds to the preference-space model. Combined with Equation ( 9), utility (without considering any socio-demographic characteristic effects) can be rewritten in the WTPspace as follows:

U n,i = β cost,n x cost,i + K-1 k=1 W T P k,n x k,i + n,i (10)
When adding the effects of individuals' characteristics, it is necessary for all γ k,a (which measure the effects of individuals' characteristics on attribute preferences) to be divided by the cost attribute's coefficient β cost,n in order to generate their estimation in the WTP-space. 5 The newly formed γ k,a in the WTP-space therefore measure the effect of the individuals' characteristic z n,a on the WTP for attribute k. It is possible to consider these effects as being directly included in individuals' WTP so as to facilitate model interpretation.

The WTP is now defined as follows:

W T P k,n = W T P base,k,n + A a=1 γ k,a z n,a with W T P base,k,n = β k,n β cost,n (11)
Equation ( 10) becomes:

U n,i = β cost,n x cost,i + K-1 k=1 W T P base,k,n + A a=1 γ k,a z n,a x k,i + n,i (12)
3 Econometric models

Different econometric models can be used to analyze discrete choice data. Indeed, the heterogeneity of preferences can be modeled to different degrees, depending on (i) assumptions made on the distribution of the column vector of parameters β n = (β 1,n,i , ...β K,n,i ) and (ii) the way socioeconomic variables are introduced. This Section presents Conditional Logit (CL), Random Parameter Logit (RPL), Latent Class (LC) and Integrated Choice and Latent Variable (ICLV) models.

Conditional Logit models

The CL model is a cornerstone for analyzing discrete choice data and has been widely used in DCE. This model has however several well-known limits. An important drawback is that it assumes homogeneous preferences across respondents, meaning that the probability that an agent n chooses alternative i in a choice set S, is considered fixed across all individuals (β n = β, Γ n = Γ ∀ n), while we can expect the preferences to vary. Two other important drawbacks are uncorrelated unobserved components and the hypothesis of the independence of irrelevant alternatives (IIA).

IIA implies that the relative probabilities of two options being chosen are unaffected by the introduction or removal of other alternatives. The (multinomial) logit probability that a respondent n chooses a particular alternatives i is6 :

P n,i = e Vn,i j e Vn,j (13)
The hypothesis of irrelevant alternatives implies that the relative probabilities of two options, here i and h, being chosen are unaffected by the introduction or removal of other alternatives. This is an underlying assumption of the CL as shown by Equation ( 14):

P n,i P n,h = e V n,i j e V n,j e V n,h j e V n,j = e Vn,i e V n,h (14)
If the IIA property is violated then the CL model does not fit the data. Results will be biased, leading to unrealistic predictions, and hence a discrete choice model that does not require the IIA property should be used. Moreover, it does not uses the information that the same individual make several choices, thus cannot be as accurate as a model taking into account this supplementary information.

Random Parameter Logit models

Compared to the CL model, the RPL model [START_REF] Mcfadden | Mixed MNL models for discrete response[END_REF][START_REF] Train | Discrete Choice Methods with Simulation[END_REF] relaxes the IIA hypothesis and is able to take the heterogeneity of preferences into account. Indeed, the preferences parameters β are allowed to vary randomly across respondents, in allowing for the fact that different decision-makers may have different preferences:

β n = β m ∀ n = m ; (n, m) ∈ 1; N 2 .
As such, conditional on the individual-specific parameters and error components, we can define the (multinomial) logit probability7 that respondent n chooses a specific alternative i for a given β:

P n,i |β = L n,i (β) = e Vn,i(β) j e Vn,j (β) (15)
Without taking into account the cross-effects of socio-demographic characteristics, the individual specific utility is simply modeled by:

U n,i = β n X i + n,i (16)
where n,i ∼ iid extreme value type I, and β n ∼ g(β|Ω).

Since β n is not known, the unconditional choice probability of person n choosing alternative i is the integral of P n,i |β over the distribution of β:

P n,i = L n,i (β)f (β|Ω)dβ (17)
f (β|Ω) is the density describing the distribution of preferences across individuals. Ω is a vector of the true parameters of the taste variation, e.g. fixed parameters of the distribution representing the mean and standard deviation of β n within the population. The true distribution remains unknown, so, in theory, any distribution could be applied [START_REF] Hensher | The Mixed Logit model: The state of practice[END_REF]. A common assumption made in the literature is to assume that random parameters are normally distributed, but the one associated with the cost attribute.8 The latter is usually kept fixed or supposed to follow a log-normal distribution in order to avoid a "wrong" sign (i.e. negative) for a share of respondents.

The choice probability in equation ( 17) cannot be calculated exactly because the integral does not have a closed form solution in general. This integral is approximated through simulations. For a given value of the parameters Ω, a value of β is drawn from its distribution. Using this draw, the logit formula in ( 15) is calculated. This process is repeated for many draws, and the mean of the resulting L n,i (β) is taken as the approximate choice probability yielding equation (18):

SP n,i = 1 R R r=1 L n,i (β r ) ( 18
)
where R is the number of draws of β, and SP is the simulated probability that an individual n chooses alternative i.

With a RPL model, one can consider that parameters are random but not independent between them. [START_REF] Train | Discrete Choice Methods with Simulation[END_REF] describes how the simulation and the estimation of the parameters can be done when preference parameters are assumed to follow a random normal joint distribution β ∼ N (µ, Σ) where µ is the mean vector of parameters and Σ the covariance matrix for each individual.

Latent Class models

Another way to relax the hypothesis of the IIA and to take account for the heterogeneity in respondents' preferences is to estimate a LC model. In the latter, each respondent is sorted into a number of classes C in which preferences are assumed to be homogeneous with respect to attributes. Preferences are allowed to be heterogeneous between each latent class segment c ∈ C. 15), the logit probability that respondent n prefers a specific alternative i over alternatives j is no more defined for a given β but becomes conditional on class c. Indeed, the β's are now assumed to follow a discrete distribution and belong to one class c of C classes. Thus, the conditional probability that respondents who are members of class c choose alternative i is:

Compared to Equation (

P n,i |β c = e Vn,i(βc) j e Vn,j (βc) ∀ c ∈ 1, ..., C (19)
where β c is the vector of preferences parameters specific to each class c, representing the average importance of each attribute for respondents belonging to c.

The unconditional probability of individual n selecting choice option i can be expressed as:

P n,i = C c=1 (Π n,c P n,i |β c ) = C c=1 Π n,c e βcXi j e βcXj ( 20
)
where Π n,c is the probability of membership of respondent n in class c:

Π n,c = e φcZn h e φ h Zn (21)
where Z n is the vector of psychometric constructs and socioeconomic characteristics, and φ is the vector of parameters associated to Z n [START_REF] Boxall | Understanding Heterogeneous Preferences in Random Utility Models: A Latent Class Approach[END_REF]. According to Equation ( 21), the probability of belonging to a class c with specific preferences is probabilistic, and depends on the social, economic and attitudinal characteristics of the respondents. Combining Equation ( 20) and Equation ( 21), it comes that the LC model assumes that respondent characteristics affect his choice indirectly through their impact on segment membership. Note that φ c includes C -1 class membership parameters with φ C being normalized to zero for identification. All other coefficients φ c are thus interpreted relative to this normalized class.

Integrated Choice and Latent Variable models

General Mathematical specification

Hybrid models are composed with a Discrete Choice Model (DCM) and a Structural Equation Model (SEM). The DCM is very similar to the previous models, the difference being in the introduction of the Latent Variables into the deterministic utility function V :

U n,i = V (X i , Z n , LV n ) + n,i (22)
where

LV n = (LV n,1 , ..., LV n,Q ) is a (Q × 1)
vector containing the Q latent variables for individual n and n = ( n,1 , ..., n,I ) is the (I × 1) vector of iid error terms, for the I alternatives, following a extreme value distribution type I with 0 mean and a covariance matrix Σ ,n .

The structural equation of the SEM gives information on the distribution of the latent variables given the observed socio-economic characteristics:

LV n,q = S(Z n ) + v n,q (23)
where

v n = (v n,1 , ..., v n,Q ) is the (Q × 1) vector of iid error terms for the Q latent variables, i.e. v n ∼ N (0, Σ v,n ), with Σ v,n the covariance matrix.
The measurement equation of the discrete choice part of the model is the same as the other DCM:

Y n,i = 1 if U n,i > U n,j ∀ j = i ; (i, j) ∈ 1; I 2 0 otherwise (24)
The measurement equation of the SEM uses the values of the attitudinal indicators as dependent variables, and explain their values with the help of the latent variables:

I * n,l = M (LV n ) + w n,l (25)
where I * n is a (L × 1) vector containing the L attitudinal indicators used to measure the latent variables. w n = (w n,1 , ..., w n,L ) is the (L × 1) vector of iid error terms for the L attitudinal indicators, i.e. w n ∼ N (0, Σ w,n ), with Σ w,n the covariance matrix.

The attitudinal indicators are almost always measured with ordinal variables based on Likertscale items. The most simple option is to consider these ordinal variables as continuously normally distributed [START_REF] Bouscasse | Integrated choice and latent variable models: A literature review on mode choice[END_REF]. However, a more accurate way to take into account the ordinal nature of the indicators is to add a threshold model equation to the measurement model of the SEM, which is usually done in the SEM literature [START_REF] Kamargianni | Investigating the subjective and objective factors influencing teenagers' school travel mode choice -an integrated choice and latent variable model[END_REF][START_REF] Giansoldati | Does electric car knowledge influence car choice? Evidence from a hybrid choice model[END_REF]. Therefore I * n,l is a latent variable coming from a threshold model equation:

I n,l =        1 if -∞ < I * n,l ≤ τ 1,l 2 if τ 1,l < I * n,l ≤ τ 2,l ... D if τ D-1,l < I * n,l ≤ ∞ (26)
I * n,l is the exact but unobserved dependant variable and D is the total number of categories (D could vary between all attitudinal indicators L, but here was supposed constant to simplify notations). τ d,l are the threshold for I * n,l that give the probability of observing the answer I n,l . The probability that the answer I n,l = d is equal to

P (I n,l = d) = P (τ d-1,l < I * n,l ≤ τ d,l ) = F v (τ d,l ) -F v (τ d-1,l ) (27)

Specification for a non labeled choice experiment

To simplify notations, we consider here that V depends only on the attribute alternatives X i , so that the individual characteristics Z n are only included in the Structural equation of the SEM. A linear specification is adopted and the RUM equation becomes:

U n,i = β n,i X i + ξ n,i LV n + X i Λ n,i LV n + n,i (28)
The matrix Λ n,i of size (K × Q) is composed of coefficients λ n,i,k,q , capturing the cross-effect of latent variable LV n,q on attribute x i,k .

In the case of non labeled DCE, the previous coefficients do not depend on the alternative i. Therefore they can be simplified as β n , ξ n and Λ n . ξ n LV n can be removed because constant between the alternatives i. This leads to the following equation:

U n,i = β n X i + (X i ) Λ n LV n + n,i (29)
Like in ICLV empirical models found in the literature, structural (23) and measurement (25) equations are usually chosen as linear specifications to limit the number of parameters and enable convergence of the algorithm.

LV n,q = γ n,q Z n + v n,q (30)

I * n,l = ν n,l LV n + w n,l (31)
To enable better understanding of psychological mechanisms, attitudinal indicators are assumed to be explained by only one latent variable: the vector ν n,l is composed of Q -1 zeros and only one non-zero coefficient for the corresponding latent variable. This assumption is made in every empirical ICLV studies we consulted.

Maximum likelihood techniques are used to estimate the unknown parameters. Assuming the error components ( n , v n , w n ) are independent, the joint probability of the observable variables I n and Y n , conditional on the exogenous variable X and Z n is:

f Y,I (Y n , I n |X, Z n ; β n , Λ n , γ n , ν n , Σ ,n , Σ v,n , Σ w,n ) = LV P (Y n |X, Z n , LV n ; β n , Λ n , Σ ,n ) f LV (LV n |X, Z n ; γ n , Σ v,n )f I (I n |LV n ; ν n , Σ w,n )dLV n (32)
with f LVn and f I * the density functions depending on the error terms γ n and ν n . The first term of the integral corresponds to the choice model, the second term corresponds to the measurement equation, and the third term to the structural equation.

Conclusion

Models presented in this working paper belong to the Generalized Extreme Value (GEV) models. This large class of models is adapted to the case where the variable to be explained is a qualitative variable whose cannot be ranked in relation to each other. They aim at explaining people's choice where the choice set is exhaustive, made up of finite number of mutually exclusive alternatives described by attributes. Although seemingly restrictive, these conditions are found in many cases in everyday life such as mode of transport and/or car ownership, home location, food choices, medical or voting decisions, etc. Despite an RP resurgence, most discrete choice applications are currently based on SP/SC data where respondents have to choose between different hypothetical scenarios. In a DCE survey respondents are faced, by construction, with multiple choice situations which allow researcher to obtain more information per respondent. Analyst can control settings, in the manner of experimental sciences, by (i) choosing not only attributes but also their levels and (ii) designing SP surveys that encourages trade-off behaviour. The latter is an area of research in its own right. It aims to obtain as much information as possible with as few choice situations as possible. Indeed, respondents cannot be offered all the possible combinations between the different levels of the different attributes as they would get tired and would eventually not respond with sufficient concentration.

Compared to RP data, where real world choices are analysed, SP methods produce data that can look at hypothetical choice scenarios by proposing to choose alternatives that do not exist yet.

It is a strength as it is one of the few methods that allow for an ex ante evaluation of public policies that do not yet exist. It may also be a curse as respondents make fictional choices rather than real ones. This drawback is referred to in the literature as hypothetical bias. This is a critical issue: if surveys are not taken seriously by respondents, it eventually yields to non reliable data. Hypothetical bias arises when choice tasks are difficult to complete or choices are not incentivecompatible. One solution to the latter is to emphasize consequentiality of the choices (referred to as "introducing consequentiality") -and to control for -in the survey, which is a necessary condition for incentive compatibility (see [START_REF] Johnston | Contemporary guidance for stated preference studies[END_REF].

P n,i = +∞ s=-∞
exp -e -s j e -(Vn,i-Vn,j ) e -s ds (42)

We now define t = e -s . The expression -e -s ds therefore gives dt and note that t approaches zero (resp. positive infinity) if s tends to infinity (resp. negative infinity) as:

P n,i = 0 t=+∞ exp -t j e -(Vn,i-Vn,j ) (-dt) (43) that is to say: Vn,i-Vn,j ) dt (44

P n,i = +∞ t=0 exp -t j e -(
)
This expression is now easy to integrate and allows us to obtain expression in equation (15).

P n,i = exp( -t j e -(Vn,i-Vn,j ) )

j e -(Vn,i-Vn,j ) +∞ 0 (45) Vn,i-Vn,j ) (46)

P n,i = 0 - 1 -j e -(
P n,i = 1 j e -(Vn,i-Vn,j ) = e Vn,i j e Vn,j

where P n,i , the (multionomial) logit probability, only depends on observable components.

The CL model is estimated using maximum likelihood procedures. The probability that a respondent n chooses a particular alternative is i (P n,i ) yn,i with y n,i = 1 if the alternative i is chosen and zero otherwise. Assuming the independence in choices of each respondent, the likelihood and log-likelihood functions are given by:

L(β) = N n=1 i (P n,i ) yn,i (48)

The non-use value of non-market goods is the existence value or bequest value. It is useful when valuing, for instance, issues linked to biodiversity and climate change.

It should be noted that the derivative of the unobserved part of the utility function is supposed to be zero with respect to both attributes.

See Appendix B for attributes modeled as dummy variables and specific cases of both RPL and LC models.

The cost attribute is considered to be the last one, for purpose of consistency with the previous notations.

We have opted to maintain the same variable for the effect of the individuals' characteristics in both the preference and WTP-space in order to simplify notations.

See Appendix A for more details.

As the error term is assumed to be independent and identically distributed (iid) Type I Extreme Value. Note that Appendix A details calculation to obtain this probability.

Random parameters are generally supposed to be normally distributed in the RPL model because it is the most easily applied distribution allowing for both negative and positive preferences.

Appendices

Appendix A The multinomial logit model Assuming n,i being iid and following a type I extreme-value distribution, i.e. a standard Gumbel distribution, we thus specify a conditional logit model (or multinomial logit model).

The cumulative distribution function F and the density function f of each n,i are given by:

Since the unobserved components are independent, we can multiply Eq. ( 6) to obtain the probability of individual n choosing alternative i, conditional on n,i :

The non conditional probability for an agent n to choose the alternative i is therefore the integration of P n,i | n,i over the distribution of n,i :

By replacing n,j with s, equation (37) becomes:

e -e -(V n,i -V n,j +s) e -s e -e -s ds (38)

As V n,i -V n,i = 0, we have:

and the last term can be introduced into the product,

By removing the first exponential from the product, we obtain:

Appendix B Computing WTP/WTA: some complements

For attributes modeled as dummy variables, the W T P l k,n associated with attribute k and category l is

The variation of utility associated with a change of the attribute k from the status quo to category l is a measure of the willingness to pay to do so.

Since the RPL model assumes attributes' coefficients are randomly distributed, we then have the following result in this case:

Hole ( 2007) explains "since the logit model is typically estimated using maximum likelihood, which implies that the coefficients in the model are asymptotically normally distributed, it is reasonable to assume that WTP is given by the ratio of two normally distributed variables when the model is estimated using a large sample. The distribution of the ratio of two normally distributed variables has been derived by [START_REF] Fieller | The distribution of the index in a normal bivariate population[END_REF] and [START_REF] Hinkley | On the ratio of two correlated normal variables[END_REF], who show that the distribution is approximately normal when the coefficient of variation of the denominator variate (in this case β cost,n ), is negligible."

For a LC model, the WTP for the individual n for a variation of the attribute k can be computed per class as: For each models, the estimated standard deviations and confidence intervals around the mean of the WTP estimates may be obtained using different methods, depending whether or not the previous property is verified. The four main methods are i) the delta, ii) the Fieller [START_REF] Fieller | Some problems in interval estimation[END_REF], iii) the Krinsky and Robb (Krinsky andRobb, 1986, 1990) and iv) the bootstrap methods. In calculating a WTP, it is important that both parameters used in the calculation be statistically significant, otherwise no meaningful WTP measures can be established.


Llyes Hadbaoui

Abdelhakim Senoussi

Johann Huguenin

IAV A Araba

Hassan Ii

Morocco A López-Francos

Iamz-Ciheam Spain

Feeding practices in sheep stockbreeding in Algeria steppe (M'Sila region): Current situation and prospects for improvement

Keywords: Llyes Hadbaoui, Abdelhakim Senoussi, Johann Huguenin. Feeding practices in sheep stockbreeding

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

S3-12

Feeding practices in sheep stockbreeding in Algeria steppe (M'Sila region): Current situation and prospects for improvement , 2016). For the past 50 years, it knowing many changes and transformations resulting from the stockbreeding management that contributed to a decline in forage productivity and the surface of steppe pastures. In this way, steppe pasture covers only 30 to 40% of the animal's feed needs.

In this situation, stockbreeders are adopting various feeding strategies. Our objective was to assess the feeding practices adopted by the herders in M'Sila region. This region has 1 million ha of rangelands and hosts 1.63 million sheep heads (DSA, 2016). One hundred herders, from the different breeding systems and different pastoral areas, were investigated by technical-organizational approaches. The analysis of field interviews revealed the presence of two feeding calendars, involving different feed sources (concentrate feeds, cultivated forages, and cereal co-products). Although the intake of concentrated feeds is a recurrent practice to meet animal needs (about 40% of ration), there are new feeding trends based on the direct consumption of green fodder (barley, oat and alfalfa). The choice of feeding strategies to be undertaken depends on the socio-economic and agricultural considerations of each herder (herd size, financial capacity, farming practice, income diversification, etc.).

Keywords. Agropastoral stockbreeding systems -Steppe rangelands -M'Sila Algeria -Feeding calendar -Feed supplement -Prospects.

Pratiques d'alimentation dans les élevages ovins en steppe algérienne (Région de M'Sila) : Réalité et perspectives d'amélioration

Résumé. Le territoire steppique algérien de 20 millions d'hectares, se délimite entre les isohyètes de 100 et 400 mm/an. Il abrite 17 millions d'ovins (MADRP, 2016)

Hadbaoui Ilyes 1,2 , Abdelhakim Senoussi 1 and Johann Huguenin 3 1

  Saharan Bioresources Laboratory: Conservation and Recovery. University kasdi Merbah Ouargla, Algeria 2 Centre for Scientific and Technical Research on Arid Regions (CRSTRA) Taouiala Station, Biskra, Algeria 3 CIRAD, UMR Selmet, CIRAD, INRA, Montpellier SupAgro (Univ Montpellier), 34398 Montpellier, France Abstract. The Algerian steppe territory covers 20 million hectares and is located between the isohyets of 100 and 400 mm/year. It is home to17 million sheep (MADRP

  . Depuis 50 ans il a connu de nombreuses évolutions et transformations engendrant des conduites d'élevage qui contribuent à une baisse de la productivité fourragère et de la surface des pâturages steppiques. Ainsi, la pâture steppique couvre seulement 30 à 40% des besoins alimentaires des animaux. Pour faire face à cette situation, les éleveurs adoptent diverses stratégies alimentaires.

	Notre questionnement s'intéresse aux pratiques et modalités d'alimentation adoptées par les éleveurs de la ré-
	gion de M'Sila. Cette région compte 1 million d'ha de parcours et 1,63 million d'ovins (DSA, 2016). Cent éleveurs
	des différents systèmes d'élevage et de différentes zones pastorales ont été enquêtés par approches technico-
	organisationnelles. L'analyse des données des entretiens de terrain révèlent la présence de deux calendriers ali-
	mentaires, avec différentes tendances fourragères (aliments concentrés, fourrages cultivés, et sous-produits de
	la céréaliculture). Bien que l'apport des aliments concentrés relève d'une pratique récurrente pour satisfaire les
	besoins des animaux (soit environ 40% de la ration), il existe des nouvelles tendances alimentaires basées sur
	les fourrages cultivés et consommés généralement en vert (orge, avoine et luzerne). Le choix des stratégies ali-

mentaires à entreprendre dépend de considérations socio-économiques et agricoles de chaque éleveur (taille de troupeau, capacité financière, pratique de l'agriculture, diversification des revenus, etc.).

Mots-clés. Systèmes d'élevage agropastoraux -Parcours steppiques -M'Sila

Algérie -Calendrier alimentaire -Complément alimentaire -Perspectives.

S. Ates, Oregon State Univ., USA G. Luciano, Univ. Catania, Italy A. Bernués, CITA-Aragón, Spain D. Martín Collado, CITA-Aragón, Spain M. Bengoumi, FAO-SNE, Tunisia Y. Mena, Univ. Sevilla, Spain H. Ben Salem, INRAT, Tunisia E. Molina-Alcaide, EEZ-CSIC, Spain I. Casasús, CITA-Aragón, Spain G. Moreno, Univ. Extremadura, Spain M. Chentouf, INRA, Morocco M. Mounsif, ENA-Meknès, Morocco P. de Frutos, CSIC, Spain M. Napoléone, INRA, France A. de Vega, Univ. Zaragoza, Spain P. Nichols, Univ. Western Australia, Australia P. Gaspar, Univ. Extremadura, Spain C. Ovalle, INIA, Chile S. Giger-Reverdin, AgroParisTech, France F. Pacheco, DRAP-Norte, Portugal I. Hadjigeorgiou, Agr. Univ. Athens, Greece C. Porqueddu, CNR-ISPAAM, Italy M. Jouven, SupAgro, France A. Priolo, Univerity of Catania, Italy A. Keli, ENA-Meknès, Morocco M. Rekik, ICARDA, Tunisia A. Kyriazopoulos, Democritus U. Thrace, Greece P. Swanepoel, Stellenbosch Univ., South Africa L. López Marco, IAMZ-CIHEAM, Spain D. Yáñez, CSIC, Spain A. Zoghlami, INRAT, Tunisia

Joint Meeting of the FAO-CIHEAM Network for Research and Development in Sheep and Goats (Subnetworks on Nutrition and Production Systems)


Christophe Moulin
email: christophe.moulin@univ-st-etienne.fr

Cecile Barat
email: cecile.barat@univ-st-etienne.fr

Mathias Géry
email: mathias.gery@univ-st-etienne.fr

Christophe Ducottet
email: ducottet@univ-st-etienne.fr

Christine Largeron
email: largeron@univ-st-etienne.fr

Cecile Barat

Ujm At

Cécile Barat

UJM at ImageCLEFwiki 2008

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The capacity of data storage increases constantly making possible the collection of large amount of information of all kinds, as texts, images, videos or combinations of them. In order to retrieve documents in such amount of data, information retrieval techniques tailored for the data types are required.

The ImageCLEF collection consists of multimedia documents made up of text and images. In this paper, we present our participation to the ImageCLEFwiki task. Our research goals are twofold: First, we aim to propose a multimedia document model combining text and image modalities adapted for multimedia retrieval. Second, we want to study the performance of our model compared to a text retrieval approach. In order to benefit from our long time experience with textual model, we develop a vector-based model composed of textual and visual terms. The textual terms correspond to words of the text. The visual terms are obtained through a bag of words approach. Local colour descriptors are extracted from images and quantized by k-means leading to an image vocabulary.

After presenting our model, we will describe the runs we submitted. Then, we will comment on the results we obtained and conclude.

Visual and textual document model

ImageCLEFwiki is a multimedia collection where documents are composed of text and image. User needs are represented by queries ("topics"), which are also multimedia (text, image and concept). Hence a multimedia document model is necessary to handle such a collection. We focus our work on combining textual and visual information without using the concept field of the topics. Before explaining our model, we will describe the collection data.

Description of the data: ImageCLEF Wikipedia collection

The ImageCLEF Wikipedia collection is composed of 151'519 multimedia XML documents and 75 multimedia topics. The documents are made up of an image and a short text. Images are in a common format (jpeg and png) and their sizes are heterogeneous. They depict either drawings, paintings or screenshots. They are in colour or in black and white. The textual part of a document is unstructured and consists of a description of the image, information about the Wikipedia user who has uploaded the image, or the copyright of the image. The average number of words per document is about 33 words.

It is also possible to get the whole Wikipedia article from which each image is extracted. The ImageCLEF Wikipedia collection provides 75 topics whose only 28 contains an image as example of the expected answer. For the text part, these queries were composed of about 2 or 3 words.

Textual representation model

One of the most known document model in textual information retrieval is the vector space model introduced by Salton and al. [START_REF] Salton | A vector space model for automatic indexing[END_REF]. This model is based on a textual vocabulary T = {t 1 , ..., t j , ...t |T | }. Each document is represented as a vector of weights d i = (w i,1 , ..., w i,j , ..., w i,|T | ) where w i,j is the weight of the term t j in the document d i . In order to calculate the weight of a term t j in a document d i , a tf.idf formula is usually applied. The tf i,j (term frequency) measures the relative frequency of a term t j in a document d i . We have used the one defined in the Okapi formula from Robertson and Jones [START_REF] Robertson | Okapi at trec-3[END_REF]:

tf i,j = (k 1 + 1) * n i,j n i,j + k 1 * (1 -b + b * |di| davg )
where k 1 = 1.2 and b = 0.75 are two constants empirically defined in the Okapi formula, n i,j is the occurrence of the term t j in the document d i , |d i | is the size of the document d i and d avg is the average size of all documents in the corpus. The idf j (inverse document frequency) measures the discriminatory power of a term t j and is defined as [START_REF] Robertson | Okapi at trec-3[END_REF]:

idf j = log |D| -df j + 0.5 df j + 0.5
where |D| is the size of the corpus and df j is the number of documents in which the term t j occurs at least one time.

The weight w i,j is then obtained by multiplying tf i,j and idf j . This weight is higher when the term t j is frequent in the document d i but rare in the others. In our case, the size of our vocabulary T is 217'323 after applying a Porter stemming. The indexing has been performed with the Lemur software1 .

Visual representation model

In order to combine the visual information with the textual one, the images are also represented by a vector of visual words. It is therefore necessary to create a visual vocabulary V {v 1 , ..., v j , ..., v |T | }.

Our method consists in partitioning all images into 16x16 grids, a minimum of 8x8 pixels being required for each cell. It leads to about 256 cells per image, or about 38 million over all images.

For each cell, we compute a feature vector that contains the colour properties of the region. The vector is a 6 dimensional vector. The 6 dimensions correspond to the mean and the standard deviation for We apply a k-means algorithm over 4 millions of cells randomly selected within the 38 millions of cells to obtain 2'000 visual terms, which correspond to our visual vocabulary V . 2'000 has been chosen arbitrarily while 4 millions corresponded to the maximum number of cells we could compute due to the complexity of the k-means. Each visual term represents a cluster of feature vectors.

Then, each new image can be represented using a vector of visual terms. It is decomposed into a 16x16 grid and the local features are computed. Each cell is then assigned to the closest visual term from our visual vocabolary V , using the euclidean distance.

Experiments

Using the model described in the last section, we present our approach for multimedia documents retrieval from multimedia queries. Then we describe the runs we submitted to ImageCLEFwiki in order to evaluate our model.

Queries and matching

As mentionned before, the ImageCLEFwiki topics are composed of text, image and concept modalities. However, our model is designed to take only into account text and image modalities. Our retrieval approach consists in computing a similarity score between each document d i and the query q k using the Okapi method [START_REF] Robertson | Okapi at trec-3[END_REF]. Documents are then ranked according to their scores. The following expression is used to compute the score:

score(q k , d i ) = uj ∈q k tf i,j * idf j * qtw k,j
where qtw k,j is defined as:

qtw k,j = (k 3 + 1) * n k,j k 3 + n k,j
where k 3 = 7 is a constant defined in the okapi formula and where n k,j represents the occurrence of the term u j in the query q k . This score is higher when the term u j is frequent in the document d i but rare in the others and is weighted by the occurrence of u j in the query.

Let us insist on the fact that the term u j can be either a textual term t j or a visual term v j and that queries can be composed of textual terms only, visual terms only or both textual and visual terms which allows to perform text only queries, image only queries or multimedia queries.

The textual terms used for queries are those provided with topics. When visual terms are used, they are extracted either from the topic images or from the collection images. This will be detailed in the following section.

Submitted runs

We have submitted 6 runs to ImageCLEFwiki 2008, labelled from run 01 to run 06. Part of them are automatic (auto), others required manual selection of some relevant documents (man). Thanks to these 6 runs, we aim to study several aspects of our model, as the choice of the visual words and local features, the combination of textual and visual words for a query and the performance improvements obtained when adding visual information to a pure textual model. All these runs are summed up in table 1.

run name first query (R 1 ) run type R 1 use second query (R 2 ) results LaHC run01 t auto - - R 1 LaHC run02 t auto v 10 v t R 2 LaHC run03 t man v 100 v t R 2 LaHC run04 t auto man - v 100 t+ v q if i q exists v t else R 2 LaHC run05 t man v 100 t+v t R 2 LaHC run06 t auto v 10 v t R 1 ∩ R 2
• t: text only query: uj ∈ q ∩ T • R1: first query results (baseline results); R2: second query results

• auto: automatic run; man: selection of relevant images by user

• v10: automatic selection of the first 10 results from R1

• v100: manual selection of the relevant documents in the first 100 results from R1

• iq: query image

• vq: visual words extracted from the query image

• vt: v10 or v100

Table 1: Summary of the runs

We define a baseline that corresponds to a pure text model. It uses only textual terms for the query and scoring of documents. This run is run 01 and its results are noted R 1 . We did not use neither feedback nor query expansion for this automatic run.

All other runs exploit both textual and visual information of documents. They consist in two successive queries. The first one corresponds to the baseline (R 1 ), while the second one is a visual or textual and visual query, either automatic or manual (R 2 ). R 2 is automatic when visual terms are extracted automatically from some images: we chose to select all the visual words of the top 10 retrieved documents issued from the baseline (v 10 ), assuming these results as relevant. R 2 is manual when the user is asked to select relevant documents among the first 100 results of R 1 (v 100 ). There is no limit in the number of selected documents.

Run 02 and run 06 are automatic and only visual runs. For run 02, all the results of the second query are given as answer (R 2 ) while for run 06 an intersection is performed between the results of the first and second query (R 1 ∩ R 2 ). This intersection is interesting as it emphasises the gain of the visual information use. Indeed, if the intersection is not null, it means that the visual query leads us to find documents which were not in the baseline results.

Run 03, run 04 and run 05 are manual runs. Run 03 is visual, while run 04 and run 05 are multimedia. For run 03, all the visual words of the selected images are used to perform the second query. For run 04 and run 05 we did a query expansion in order to analyse the combination of textual and visual information (t+ v 100 ). We kept the textual words of the initial query and added some visual words. For run 05, these words come from the manual selected images. The run 04 proceeds as run 05, unless a query image is provided for the considered topic. In that case, the second query is composed of the visual words extracted from the query image. Thanks to these two runs, we can study the influence of the number of relevant images used for queries.

Results

Table 2 shows that our best textual run (LaHC run01) ranks us 6th on 12 participants. This run is ranked 22th on the whole 77 runs. These results provide a basis to evaluate our approach on different aspects: the choice of the visual words, the way to combine textual and visual words and the performance improvements obtained when adding visual information to a pure textual model.

Rank

Concerning the choice of the visual words, results from run 02 and run 03 show us that they are meaningful in some particular cases. Indeed, using only the visual words, we are able to retrieve relevant documents which were not used for the query. For example, in run 03, we retrieved 42 relevant documents for the query "blue flowers" whereas we had just selected 9 images. However there are some topics for which the results are bad. Only a third of the topics leads to more relevant documents than the number used for the query. This is due to the visualness of the query. It is obvious that for a query like "blue flower", the visual information is more useful than for a query like "peace anti-war protest". Moreover, we observe from the results of run 04 and run 05 that only one image is insufficient to obtain good results. For topics which provide an image query, the results are always worth than the obtained results with several selected images. This can be explained by the fact that image query were not enough representative. Furthermore it is obvious that one image can not be enough expressive compared to several relevant images. The conclusion about the combination of textual and visual words is surprising. The comparison between run 03 and run 05 tells us that adding directly the text to the same visual words leads to worth results. The explanation is not obvious, but it seems to be that the proportion between textual and visual parts is unbalanced. Where we only have 2 or 3 textual words, we have several thousand of visual words.

Results also show that visual information bring complementary relevant documents that were not found with the text query. These initial results are promising and encourage the development of our multimedia model. The comparison between run 01, run 02 and run 06 inform us that 214 new relevant documents were retrieved with the visual query only over all topics. To take an example, 32 documents were found with run 01 and 30 with run 02. The intersection of both results coming from run 06 gives 13 shared documents. Thus, 17 relevant documents were retrieved with the visual information only.

Conclusion

For future work, as our local image features are basic, other features such as texture and edge information would surely lead to performance improvements. We also plan to automatically select the number of visual words using machine learning approaches. Finally, we aim to combine more efficiently textual and visual information.

  3 * 255 where R, G and B are the red, green and blue components of the cell.

Table 2 :

 2 Best textual run of each participantAs we can see on table 3, our best text+image run (LaHC run03) ranks us 4th on 7 participants. This run, ranked 57th on the whole 77 runs, is outperformed by our best textual run (LaHC run01).

		Participant	Run	Feedback/Expansion MAP	P@10
	1	upeking		zzhou3		QE	0.3444 0.4760
	3	ualicante		IRnNoCamel		NOFB	0.2700 0.3893
	4	cea		ceaTxt		QE	0.2632 0.4427
	11	sztaki		bp acad textonly qe		QE	0.2546 0.3720
	13	cwi		cwi lm txt		NOFB	0.2528 0.3427
	22	curien		LaHC run01		NOFB	0.2453 0.3680
	30	chemnitz		cut-txt-a		NOFB	0.2166 0.3440
	44	imperial		SimpleText		NOFB	0.1918 0.3240
	48	irit		SigRunText		NOFB	0.1652 0.2880
	52	ugeneva		unige text baseline		NOFB	0.1440 0.2053
	56	upmc-lip6 TFUSION TFIDF LM		NOFB	0.1193 0.2160
	70	utoulon		LSIS TXT method1		NOFB	0.0399 0.0467
	Rank Participant		Run		Feedback/Expansion MAP	P@10
	10	sztaki		bp acad avg5		NOFB	0.2551 0.3653
	27	chemnitz		cut-mix-qe		QE	0.2195 0.3627
	53	imperial		ImageText		NOFB	0.1225 0.2213
	57	curien		LaHC run03		FB	0.1174 0.2613
	62	upmc-lip6		TIFUSION LMTF COS		NOFB	0.1050 0.2267
	68	upeking		zhou ynli tliu1.1		FBQE	0.0603 0.0040
	71	utoulon	LSIS4-TXTIMGAUTONOFB	NOFB	0.0296 0.0307

Table 3 :

 3 Best text+image run of each participant

http://www.lemurproject.com

LIMA project: http://liris.cnrs.fr/lima/

WI project: http://www.web-intelligence-rhone-alpes.org/

Acknowledgements

This work is supported by the LIMA project 2 and the Web Intelligence project 3 which are 2 Rhône-Alpes region projects.


Marko Horvat
email: marko.horvat@fer.hr

Frane Anzulović

Dominik Kisić

Kristian Jambrošić

THE ACCURACY AND PRECISION OF LOCALISATION OF INSTRUMENTS AND VOCALS FOUND IN A TYPICAL POP OR ROCK BAND IN AN AMBISONICS-BASED VIRTUAL ENVIRONMENT

This paper aims to explore the degree of accuracy and precision of the localization of instruments and vocals that are typically found in contemporary pop or rock bands, when the mix of a recording of such a band is to be made for an Ambisonics-based virtual environment, rather than for the usual two-channel stereo reproduction. For this purpose, a dry recording of a single song was acquired, containing the tracks with drums, bass guitar, lead guitar, keyboards and vocal. A short excerpt of the recording and, consequently, of the corresponding tracks was cut out from the original recording. A virtual environment was created based on Ambisonics encoding. Within this environment, seven positions in the front horizontal half-plane were determined for the placement of the source. Each of the shortened tracks was used as a source signal on all seven positions and the encoding was made in first-, second-and third-order Ambisonics. A series of listening tests was made, in which the subjects were asked to localize the source with its content reproduced to them through a loudspeaker system capable of reproducing Ambisonics signal up to the third order. The results of the tests were analysed in the context of how accurately and precisely the listeners could localize the sources, with regard to source position, the content of the source, and the order of the encoded signal.

INTRODUCTION

Ambisonics is a form of spatial audio (re)production proposed by Gerzon [START_REF] Gerzon | Periphony: With-Height Sound Reproduction[END_REF][START_REF] Gerzon | Practical Periphony: The Reproduction of Full-Sphere Sound[END_REF]. What makes if different from other systems designed for this purpose is that the audio information is processed in two stages; i.e. the encoding stage and the decoding stage. In the encoding stage, a sound environment can be recorded or synthesized through a set of functions known as spherical harmonics. The highest order of spherical harmonics used in the encoding stage determines the final encoding order of the Ambisonics-encoded material. As a rule, higher order encoding should improve localization accuracy and precision. On the other hand, lower order encoding introduces diffusivity into the overall sound image, and the localization accuracy and precision deteriorates, causing a spatial blur [START_REF] Carpentier | Ambisonics Spatial Blur[END_REF]. In the decoding stage, the audio data can be converted and adapted for practically any kind of reproduction system with a reasonably regular configuration, ranging from mono and stereo systems to 5.1, 7.1 or similar setups, all the way to full-3D reproduction. Decoding for binaural listening is available as well. Therefore, Ambisonics offers unprecedented flexibility, compared to all other audio (re)production systems. The highest order of the Ambisonics audio a sound reproduction system will be able to reproduce directly depends on the number of available loudspeakers. With N being the desired Ambisonics order, the minimum number of loudspeakers for 2D reproduction is 2N+1, and for 3D reproduction this number should be (N+1) 2 . To date, many different subjective assessments regarding the performance of Ambisonics systems have been made. In [START_REF] Horsburgh | Subjective Assessments of Higher Order Ambisonic Sound Systems in Varying Acoustical Conditions[END_REF], source perception in higher-order Ambisonics systems is investigated in light of localization accuracy, perceived source width and the overall realistic impression. The stimuli includes speech, pure tones and a musical environment synthesized with multiple instruments. The results speak in favour of third-order systems, regarding their quality, size and localization. In [START_REF] Bertet | Investigation on Localisation Accuracy for First and Higher Order Ambisonics Reproduced Sound Sources[END_REF], localization accuracy was investigated from the viewpoint of different microphones capable of recording first-and higher-order Ambisonics audio. The results show that the localization accuracy depend both on the Ambisonics order and on the direction of arrival of a sound. In [START_REF] Frank | Localization Experiments Using Different 2D Ambisonics Decoders[END_REF], localization has been investigated using Ambisonics decoders of different orders and with different methods of spatial smoothing. The listening test was made both in and out of the sweet spot of the reproduction system. The results show an improvement of localization as the order increases, and better localization in the sweet spot. The method of spatial smoothing also has an influence on localization, with the in-phase decoder as the worst option. This paper presents the investigation of localization accuracy and precision from the possible viewpoint of a composer or a mixing engineer, whose task is to place instruments and vocals into an Ambisonics-encoded auditory scene. In light of this, contemporary music has been chosen as the stimuli for the experiment, rather than hypothetic signals such as pure tones, pulses or broadband noise. Namely, the samples of instruments and vocals typically found in a pop or rock band were used. The purpose of the research is to examine whether the position of real musical sources in a spatial audio mix is perceived correctly, i.e. at the positions they were placed to by the mixing engineer.

The independent variables in the experiment are the type of source, the position of the source, and the Ambisonics order of the encoding used to place the source to the required position. The dependent variable is the perceived azimuthal position of the source, as reported by the listeners.

THE LISTENING TEST

Setup

The listening experiment was conducted in the Auralization laboratory of the Faculty of electrical engineering and computing in Zagreb. The laboratory hosts a spatial audio reproduction system based on Ambisonics coding and decoding scheme.

The listening room, in which the loudspeaker system is mounted, is separated from the control room, where the rest of the hardware is located. A two-way audio connection has been established between the two rooms, consisting of a loudspeaker and a microphone in the listening room (used by the listener), and a microphone and a pair of headphones in the control room (used by the operator). This setup completely prevents the formation of acoustic feedback. At this time, no video communication is installed.

The system consists of 16 loudspeakers in a 4-8-4 configuration. The horizontal plane is covered with eight evenly spaced loudspeakers (45° apart, starting at 0° as the frontal direction). The extension to three dimensions is achieved by placing four loudspeakers above the horizontal plane, namely at the elevation of 40° and the azimuths of 45°, 135°, -45° and -135°. The lower halfspace is covered by another four loudspeakers at the elevation of -40°, at azimuths stated above. All loudspeakers are positioned on a sphere with the radius of 1.65 meters. The described configuration is capable of reproducing third order Ambisonics audio in the horizontal plane, and second order Ambisonics audio in three dimensions.

Due to inherent limitations of the spatial audio reproduction system, the "sweet spot" of the system is small, allowing only one listener at a time to take the test. Therefore, only one fully rotatable barstool is permanently installed in the room. The barstool was chosen instead of a regular chair for its height, so that the listener's head would ultimately be located in the horizontal plane.

The software part of the system consists of REAPER DAW software [START_REF]REAPER software[END_REF] with IEM Plug-in Suite [START_REF]IEM Plugin Suite[END_REF], an open source package dedicated to Ambisonics processing.

The setup of the experiment is shown in Figure 1. Only the horizontal plane is shown, because the virtual sources were placed in the horizontal plane only. The green circles represent the positions of loudspeakers in the horizontal plane. The blue circles show the positions (azimuths) of virtual sources, namely, 90° to the left, then 60° and 30° to the front left, then 0° straight in front of the listener, then -30° and -60° to the front right, and, finally, -90° to the right of the listener. As shown in the figure, the defined setup yields three virtual sources that coincide with actual loudspeakers (90°, 0°, and -90°), and four that do not (60°, 30°, -30°, and -60°).

To facilitate the assessment of the perceived azimuthal direction of the virtual source, a paper tape (in red) was stretched along the circumference of the horizontal circle that contains the loudspeakers. To avoid a skewed distribution of the reported azimuths for the two outermost virtual sources, the tape covered the extended azimuthal range from 120° to -120°. The markings on the tape were made in 1° resolution, and were made to resemble a measuring tape tool; i.e. the 10°-marks were made across the full width of the tape, the 5°-marks were made across the 2/3 rd of the width, and the rest were made across the 1/3 rd of the width of the tape. To avoid confusing the listeners, the values indicated on the tape ranged from 0° at azimuth 120° to 240° at azimuth -120°, thus increasing from left to right and avoiding negative values altogether. The values reported by the listeners were transformed into the standard coordinate system in the data processing stage.

Figure 1.

The positions of the loudspeakers in the horizontal plane of the reproduction system (green), the positions of the virtual sources (blue); the range of azimuths used by the listeners (red)

Stimuli

An excerpt from a raw recording of a rock band was found in a freely available online database [START_REF]Cambridge Music Technology -The 'Mixing Secrets' Free Multitrack Download Library[END_REF]. Originally, it consisted of 12 tracks. Seven tracks covered the drums, and the rest were the bass guitar, two electric guitars, the piano, and the vocal. The instruments and the vocal were recorded one by one, so their recordings were quite clean, with no spill coming from other instruments. To yield the final stimuli that was used in the listening experiment, the drums tracks were mixed together to the final mono downmix, and all the other tracks were taken "as is". Only one electric guitar was selected. In the end, five different mono source samples were obtained: the bass guitar, the drums, the electric guitar, the piano, and the vocal. All samples were cut to 10 seconds, with a short fade-in and fade-out.

To form a test sequence for the listening test, each of the five source samples was encoded as the first, second or third order Ambisonics signal at one of seven defined positions in the horizontal plane. Therefore, the test sequence consisted of 105 encoded samples. Five additional training samples were added to the beginning of the test sequence. The test sequence was formed so that each 10-second sample was followed by 10 seconds of silence, thus yielding the total length of the test sequence of 37 minutes. To avoid systemic errors, the samples were randomly ordered within the test sequence, and three different sequences were made.

Listeners

A total of 20 listeners took part in the listening test, of which 5 were women, and 15 were men. All listeners have a background in audio, sound, and music. The overall average age of the whole group is 32 years. However, it should be mentioned that 11 listeners were students in the age range from 20 to 23 years old, while the rest of the group consisted of faculty employees, most of which are in their late 30s and 40s.

Initially, all listeners reported to have normal hearing. To test their clams, the responses of each listener were tested by assessing the histogram of all the reported values. The values were first transformed to yield the deviation of the reported azimuth from the real azimuth of a virtual source. For a person with normal hearing, the expectation was to obtain a symmetric distribution around zero, as shown in Figure 2.

Figure 2.

The distribution of the deviation of the reported azimuth from the real azimuth of the virtual source for a listener with normal hearing However, during this analysis, the histogram was obtained for one of the listeners as shown in Figure 3, suggesting that the perception of direction for this listener is consistently shifted to the right by an average of 12°. The listener ultimately reported to have tinnitus in their left ear, and their test results were discarded.

Figure 3.

The distribution of the deviation of the reported azimuth from the real azimuth of the virtual source for a listener with impaired hearing

Procedure

The listeners were tested individually. At the beginning of the test, the listeners were informed on the nature of the experiment. It was asked of them to try to locate the presented sound sources as accurately as possible, using the tape measure with marked angles.

They were advised to close their eyes while performing the task and to point their hand in the desired direction. This way, all the distractions in the listening room were avoided. They were allowed to rotate on the barstool. Once they were sure that the chosen direction is the right one, they could open their eyes and simply read the angle value on the tape in the direction in which their hand was pointing. To allow the listeners to concentrate solely on their task, they were not required to record the values by themselves in any way, e.g. by writing them down, thus leaving their hands free. Instead, they were asked to use the audio communication system by uttering the values during the 10-second period of silence between two sound samples. The operator in the control room constantly monitored the open audio link to the listening room, and repeated the values back to the listener. If the read-back value was the equal to the one originally uttered by the listener, the operator wrote it down into the test form. If not, the listener warned and corrected the operator, and the correct value was written down.

Before engaging in the test, five training samples were presented to the listeners, thus familiarizing them with the test procedure. Simultaneously, the listeners were acquainted with all five types of instruments/vocals that would be repeatedly used in the rest of the test.

The entire test was performed with no breaks. Afterwards, the listeners gave the feedback to the operator, claiming without exception that the test was not physically or mentally demanding.

RESULTS AND DISCUSSION

As 20 listeners evaluated 105 different sound samples (5 instruments x 7 source positions x 3 Ambisonics encoding orders), a total of 2100 observations was collected.

The raw observations were azimuthal angles in the interval [0°, 240°], increasing clockwise. To match the standard coordinate system (x-axis points to the front, yaxis points to the left), the data was transformed to match the original azimuth range of [-120°, 120°], that increases in the counterclockwise direction.

The data was further transformed so that is reflects the deviation of the perceived azimuth from the assigned azimuth of a source. A negative value of the deviation implies that the listener perceives the sound source to the left of its "true" position, whereas a positive value of the deviation implies that the perceived position of a source is to the right of its assigned position.

The defined dependent variable that represents the described deviation was used to test both the accuracy and the precision of the localization of the sound sources, as displayed by the listeners. The mean value was taken as the measure of accuracy, and the variance as the measure of spread was taken as the descriptor of precision. Ideally, if all the listeners were able to locate the sound source perfectly, the mean value of the deviation would be zero, and the variance would be equal to zero as well.

Normality of the data

In the initial stage, the data was checked for normality by means of the Lilliefors normality test. Not only was each of the 105 test cases checked, but also the data grouped by instrument (5 groups), azimuth (7 groups), and order (3 groups). The normality tests revealed that the observations are normally distributed in about 65% of the tested cases and groups. The ones that did not pass the test were inspected visually by means of histograms, and it was determined that there are no severe violations of the normal distribution.

Accuracy of localization

In the second stage, the localization accuracy was tested by means of a one-sampled t-test. First the two-tailed test was made, with the null hypothesis was that the sample mean of the deviation for each of the 105 cases was not significantly different from zero, against the alternative hypothesis that it is. The results of this test reveal that the null hypothesis can be rejected at the 0.05 significance level in 47 out of 105 cases. Next, both lower tail and upper tail one-tailed tests were made as well. The null hypothesis was the same as before. The alternative hypothesis for the lower tail test was that the sample mean is lower than zero, indicating a significant perceived shift of the source to the left of its true position. Similarly, the alternative hypothesis for the upper tail test was that the sample mean is larger than zero, indicating a significant perceived shift of the source to the right of its true position. To match the two-tailed test, these tests were made at the 0.025 significance level. The results show a statistically significant shift to the left in 14 cases, and to the right in 33 cases. The count of significant deviations from zero, made according to instrument, is the following: 8 cases for bass and vocal, 10 cases for drums and piano, and 11 cases for the guitar. The same count made taking into account the direction of arrival reveals the following: one case for 90°, 11 cases for 60°, 15 cases for 30°, 9 cases for 0° degrees, 2 cases for -30°, 6 cases for -60°, and 3 cases for -90°. The count made according to the encoding order is the following: 12 cases for first order, 17 cases for second order, and 18 cases for third order.

The results indicate that the accuracy of localization does not depend much on the type of sound, i.e. the instrument/vocal. On the other hand, the direction of arrival of sound, as well as the encoding order seem to have an influence on the localization accuracy.

The surprising finding is that the first order Ambisonics encoding is apparently the most accurate one, and that the source placed right in front of the listener, i.e. at 0°, will significantly deviate from that direction, as perceived by the listener. Additionally, all cases with sources placed at 30° show a significant deviation to the right from that direction, while most of the cases with sources at 60° indicate a significant deviation to the left.

To explain these findings, the accuracy indicated by the means has to be observed in light of precision indicated by the standard deviations. Specifically, if the spread of data represented by standard deviation is large, even a large deviation from the true direction might turn out not to be statistically significant. On the other hand, a small spread of data indicated by a small standard deviation will make even a rather small deviation from the true direction statistically significant.

To accompany the results of the t-test, the sample means for all 105 cases are shown in Table 1, and the corresponding standard deviations are shown in Table 2.

The data confirms that the frontal direction generally has the smallest spread of observations, yielding a number of cases with small mean deviation declared as statistically significant. The comparison of standard deviations by encoding order reveals that the first order is indeed the least precise, having the largest spread of observations. As stated above, this can lead to a misleading conclusion that it is the most accurate.

Given the fact that the sources placed at 30° subjectively significantly deviate from that direction in all cases leads to the conclusion that there is a possibility that some kind of systemic error has occurred. Detailed examinations of the listening setup have found nothing that supports this conclusion. Table 1. The means of the deviation from true direction of arrival; values in bold indicate statistically significant cases.

Instrument

Precision of localization

In the third stage of statistical analysis, attention was turned to evaluating the precision of localization. To achieve this, the hypothesis of equal variances was tested by means of the Levene's test for equality of variances. This test is reasonably tolerant to deviations from normality and to the existence of outliers. The null hypothesis is that the variances of all tested groups are equal. In this particular case, the null hypothesis states that the independent variable (instrument, direction, or encoding order) has no influence on the precision of localization indicated by the variance of the data. The alternative hypothesis is that the variances are not equal, i.e. that at least one group is significantly different in this respect from the other ones. The Levene's test was performed on the data grouped by instrument, direction, and encoding order. The results of this test are shown graphically in Figures 4, 5 and 1. The Levene's test on the data grouped by instrument did not prove to be statistically significant at the 0.05 level, indicating that the null hypothesis that the variances of these five groups are equal cannot be rejected. Therefore, the conclusion of the test is that the type of intrument does not have an influence on the precision of its localization in the given setup. The boxplot shown in Figure 4 supports this conclusion. The Levene's test on the data grouped by direction was statistically significant at the 0.05 level, indicating that the assigned direction has a strong influence on the precision of localization. The boxplot shown in Figure 5 reveals that all three assigned directions that coincide with the directions of actual loudspeakers (90°, 0° and -90°) have a much smaller spread of observations than the remaining four directions (60°, 30°, -30° and -60°). This finding confirms the theory that lies behind Ambisonics systems and their predicted behaviour. The Levene's test on the data grouped by encoding order was statistically significant at the 0.05 level, indicating that the encoding oder has an influence on the precision of localization. The boxplot shown in Figure 6 the spread of observations is smaller for a higher encoding order, indicating that the precision of localization indeed increases with increasing encoding order. This finding is also in agreement with the behaviour of Ambisonics-based systems.

Limitations

The limitations of the study are mainly connected to the data collection method and to the limitations imposed by the laboratory setup. The listeners were asked to close their eyes and to open them again only after they had chosen the perceived direction of the source. Nevertheless, they often reported that them being able to see the loudspeakers was a bit distracting. Therefore, to improve similar estimations in the future, a cylinder of suffcicient size, made of acoustically transparent fabric, shall be placed around the listener position. In this manner, a visual barrier will be created that will prevent the listeners to see any part of the sound reproduction system. The data collection method that utilizes the paper tape with indicated angles can only be used in the horizontal plane, with no possibility of extending it to all three dimensions by assessing the elevation as well. Therefore, a more adequate method needs to be devised before any localization experiments are to be made in full 3D space. Finally, many listeners reported the sound coming from above the horizontal plane, even though the position of the sources was set to horizontal plane only. As it was beyond the scope of this study to study this effect, and the laboratory setup was not equipped for assessing the perceived elevation of the source, the listeners were advised to report only the perceived azimuth.

CONCLUSIONS

The research presented in this paper was focused on examining whether or not accuracy and precision of source localization in an Ambisonics-based mix is affected by the type of the source, its assigned direction and the Ambisonics order used in the encoding stage. The hypothetical mix is made of sounds of instruments and vocals typically found in a rock band, and the sources are placed in the front half of the horizontal plane.

The results reveal that the bass guitar and the vocal are the sounds that are the easiest to locate accurately, followed by the drums and the piano, while the electric guitar is the one most difficult to place to its intended position.

The analysis of accuracy dependent of the true position of the source yielded a surprising result: the sounds coming from the loudspeakers directly to the left and to the right of the listener were the one that were located the most accurately. In general, the sounds coming from the front right quadrant were located accurately, while the sounds coming from the left quadrant showed considerable deviation of their perceived direction from the assigned one. Moreover, the sounds coming from the direction 30° had the tendency of deviating to the right, i.e. towards the frontal direction of 0°, whereas the sound coming from the direction 60° had the tendency of deviating to the left, i.e. towards the left side direction of 90°.

The accuracy dependent on the encoding order yields another surprising finding that the first order seems to yield the most accurate localization.

The precision of localization does not yield a statistically significant difference between the instruments involved in the experiment.

Figure 4 .

 4 Figure 4. The boxplot of the data grouped by instrument to complement the corresponding Levene's test

Figure 5 .

 5 Figure 5. The boxplot of the data grouped by direction to complement the corresponding Levene's test

Figure 6 .

 6 Figure 6. The boxplot of the data grouped by encoding order to complement the corresponding Levene's test

Table 2 .

 2 6. The standard deviations of the deviation from the true direction of arrival; values in bold indicate statistically significant cases to accompany the means in Table

	Instrument

10.48465/fa.2020.1136

e-Forum Acusticum, December 7-11, 2020

ACKNOWLEDGMENTS

The authors would like to thank the creators of the IEM Plugin Suite as an invaluable resource and tool, and to express their gratitude to Mike Senior of Cambridge Music Technology for creating the vast database of raw audio material and making it freely available.

The influence of the assigned direction of the source on the precision of localization is clearly visible. The source positions that coincide with actual loudspeakers were located with significantly better precision than those placed in directions that do not coincide with any loudspeakers.

The encoding order has a significant influence on the precision of localization, in the sense that the precision improves with increasing encoding order.

Future research will strive towards extending the possible directions of arrival to all three dimensions. For such an experiment, a reliable way of assessing both the azimuth and the elevation needs to be devised.


Paweł Łojek
email: lojek@agh.edu.pl

Andrzej Gołaś

LINEARIZED EULER EQUATIONS FOR TIME ANALYSIS OF WAVE PROPAGATION IN DUCTS

This paper presents the use of linearized Euler equations to solve acoustic wave propagation problems in ducts and mufflers with and without the presence of mean flow. The software developed by the authors is used to solve the equations using the finite element method. In order to be able to model a wider range of cases, a non-reflective boundary condition -the buffer zone -was also implemented. In addition, the tool allows to take the mean nonuniform flow velocity field into account. The tool was developed using open-source FEniCS library and Python programming language. In this paper, developed software is used to solve the acoustic wave propagation for 2 different reflective mufflers. The results were compared to the analytical model of insertion loss of the mufflers. The tool can be used to estimate the insertion and transmission loss of newly developed mufflers with complex geometry.

INTRODUCTION

Acoustic mufflers are used in many branches of industry to reduce ubiquitous noise. Examples include exhaust and ventilation noise. Behavior of reflective mufflers is usualy computed using one-dimensional wave equation,or estimated by means of experimental research [START_REF] Gołaś | Acoustics of Ducts and Mufflers[END_REF]. Because of the simplification of the equations or limited research capabilities, it does not give a complete description of the phenomena occurring in the mufflers. In this study, we attempt to describe the phenomena of acoustic wave propagation and dampening in mufflers by means of linearized Euler equations.

Linearized Euler equations are one of the methods to describe sound propagation in time domain. They are related to wave equation, but contrary to it, they describe the propagation in terms of acoustic velocity and pressure. This property of linearized Euler equations can be used to couple the acoustic field to flow field and allows to study the influence of flow on wave propagation.

In this study, we present the tool developed using FEn-iCS library, allowing for solving linearized Euler equations, as well as the results of wave propagation in mufflers obtained with it. The FEniCS is finite-element computing platform based on Galerkin method and is decribed in detail in [START_REF]Automated Solution of Differential Equations by the Finite Element Method[END_REF]. This library allows for quick implementation of partial differential equations written in variational form and for solving them by finite element method. Developed solver is described in greater detail in [START_REF] Łojek | Finite element modeling of a flow-acoustic coupling in unbounded domains[END_REF].

We conducted simulations of wave propagation in mufflers for various frequencies. For one of the mufflers, we additionally investigated the influence of the mean flow velocity on wave propagation.

Modeling of the free outflow from the muffler was possible by implementation of a non-reflecting boundary condition, called buffer zone, which significantly increased the capabilities of the tool.

The developed code is easily scalable, and can be used both on personal computer and HPC cluster. All of simulations described in this work was conducted on Prometheus HPC cluster.

THE SOLVER

Mathematical model

To describe the propagation of the acoustic wave, the linearized Euler Equations were used. They were chosen instead of the wave equation due to the possibility of coupling them with the mean flow. The equations can be derived from Navier-Stokes equations by assuming inviscid, homoentropic flow [START_REF] Mechel | Formulas of Acoustics[END_REF] and by linearizing the flow around mean values [START_REF] Łojek | Finite element modeling of a flow-acoustic coupling in unbounded domains[END_REF] . The most basic form of linearized Euler equations (without mean flow) is given by [START_REF] Åbom | An introduction to Flow Acoustics[END_REF]:

∂p ∂t + c 2 ρ ∂v j ∂x j = 0, (1)
ρ ∂v i ∂t + ∂p ∂x i = 0 (2)
p = p + p , v i = v i +v i , ρ = ρ + ρ , (3)
where: p, v, ρ -mean variables, p , v i , ρ -acoustic variables (pressure, velocity in i-th direction, density).

The equations (1,2) do not take mean flow into account. If the presence of the flow was assumed, the equations would take more complex form. For two dimensional case, it can be written [START_REF] Povitsky | Numerical study of wave propagation in a non-uniform flow[END_REF]:

∂p ∂t = -ρc 2 ∂u ∂x + ∂v ∂y -U ∂p ∂x -V ∂p ∂y -p ∂U ∂x + ∂V ∂y (4) ∂u ∂t = - 1 ρ ∂p ∂x -U ∂u ∂x -u ∂U ∂x -V ∂u ∂y -v ∂U ∂y - p ρc 2 U ∂U ∂x + V ∂U ∂y (5)
∂v ∂t = - 1 ρ ∂p ∂y -V ∂v ∂y -v ∂V ∂y -U ∂v ∂x -u ∂V ∂x - p ρc 2 U ∂V ∂x + V ∂V ∂y (6)
where:

v 1 = u, v 2 = v, v 1 = U, v 2 = V, p = p, p = 0, x 1 = x, x 2 = y.
Both versions of linearized Euler equations, with (45,6) and without the flow (1,2) was implemented in the solver.

Non-reflective boundary condition

In order to model the propagation in unbounded domains, or domains with inlet/outlet, non-reflective boundary condition (NRBC) was also implemented. Among many NR-BCs, the buffer zone model, similar to the ones described in [START_REF] Richards | The evaluation of non-reflecting boundary conditions for duct acoustic computation[END_REF][START_REF] Gill | Towards an effective non-reflective boundary condition for computational aeroacoustics[END_REF] was chosen.

Buffer zone was introduced into the governing equations by adding a damping term as follows:

∂φ ∂t = L(φ) + σ(x i )φ (7)
where φ -flow variable, L(φ) -spatial operator, σ(x i )damping term.

The damping term can be described as piecewise function:

σ(x) = -α 1 -x L β x > L, 0, x < L. (8)
where: α -strength of damping, L -start of the buffer zone, β -shape of the damping curve.

Flow modeling

Linearized Euler Equations can take the mean flow into account. Navier-Stokes equations were used to model the flow inside the domain and to compute mean values of velocity. Because we were interested only in mean values, the incompressible, time-independent Navier-Stokes equations were used [START_REF] Hirsch | Numerical Computation of Internal and External Flows: The Fundamentals of Computational Fluid Dynamics[END_REF]:

ρ ∂u j u i ∂x j -ν ∂u j ∂x j + ∂p ∂x i = 0 (9)
where: ν -kinematic viscosity.

Numerical methods

The linearized Euler Equations with damping term had to be discretized, both in time and space.

Adams-Moulton method of the first order [START_REF] Czajka | Engineering methods of numerical analysis and design of experiment[END_REF], given by eq [START_REF] Czajka | Engineering methods of numerical analysis and design of experiment[END_REF], was used for temporal discetization. For the equation in form ∂x/∂t = f (x, t), one can write

x [n+1] = x [n] + ∆tf x [n+1] , t n+1 (10)
where:

x (t n ) -x [n] , x (t n+1 ) -x [n+1] , ∆t = t n+1 -t n
The FEniCS software was used for spatial discretization and solving the equations. It requires writing equations in variational (weak) form in order to solve them. For linearized Euler Equations without mean flow [START_REF] Gołaś | Acoustics of Ducts and Mufflers[END_REF][START_REF]Automated Solution of Differential Equations by the Finite Element Method[END_REF] discretized in time, it can be written:

Ω p [n+1] q dx - Ω p [n] q dx + ∆tc 2 ρ Ω ∇ • u [n+1] i q dx + σ(x) Ω p [n+1] q dx = 0 in Ω, (11)
Ω v [n+1] i w dx - Ω v [n] i w dx + ∆t ρ Ω ∇p [n+1] w dx + σ(x) Ω u [n+1] i w dx = 0 in Ω. ( 12
)
The equations that take mean flow into account were not shown here, because of their complexity, however they have been implemented and used in the simulations.

The Navier-Stokes equation given by ( 9) also had to be written in variational formulation [START_REF] Blechta | FEniCS hands[END_REF]:

Ω ∂u j ∂x j u j w j dx + Ω ν ∂u j ∂x j ∂w j ∂x j dx - Ω p ∂w ∂x i dx - Ω ∂u j ∂x i qdx = 0 (13)

The algorithm

The entire simulation process can be divided into several stages:

• mesh preparation using GMSH software,

• mesh conversion from .msh to .xml format used by FEniCS,

• flow simulations with (13) , export of mean velocity results (optional step),

• acoustics simulations with [START_REF] Blechta | FEniCS hands[END_REF][START_REF] Sreenath | Evaluation of noise attenuation due to exhaust mufflers[END_REF], with or without mean flow,

• analysis of result with Paraview and Python programming language.

Both acoustic and flow simulations are performed on the same mesh, to ensure compatibility.

ANALYTICAL MODEL

In this study, two models of the mufflers was analysed, i.e. two-chamber and side-inlet side-outlet muffler. These two

d 1 d 2 l 1 l 3 l 1 l 2 l 2
IL = 20 log Z n+1 Z n+1 + Z 0 v n+1 v 0 = 20 log p n p 0 ( 14
)
For the two-chamber muffler, shown in fig. 1, the algorithm derived by Munjal [START_REF] Munjal | Velocity ratio in the analysis of linear dynamical systems[END_REF][START_REF] Munjal | An algebraic algorithm for the design and analysis of linear dynamical systems[END_REF] was used to calculate velocity ratio v n /v 0 and thus the insertion loss. Calculated insertion loss was shown in the fig. 4.

For the SISO muffler, shown in the figure 3, the acoustic equivalence of SISO and extended-inlet/extended-outlet (EIEO) was used [START_REF]Plane wave analysis of side inlet/outlet chamber mufflers with mean flow[END_REF][START_REF]Tuning a two-chamber muffler for wideband transmission loss[END_REF]. The plane wave equations for EIEO was used to derive the insertion loss equation.

Calculated insertion loss was shown in the fig. 4. Computational mesh was based on models shown on fig.

1 and 3 with dimensions as in table 1. Two-dimensional mesh was prepared and created with GMSH, open-source finite element mesh generator. For both cases, the mesh consisted of triangular elements, created with Frontal-Delaunay algorithm. Two-chamber muffler mesh consisted of 75368 nodes and 150754 elements, while SISO muffler mesh -of 54729 nodes 109468 elements. Fragment of SISO muffler mesh was shown in fig. 5.

Boundary Conditions

For both flow and acoustic simulations, three named boundaries was defined, inlet, outlet and wall. The inlet and outlet boundaries was marked in the figs 6 and 7. On all unmarked boundaries, the wall boundary condition was assumed. s ) The boundary conditions were defined separately for pressure and velocity. Pressure Neumann boundary conditions were defined on each boundary:

Flow simulations

∂p ∂n = 0 (15)
For velocity field, Neumann boundary condition was defined for outlet:

∂u i ∂n = 0 (16)
and Dirichlet boundary conditions were defined for walls:

u i = 0 (17)
Dirichlet boundary condition was also defined for inlet, as parabolic velocity distribution. For SISO muffler, velocity was defined as:

u 1 = 0, u 2 = 400x 2 1 -90x 1 + 3, (18)

Acoustics simulations

For acoustic simulations, the medium in which wave propagates, was also air (density ρ = 1.2 kg m 3 , speed of sound c 0 = 340 m s ). The simulated time in each simulation was equal to t = 0.02s, and the time step was equal to ∆t = 1e-5. For inlet, pressure was defined as Dirichlet boundary condition given by sine function:

p = A sin(2πf t) (20)
where: A = 1 -amplitude, f -constant frequency, ttime. A series of simulations were carried out for two dampers, in which the frequency of pressure changed, it assumed the following values: 340, 500, 1000. For velocity field, Neumann boundary condition was defined for inlet and outlet

∂u i ∂n = 0 (21)
and Dirichlet boundary conditions were defined for walls:

u i = 0 (22)
In acoustics simulations, to model anechoic termination, the buffer zone was used. For both mufflers, the coefficents of eq. ( 8), was equal to α = 25, β = 3. The zone was 0.2m long from the outlet. asd 5. RESULTS

Flow simulations

Flow simulations were conducted with FEniCS software, using MUMPS solver (MUltifrontal Massively Parallel sparse direct Solver) and Newton method. For simulations, absolute convergence criterion was 1•10 -4 and relative one was equal to 1 • 10 -5 . Flow through SISO muffler converged after 180 iterations. The velocity field distribution is shown in the fig. 8.

f [Hz] T I T C A [dB] T I T C S [

Acoustic results

The acoustic simulations were conducted with FEniCS, with GMRES (generalized minimal residual method) method and ILU (incomplete LU factorization) preconditionar.

Velocity and pressure fields shown in the following sections, was generated for time of simulation equal to t = 0.02s.

Results in tables 2 and 3 were calculated by Fourier analysis of the signal recorded in the intake and exhaust pipes.

Two-chamber muffler

The results for two-chamber muffler was shown in tab. 2, where A subscript stands for analytical, and S for simulation insertion loss of two-chamber muffler.

The pressure distributions for various frequencies were shown on figs 9 and 10. The velocity distributions for various frequencies were shown on figs 11 and 12. Also, the spectrum of pressure measured in inlet/outlet pipe was compared and shown on fig. 13.

Side-inlet/side-outlet muffler

The results for SISO muffler was shown in tab.

DISCUSSION AND CONCLUSIONS

In this study, we we presented the capabilities of the developed solver. It was used to calculate the wave propagation in the mufflers without and in the presence of flow. Initial tests of the solver, described in [START_REF] Łojek | Finite element modeling of a flow-acoustic coupling in unbounded domains[END_REF], showed that it shows good agreement with the analytical formulas. The comparison of the simulation results with the results obtained from the analytical formulas shows that this tool allows to obtain results similar to the analytical ones.

The discrepancies in this study may be due to simplifications assumed by the analytical formulas or the fact that a two-dimensional system was modeled. That does not diminish it's value, developed tool still can be used to model various acoustic systems.

Mean flow field obtained with Navier-Stokes equations were used as input data for the flow-acoustic coupling simulations. The developed tool allows for the analysis of the impact of mean velocity on wave propagation, but it must be noted, that in current form, it does not take acoustic sources into account.

The use of a FEniCS software allowed to focus on issues related to the phenomena itself, instead of developing numerical methods code allowing for pre-and postprocessing and for the solution of matrices resulting from the discretization of the equations. Currently there are very few tools, both paid and open source, that allow solving linearized Euler equations, also in the presence of average flow.

The tool is constantly being developed, currently work is underway on the possibility of coupling it with Open-FOAM software. This will allow for a more accurate analysis of flows than the relatively simple implementation of

ACKNOWLEDGMENTS

This research was supported in part by PLGrid Infrastructure.

Figure 1 .Figure 2 .

 12 Figure 1. Schematic of two-chamber muffler

1 Figure 3 .

 13 Figure 3. Schematic of side-inlet side-outlet muffler

Figure 4 .

 4 Figure 4. Insertion loss for EIEO muffler (dimensions as in 1)

Figure 5 .Figure 6 .Figure 7 .

 567 Figure 5. Part of computational mesh of two-chamber muffler

For

  fluid flow simulations, flow of the air in 20 • C was modeled (density ρ = 1.2 kg m 3 , kinematic viscosity ν = 2 • 10 -5 m 2

Figure 8 .

 8 Figure 8. Mean velocity distribution for flow inside SISO muffler

Figure 9 .

 9 Figure 9. Pressure field of wave propagation in twochamber muffler for inlet pressure frequency f = 340Hz

Figure 10 .

 10 Figure 10. Pressure field of wave propagation in twochamber muffler for inlet pressure frequency f = 1000Hz

Figure 11 .Figure 12 .Figure 13 .Table 3 .Figure 14 .Figure 15 .Figure 16 .

 1112133141516 Figure 11. Velocity field of wave propagation in twochamber muffler for inlet pressure frequency f = 340Hz

Figure 17 .

 17 Figure 17. Velocity field of wave propagation in twochamber muffler for inlet pressure frequency f = 1000Hz

Figure 18 .

 18 Figure 18. SPL measured at inlet (dashed line) and outlet, without (solid line) and with (dotted line) for different inlet pressure frequencies for SISO muffler

Table 1 .

 1 Mufflers

dimensions (in meters) Muffler dimensions were shown in table 1

4. SETTINGS OF THE SIMULATIONS 4.1 Geometry

Table 2 .

 2 Comparsion of analytical and simulation results for two-chamber muffler

				dB] Error [dB]
	340	-0.4	0.36	0.26
	500	37	14.58	22.42
	1000	11	9.6	1.4

10.48465/fa.2020.1144

e-Forum Acusticum, December 7-11, 2020