A statistical theory for the measurement and estimation of Rayleigh fading channel

In this paper, we propose a statistical theory on measurement and estimation of Rayleigh fading channels in wireless communications and provide complete solutions to the fundamental problems: What is the optimum estimator for the statistical parameters associated with the Rayleigh fading channel, and how many measurements are sufficient to estimate these parameters with the prescribed margin of error and confidence level? Our proposed statistical theory suggests that two testing signals of different strength be used. The maximum likelihood (ML) estimator is obtained for estimation of the statistical parameters of the Rayleigh fading channel that is both sufficient and complete statistic. Moreover, the ML estimator is the minimum variance (MV) estimator that in fact achieves the Cramer-Rao lower bound.

Introduction
In mobile radio channels, the Rayleigh distribution is commonly used to describe the statistical nature of the received envelope of a flat fading signal, or the envelope of an individual multipath component. Flat fading is often associated with the narrow band channel. By assuming that the real and imaginary parts of the channel gain are independent Gaussian random variables with zero mean and equal variance, the amplitude of the channel gain or PL (path loss) becomes a Rayleigh random variable. For wide band channels, multipath gains are typically assumed to be uncorrelated scattering (US) and each path gain is wide-sense stationary (WSS) which are termed as WSSUS channels [2]. By assuming that the real and imaginary parts of each path gain are independently Gaussian distributed with zero mean and equal variance, the amplitudes of the multipath gains become independent Rayleigh random variables. Theoretically a Rayleigh random variable is uniquely specified by its 2nd moment that is the sum of the variances of its two independent Gaussian components. In practice the statistics of the Rayleigh fading channel are incomplete without additional knowledge of the noise power and SNR (signal-to-noise ratio). That is, the statistical parameters of the channel power, noise power, and SNR together characterize the Rayleigh fading channel completely, giving rise to the measurement and estimation problem for these statistical parameters.
Estimation for the 2nd moment of a Rayleigh random variable has practical importance in channel modeling and estimation [3,9,12], and in radio coverage, location, and measurement [4,8]. For these reasons statistical estimation of Rayleigh fading channels has been studied and reported in the existing research literature. The early work of Peritsky [13] shows that the simple averaged square of the signal strength based on i.i.d. (independent and identically distributed) measurement samples is both an ML (maximum likelihood) and MV (minimum variance) estimator, and such an averaged square is a sufficient and complete statistic. However his results focus only on the noise-free case and have limited applications. The same problem has been investigated in [1,6,7,10,11,19,20] that encompass the Ricean and Nakagami fading distributions, as well as MIMO fading channels. But there lack optimum estimators that achieve the Cramér-Rao lower bound.
In this paper we propose a statistical theory on measurement and estimation of Rayleigh fading channels. Our contributions include derivation of explicit a priori bounds on the measurement sample size that achieve the prescribed margin of error and the confidence level, and discovery of the optimum estimator that is both an ML and MV estimator, and achieves the Cramér-Rao lower bound. These results complement the existing work reported in the literature, In the noise-free case our sample size bound resembles the one implicitly indicated in [5] which asserts that, to estimate the binomial probability with the prescribed margin of absolute error ε and confidence level 1 − δ, it suffices to have a sample size greater than ln( 2 δ ) 2ε 2 . One notable difference is that our sample size bound is derived for the relative error while the bound in [5] is derived for the absolute error in estimation of the binomial probability. Furthermore an interval estimate similar to that in [13] is obtained with much simpler calculation.
In applications to Rayleigh fading channels, noisy measurement samples have to be taken into consideration. Our proposed statistical theory suggests that two testing signals of different strength be used. The use of two different testing signals enables us to extend the sample complexity results from the noiseless case to the noisy case which solves the sample complexity problem not only for the 2nd moment or the mean channel power, but also for the mean noise power and for the SNR. The explicit a priori bounds resemble to those in the noise-free case of which the sample size is roughly inversely proportional to the square of the margin of error and is linear with respect to the logarithm of the inverse of the gap between the confidence level and 1. Our results show that a typical margin of error can be achieved with near certainty and modest sample size. More importantly an optimum estimator is obtained for noisy Rayleigh fading channels that is both an ML and MV estimator. It inherits the same sufficient and complete statistic property from that in the noiseless case [13], and it achieves the Cramér-Rao lower bound.
The content of the paper is organized as follows. In Section 2 we first present the Rayleigh fading channel and its associated measurement and estimation problem in terms of the sample complexity and the optimum estimator. The sample complexity problem is then addressed in Section 3 by establishing an explicit a priori bound on the sample size that is asymptotically tight based on noiseless measurements. In Section 4 we present our main results on parameter estimation for Rayleigh fading channels by taking the noisy measurement samples into consideration to which complete solutions are derived for measurement and estimation of Rayleigh fading channels. Numerical examples are presented in Section 5 to illustrate our results, and the paper is concluded in Section 6. The notations are standard and will be made clear as we proceed.

Rayleigh Fading Channels
In a typical urban environment, there is no LOS (line of sight) between the transmitter and the receiver. The wireless channel is characterized by multipath and follows the Rayleigh distribution. Specifically a widely used channel model is the following continuous-time CIR (channel impulse response): where α n (t), β n (t), and τ n (t) are the amplitude, phase, and time delay, respectively, associated with the nth path, and δ(·) is the Dirac delta function. For narrow band systems, the channel gain is given by The arguments {φ n (t)} are due to {β n (t)} and {τ n (t)} induced by Doppler shifts and time delays that form a set of independent random variables uniformly distributed over [−π, π]. Under the assumption that α n (t) ≡ α n , β n (t) ≡ β n , and τ n (t) ≡ τ n , H(t) becomes a WSS (wide-sense stationary) process. For large N , the central limiting theorem can be invoked to treat both H R (t) and H I (t) as independent Gaussian random processes with zero mean and equal variance [14,18]. This gives rise to the Rayleigh fading channel in light of the fact that |H(t)| 2 is a Rayleigh random variable. For wide band systems, multipath gains can be resolved up to the resolution dictated by the sampling frequency f s . An equivalent discretized CIR can be represented by in which |H k (t)| 2 is Rayleigh distributed for 0 ≤ k < L where δ t is the Kronekar delta function.
Rayleigh fading channels are widely used in wireless communications. An important and practical problem is measurement of |H(t)| 2 and estimation of E[|H(t)| 2 ], with E[ · ] the expectation, for narrow band systems, or |H k (t)| 2 and power delay profile E[|H k (t)| 2 ] for wide band systems. Such a problem is crucial in modeling and estimation of wireless channels, and in radio coverage and location. A common strategy is to transmit i.i.d. pilot symbols or testing signals to obtain measurement samples of |H(t)| 2 or |H k (t)| 2 , and then to estimate the mean. For simplicity we drop the time argument t, and focus on narrow band systems but our results are applicable to wide band systems as well by converting frequency-selective fading channels into flat fading channels using the OFDM (orthogonal frequency division multiplexing) scheme [12]. The detail is omitted.
Denote ℜ[C] as the real part of the complex number C. Let the transmitted symbol be S. In the complex form, the received signal can be written as ("=⇒" stands for "implying") where f c is the carrier frequency and V = V R + jV I is the complex additive Gaussian noise with zero mean and equal variance σ 2 V . Since the average power of the transmitted signal is P s = E[|S| 2 ] and the average power of the received signal is E[|HS| 2 ], the channel power or PL is where σ 2 H is the common variance of H R and H I . It can be seen that σ 2 H depends on the attenuation factors {α n } that decrease as the receiver is further away from the transmitter. It follows that σ 2 H and thus SNR become small when the receiver is far away from the transmitter. To predict the coverage and the performance of wireless systems, it is important to know the operating range of the SNR at different locations. The SNR at the receiver is easily found to be where σ 2 V is the common variance of V R and V I . Since P s = E[|S| 2 ], the average power for transmitting the symbol S, is computable, to obtain an estimate for the SNR it is necessary to have an estimate for the ratio of σ 2 H to σ 2 V . However it remains unclear whether or not the knowledge of SNR helps to estimate the PL that will be investigated later. Our main results to be presented show that with appropriate signaling, the variances σ 2 H , σ 2 V , and the SNR can all be estimated in satisfaction.
The objective of this paper is in fact beyond estimation of the PL that has been investigated by several researchers [7,10,11,13,19,20]. Our goal is to develop a relevant statistical theory in order to derive an optimum estimator for the statistical parameters associated with the Rayleigh fading channel, and to derive a priori sample complexity bounds given the prescribed margin of error and the confidence level. These two problems are not completely solved in the existing research literature that will be studied thoroughly in the next several sections. We will begin with the ideal case of noiseless measurements by assuming independent PSK (phase shift keying) or BPSK (binary PSK) signaling at the transmitter leading to N independent measurement samples {X i } N i=1 at the receiver. In light of [13], is an ML estimator for σ 2 H that is both statistically sufficient and complete. We thus have an equivalent problem for estimation of the 2nd moment of the Rayleigh random variable based on its N i.i.d. samples. The complete results on sample complexity in this case will be presented in the next section.
After solving the sample complexity problem in the noise-free case, we will investigate the sample complexity problem in estimation of the PL, the noise power, and the SNR, associated with Rayleigh fading channels based on noisy measurement samples. Two sets of independent are transmitted with difference in their average power: |S 1,i | 2 = P s1 and |S 2,i | 2 = P s2 for all possible i where P s1 > P s2 ≥ 0. The following two quantities turn out to be sufficient and complete statistics in statistical estimation of the Rayleigh fading channel based on noisy samples. In fact the optimum estimator is obtained that achieves the Cramér-Rao lower bound. Complete results are obtained and will be presented in a later section.

Sample Complexity in the Noiseless Case
As discussed earlier estimation of the 2nd moment for the Rayleigh fading channel in the noisefree case is equivalent to that for the Rayleigh random variable. Let X be a Rayleigh random variable with PDF (probability density function) Then X 2 = X 2 R + X 2 I where X R and X I are independent Gaussian random variables with zero mean and equal variance σ 2 = µ/2. As it is seen, the PDF of the Rayleigh random variable is completely specified by its 2nd moment µ = E[X 2 ] = 2σ 2 . To estimate µ, it is conventional to obtain N i.i.d. observations X 1 , · · · , X N and take as the estimator of µ. It is shown in [13] that the above estimate is a sufficient statistic with good performance and asymptotically unbiased with large sample size. In addition an algorithm is derived in [13] to compute the interval estimate of µ that is a function of the sample size and involves the use of the tabulated χ 2 CDF. However thus far there does not exist a formula in the literature to estimate the sample size given the interval size and confidence level which will be provided in this section. Different from the result in [13], we begin our investigation on relative error. Let S be an event, and denote Pr{S} as the probability associated with event S. Specifically we aim at solving the following problem: Let ε and δ be positive numbers less than 1; Given margin of relative error ε and confidence level 1 − δ, how large should the sample size N be to ensure Pr Our first result is the explicit formula for computing the sample size as presented next.
Theorem 1 Let ε and δ be positive numbers less than 1. Then the probability inequality (9) holds, if Remark 1 Theorem 1 shows that for small values of ε and δ, the a priori bound on the sample size in (10) can be approximated as As a result, the sample size is roughly inversely proportional to the square of the margin of relative error and is linear with respect to the logarithm of the inverse of the gap between the confidence level and 1. In addition, an important feature of Theorem 1 is that the sample size bound is independent of the value of the true second moment.
Theorem 1 does not use the margin of absolute error as the measure of accuracy. Its reason lies in the fact that if the margin of absolute error is used then the sample size is essentially dependent on the unknown second moment µ, which is to be estimated. Quite contrary, when the margin of relative error is used, the sample size can be determined without any knowledge of the unknown µ. Moreover the relative error is a better indicator of the quality of the estimator: an estimate with a small absolute error may not be acceptable if the true value is also small. Nevertheless Theorem 1 will be used to provide an interval estimate for µ given the sample size and the confidence level or provide sample complexity bound given the estimation interval and confidence level that is similar to a result in [13] but without complicated computation. The sample complexity result for the interval estimate will be presented at the end of this section after the proof for Theorem 1. First we will state and prove the following result that indicates that the a priori bound in (10) is tight as ε and δ become small.
Theorem 2 Under the same hypothesis of Theorem 1, define N (ε, δ) as an integer such that . Then for each µ > 0 and with Z δ satisfying 1 Proof. By the central limit theorem, we can write that verifies the first half of (12). Alternatively since as δ → 0. Taking the natural logarithm on both sides with appropriate rearrangement yields and therefore (12) follows that concludes the proof. 2 Proof of Theorem 1: It can be seen that Note that 2X 2 i µ has mean 1 and possesses a χ 2 distribution with 2 degrees of freedom. Because where the Chernoff bound from [5] is applied. It is easy to see that the infimum is achieved at λ = ln( N θ ) < 0 for which we have e −θ e θe λ −λN = e −θ θe . Using the following inequalities (to be proven later) Finally to prove the first inequality of (16), define It follows that the first inequality of (16) holds. Similarly to prove the second inequality of (16), define Hence the second inequality of (16) is true. The proof of Theorem 1 is now complete. 2

Remark 2
The a priori sample complexity bound in Theorem 1, while being tight asymptotically, can be further improved by taking N = N m1 to be the minimum positive integer such that in light of (13) where F χ 2 ,2N (·) is the CDF of χ 2 random variable with 2N degree of freedom.
Indeed there holds N > Although only the margin of relative error is investigated in Theorem 1 for the sample complexity problem in statistical estimation of the 2nd moment of Rayleigh random variables, a similar result can be obtained as well for the corresponding interval estimate. The probability associated with such an interval estimate will be termed as coverage probability to signify its utility in coverage analysis [4,8].
Corollary 1 (i) Let ε and δ be positive numbers less than 1. Then 19) provided that N satisfies the sample complexity bound in (10). Recall that µ N as in (8) is the maximum likelihood estimate for the 2nd moment.
Then the inequality (19) for the coverage probability holds.
Proof. We note the following chain of equalities: with respect to ε yields the unique positive root ε as in (20). Note that decreases as ε > 0 increases, and ε = 1 is the unique positive root of N = 2(3+2ε) ln( 2 δ ) 3ε 2 at N = 10 3 ln( 2 δ ). Therefore, the root ε is less than 1 for N > 10 3 ln( 2 δ ). It follows that (ii) holds true as well in light of Theorem 1. 2

Remark 3
The problem of interval estimate has been investigated in [13] that proposes a computational procedure but without an a priori bound on the sample size as in Corollary 1. Indeed it computes the exact value of Pr {C 1 µ N ≤ µ ≤ C 2 µ N } = p given 0 < p < 1 by searching for (C 1 , C 2 ) which are dependent on N and ε and satisfy 0 < C 1 < 1 < C 2 . Although in a different form, C 1 and C 2 approach to 1 in accordance with Theorem 2 as N → ∞. However the algorithm as proposed in [13] makes the use of the tabulated χ 2 CDF repeatedly in searching for an appropriate interval and sample size that admit multiple solutions for a given p because of the trade-off between the interval length and sample size. Such a search is bypassed in Corollary 1. In fact such a search is unnecessary that is completely characterized by (17). The explicit relation between N and (ε, δ) enables us to compute an interval estimate for µ a priori given the sample size and the confidence level, or compute the a priori sample complexity bound given the estimation interval and confidence level that is in sharp contrast to the existing results in the literature. Numerical result indicates that our explicit formula in Corollary 1 on the sample complexity bound is very tight as shown on the left of Figure 1 with δ = 0.01 where C 1 and C 2 are solved from (17) and On the right of Figure 1 plots 1+ε 1−ε in unit dB vs. the a priori bound N that is a function of ε, given δ.

Measurement and Estimation with Noisy Samples
Recall that in the noisy case two sets of i.i.d. BPSK symbols are transmitted with difference in their average power satisfying P s1 > P s2 ≥ 0. The problem is statistical estimation for σ 2 H , σ 2 V , and the SNR, as well as their associated sample complexity bound based on i.i.d. observations of the received signal X = HS + V which are recorded as , respectively. We consider the case where N 1 = N 2 = N with X 2 1,N and X 2 2,N as defined in (6). Two different subjects will be investigated.

Interval Estimate
We will demonstrate not only how to construct a simultaneous confidence interval based on the noisy measurements but also how tight these interval estimates are. Our main results begin with the following theorem that provides a simultaneous interval estimate and the solution to the sample complexity problem with an a priori bound.
Theorem 3 Let ε and δ be positive numbers less than 1. Denote ∆P s = P s1 − P s2 . Then Proof. For BPSK the transmitted symbol S is real with the transmitting power P s = |S| 2 , and thus where A and B are independent Gaussian variables with zero mean and unity variance. Hence observations of X, and with X 2 1,N and X 2 2,N as defined in (6), Applying (i) of Corollary 1 yields the interval estimate for an appropriately chosen sample size N . In fact the above holds, if in light of Corollary 1. The fact that Hence the inequality in (22) holds provided that N > 2(3+2ε) ln( 4 δ ) 3ε 2 as in the statement of the theorem. Now by the independence of the experiment, Define the following sets of events as illustrated in Figure 2: : the joint events in the probability expression (21) hold (23) We will show that S 1 ∩ S 2 ⊆ S 0 leading to Pr{S 0 } ≥ Pr{S 1 ∩ S 2 } > 1 − δ thereby concluding the proof. Figure 2: Illustration of the sets of events Indeed suppose that (σ 2 H , σ 2 V ) ∈ S 1 ∩ S 2 represented by the shaded area in Figure 2. Then The above inequalities can be combined to eliminate σ 2 V yielding On the other hand for (σ 2 H , σ 2 V ) ∈ S 1 ∩ S 2 , there hold The above inequalities can be combined to eliminate σ 2 H yielding (σ 2 We can thus conclude that S 1 ∩ S 2 ⊆ S 0 = S H ∩ S V represented by the rectangular area in Figure  2. The proof is now complete. 2 The proof of Theorem 3 shows that two signal sets of different strengths are necessary in order to resolve two different interval estimates for σ 2 H and σ 2 V , respectively. Clearly signals of more than two different strengths become necessary in order to acquire more than two independent parameters in statistical estimation of wireless channels. Theorem 3 also shows that for N > In light of Remark 2 and the definition of F χ 2 ,2N (·), the sample complexity bound in Theorem 3 can be improved by taking N = N m2 , the minimum positive integer such that Recall that SNR = 10 log 10 E[|S| 2 ] + 10 log 10

, it is both important and of independent
interest to obtain an interval estimate for that is presented in the next result.
Corollary 2 Let ε, δ be positive numbers less than 1, and P s2 = 0. Suppose that N is the minimum positive integer such that (28) holds. Then, Proof. In reference to the sets of the events as illustrated in Figure 2, the hypothesis P s2 = 0 implies Then some algebraic manipulations yield that As a result S 1 ∩ S 2 ⊆ S H,V and thus Pr{S H,V } ≥ Pr{S 1 ∩ S 2 } ≥ 1 − δ, if N satisfies the (28) by Remark 3. Therefore the corollary is true. 2 In the case P s2 = 0, the result in Corollary 2 needs to be replaced by a more complex expression: if the positive integer N satisfies (28), and if 1+ε 1−ε P s1 X 2 2,N − P s2 X 2 1,N > 0 and 1−ε 1+ε P s1 X 2 2,N − P s2 X 2 1,N > 0. The above reduces to that of Corollary 2 in the case P s2 = 0.
We investigate next whether or not the interval estimates presented in Theorem 3 and Corollary 2 are tight. Such an issue is clearly dependent on the average transmitting powers P s1 and P s2 . Recall that P s1 > P s2 ≥ 0. Intuitively the larger P s1 and the smaller P s2 are, the tighter the interval estimates are that is validated by the following result.
Theorem 4 Let ε, δ be positive numbers less than 1, and N = N m2 − 1 ≥ N m1 with N m1 and N m2 the minimum positive integers satisfying (17) and (28), respectively. Then for sufficiently large P s1 with constant P s2 . In addition for sufficiently small P s2 with constant P s1 . Finally 1− δ 2 is also an upper bound for the confidence level in Theorem 3 and Corollary 2 with sufficiently large P s1 and sufficiently small P s2 , if N = N m2 − 1.
Proof. In light of Remark 4 and the hypothesis on N , there holds H,N is the limit of σ 2 H,N , a point estimate for σ 2 H to be investigated later: The last inequality follows from the hypothesis on N m2 > N m1 and from the fact that P s1 → ∞ corresponds to the noise-free case in Section 3. Recall S H as defined in (26) that can be alternatively described as where δP s = 1 − P s2 P s1 . It is easy to see that as P s1 → ∞, It follows that the lower bound in (31) is true by taking P s1 sufficiently large. On the other hand It follows that there exists P s1 > 0 such that is large enough that concludes the upper bound in (31). The proof for (32) is similar. Indeed by applying the result from the previous section and using the hypothesis N m1 < N m2 , It follows that both lower and upper bounds in (32) are true by taking P s2 > 0 sufficiently close to 0. Finally the proofs for the rest of the theorem are omitted because they are similar to those for the lower bounds in (31) and (32). 2 Remark 5 Theorem 4 shows that the conservativeness of the coverage probability can be controlled within δ 2 . Since the experimental effort grows in the order of ln( 2 δ ) for a fixed relative width of the confidence interval, a gain of δ 2 for the coverage probability can be obtained at the cost of increasing the experimental effort by It can be seen that this percentage is very insignificant for small δ. Therefore, the conservativeness of the proposed interval estimation can be made very insignificant by choosing large P s1 and small P s2 .
The hypothesis N m2 > N m1 in Theorem 4 is generically true for small δ and ε in light of the fact that N m2 ln 2 δ ≈ N m1 ln 4 δ . Because the limiting case P s1 → ∞ is equivalent to the noiseless case, there hold

Point Estimate
In light of the results in [13], X 2 ℓ,N (ℓ = 1 or ℓ = 2) is the ML estimator for P sℓ σ 2 H + σ 2 V that is both sufficient and complete statistic. Let σ 2 H,N and σ 2 V,N be point estimates for σ 2 H and σ 2 V , respectively based on X 2 1,N and X 2 2,N . Then giving rise to the following expressions for the point estimates (recall ∆P s = P s1 − P s2 > 0): Our next result shows that both σ 2 H,N and σ 2 V,N are optimum estimators for σ 2 H and σ 2 V , respectively.
H,N and σ 2 V,N as in (36) are unbiased ML and MV estimates for σ 2 H and σ 2 V , respectively that achieve the Cramér-Rao lower bound. Furthermore X 2 1,N and X 2 2,N as in (6) are sufficient and complete statistics.

Proof. It can be easily verified that
Thus the estimator in (36) is unbiased. By independence the joint PDF for {X 1,k } N k=1 and {X 2,k } N k=1 is given by by an abuse of notation where for either ℓ = 1 or ℓ = 2, Taking natural logarithm on both sides yields Setting partial derivatives of ln(f X ) with respect to σ 2 H and σ 2 V to zeros lead to the following two equations: Solving for X 2 1,N and X 2 2,N from the above two equations coincide with the expressions in (35) at σ 2 H = σ 2 H,N and σ 2 V = σ 2 V,N that in turn gives the point estimates in (36). The uniqueness of the solution implies that the points estimates in (36) are indeed ML that maximize the PDF f X (·) for every pair of sets of N i.i.d. noisy measurement samples satisfying the hypothesis.
as the parameter vector for estimation where superscript T denotes transpose. Then the corresponding FIM (Fisher information matrix) can be shown to be, after lengthy calculation, where for simplicity in notations, σ ℓ = P sℓ σ 2 H + σ 2 V is used with ℓ = 1, 2. Recall that the inverse of Fim(σ 2 H , σ 2 V ) is the smallest error covariance achievable by all unbiased estimators. In arriving to the above expression, the fact that E[G 4 ] = 3σ 4 for G distributed as N (0, σ), and the observations: are crucial. Now with the point estimates in (36), the corresponding error covariance can be shown to be where again the notation σ 2 ℓ = P sℓ σ 2 H + σ 2 V with ℓ = 1, 2 is used. Consequently there holds the identity Fim(σ 2 H , σ 2 V )Cov(σ 2 H , σ 2 V ) = I, validating the fact that the Cramér-Rao lower bound is indeed achieved by the point estimates in (36) that are indeed MV estimates. Finally the sufficient and complete statistic of X 2 1,N and X 2 2,N can be shown in the same way as in [13]. Specifically the joint PDF in (37) has the form Hence in light of the Neyman factorization criterion, X 2 1,N and X 2 2,N are sufficient statistics. The proof for the complete statistics is again the same as that in [13] that is skipped. 2 The result on the Cramér-Rao lower bound for finite N in Theorem 5 is surprising which is normally achievable only asymptotically. Next we present several results for the point estimates in (36) regarding the sample complexity bound. We begin with a corollary that follows from the results on interval estimates in Subsection 4.1.
Corollary 3 Under the same hypotheses as in Theorem 4, there holds for sufficiently large P s1 with constant P s2 . In addition for sufficiently small P s2 with constant P s1 .
Proof. We prove only (38) as the proof for (39) is similar that will be omitted. Define the set of events: Then S H ⊆ S H , and thus Pr{ S H } ≤ Pr{S H }. More importantly S H → S H as P s1 → ∞. Therefore as P s1 → ∞. A similar argument as in the proof of Theorem 4 can thus be used to conclude (38).

2
Similar to the discussion after Remark 5, we have that for N > where ε and δ are positive numbers less than 1. The next result is concerned with the point estimate related to SNR: Theorem 6 Let ε and δ be positive numbers less than 1. Let P s2 = 0. Suppose that N is the minimum positive integer such that Proof. Define two normalized random variables by Straightforward algebraic manipulations give the following chain of equations: To complete the proof, define two sets of events as follows: It can be shown that S θ ⊆ S θ . Indeed in the plane of ( θ H,N , θ V,N ), the set S θ defines a square area inside the sector area defined by the set S θ . More specifically It follows that as P s1 → ∞, there holds by independence of the events of θ H,N and θ H,N , and by the hypothesis on N . 2 , the next result follows from Theorem 6. Corollary 4 Let ε and δ be positive numbers less than 1. Suppose that P s2 = 0. Then We note that for small ε. Thus the a priori bounds for the sample complexity in Theorem 6 and Corollary 4 are roughly 4 times to that in Corollary 3. It signifies the difficulty in estimation of SNR in Rayleigh channels, and indicates that the margin of error manifested by ε is more expensive to reduce as compared with the level of confidence manifested by δ.

Numerical Simulations
This section presents numerical results to illustrate our proposed statistical theory in measurements and estimation of Rayleigh fading channels. For simplicity σ 2 H = 1 and P s1 = 1 are assumed throughout the section that can always be made true by a suitable normalization. Hence the SNR is the same as , and large P s1 reduces to large SNR. Moreover we consider only the case P s2 = 0.
Because the noiseless case has been shown in Figure 1, we begin with the noisy measurement samples under SNR = 20 dB by generating {X 1,i } N i=1 and {X 2,i } N i=1 for different sample size N . The simple averages X 2 ℓ,N as in (6) are then calculated for ℓ = 1, 2. We choose δ = 0.01 for the associated confidence level. In light of Theorem 3, there holds the joint probability where A 1 , A 2 , B 1 , B 2 are functions of ε, N , and X 2 ℓ,N for ℓ = 1, 2 which represent the interval estimates. In Figure 3 we plotted the ratios for the case P s2 = 0 with SNR = 20 dB (left) and SNR = 60 dB (right). It is interesting to observe that although X 2 ℓ,N , ℓ = 1, 2, are random for each N , the above ratios are almost deterministic owing to the small δ value used. In addition the SNR values affect little for the two ratios at large N values. We also plotted the noiseless ratio in dashed line as a comparison that should serve as a lower limit.
In Figure 4, the above is plotted against the sample size N for the cases SNR = 10 and SNR = 60.  It is commented that for both Figure 3 and Figure 4, the a priori bound N > 2(3+2ε) ln( 4 δ ) 3ε 2 , rather than the relation in (28), is used. The results are nevertheless close to each other. In Figure 5, we plotted ε vs. N (left) for the case δ = 0.05, and plotted δ vs. N (right) for the ε = 0.05. In both cases, a priori bounds are used in both the noisy and noiseless cases, plus the use of χ 2 in the noiseless case that is governed by the relation in (17). As expected, the curve based on χ 2 CDF in the noiseless case serves as a lower bound. In addition it is observed that the two curves based on the a priori bounds are almost identical for the one on left, and close to each other for the one on right. This fact indicates that our results on measurement and estimation of noisy Rayleigh fading channels are not conservative. It is further observed that the ε curves decrease at a constant slope in log-log scale, and the slope of decrease is rather slow with respect to the sample size N . This fact indicates that the reduction of ε (with fixed δ) is expensive in terms of increasing N . On the other hand the δ curves decrease at accelerated slopes with respect to the sample size N , implying that the reduction of δ (with fixed ε) is relatively cheap, especially at large sample size N . Figure 5 validates that the measurement sample size is roughly inversely proportional to the square of the margin of error ε and is linear with respect to the logarithm of the inverse of the gap δ.
The last simulation example is worked out to demonstrate the optimum estimator obtained in Theorem 5. Different from the previous two cases, M (>> 1) sets of measurement samples are taken to assess the average performance for the underlying statistical estimation. For each sample size N , {X ℓ,i (k)} N i=1 is generated for ℓ = 1, 2, and 1 ≤ k ≤ M with M = 500. The optimum estimator in (36) is used to compute estimates σ 2 H,N (k) and σ 2 V,N (k). The estimation error is then averaged to yield that are plotted together with the Cramér-Rao lower bounds: CRB H = (P s1 σ 2 H + σ 2 V ) 2 + (P s2 σ 2 H + σ 2 V ) 2 N ∆P s , CRB V = P 2 s2 (P s1 σ 2 H + σ 2 V ) 2 + P 2 s1 (P s2 σ 2 H + σ 2 V ) 2 N ∆P s versus the sample size N as shown in Figure 6 where P s2 = 0 is used. It can be seen that the RMSEs coincide very well with the CRBs, validating the optimality of the estimator in (36). Again the change of SNR affects little on estimation of σ 2 H but changes the RMSE for σ 2 V that is mainly due to the change of σ 2 V by a factor of 10. In the case M < 500, the RMSE lines are less straight and fluctuate more as M becomes smaller, but the overall trend holds.

Conclusion
Statistical estimation of Rayleigh fading channels has been investigated based on both noiseless and noisy measurement samples. Complete solutions are derived for the associated sample complexity problem and provided for the optimum estimator problem in measurement and estimation of Rayleigh fading channels. Specifically our a priori bounds on measurement sample sizes ensure the prescribed margin of error and confidence level and are contrast to the existing work reported in the literature. In dealing with the noisy measurement samples, our proposed novel signaling scheme with two different signal strengths is instrumental in extracting the statistical information on mean channel power, noise power, and SNR. Such a novel signaling scheme enables us to derive the sample complexity bounds for both interval and point estimates that are tight for the mean channel power and noise power albeit less tight for the SNR. More importantly it leads to the optimum estimator that is both an ML and MV estimator and that achieves the Cramér-Rao lower bound. The results presented in this paper constitute an independent statistical theory for measurement and estimation of Rayleigh fading channels. It should be emphasized that the sample complexity solution in the noiseless case is also instrumental without which the results for the case of noisy measurements are not possible. The numerical simulations illustrate that our proposed statistical theory is effective in statistical estimation of Rayleigh fading channels. Specifically the simulation examples indicate that our results based on the noisy measurement samples are close to that based on the noiseless measurement samples, and the signal power is not required to be high that can be compensated for by using large sample size. Currently we are investigating the Nakagami fading channel which in a special case reduces to the Rayleigh fading channel, and aim at extension of our results on the measurement sample size and optimum estimator. It is our objective to apply our results to other more general wireless fading channels and to broad applications of our proposed statistical theory.


Generalization of usual capability indices for unilateral tolerances

Capability indices are dimensionless quantities measuring the aptitude of a process to manufacture items whose characteristics must be within a specified tolerance ranges. The usual indices $C_p$, $C_{pk}$, $C_{pm}$, and $C_{pmk}$ are used for a process of normal distribution and a target located at the center of the tolerance interval. Various indices derived from the previous family allow to consider more complex situations when asymmetrical tolerances and non-normal distributions are taken into account. In this paper we study the case where a single tolerance is imposed because the shifts in the direction of this tolerance appear much more serious than in the opposite direction. We propose a family of four indices having interpretations and properties similar to those of the usual family.

Introduction
The quality of an industrial process is assured by the monitoring of one or more variables of interest. The process will be considered of good quality if it is able to produce items whose variable of interest is within the tolerances L and U (Lower and Upper tolerance limits) and close to a specified target value T. Accordingly, the capability indices used to measure the quality of the process are linked to the location and the dispersion of the supervised variable. For a target centered in its tolerance interval and a variable of normal distribution, the first work of Kane [9] was completed to lead to a family of four indices generally used. X being the supervised characteristic of normal distribution N(μ,σ), one defines When asymmetrical tolerances and normal data are concerned, many proposals have been made, the most coherent of which as the generalization of the usual family, is given by the family '' p C , " pk C , '' pm C , and " pmk C introduced by Pearn, Chen, and Lin in various articles [3,4,12,16,17].
In the case of non-normal distributions, Clements [5] proposes to generalize the indices p C and pk C replacing the mean μ by the median M and the dispersion 6σ by Up -Lp where Up and Lp are the 99.865 and 0.135 percentiles. In the same way, Pearn and Kotz [15] provide a extension of the indices pm C and pmk C . However let us note that considering the difficulty to obtain reliable estimations of Up and Lp, many other proposals have been provided the bibliography of which can be found in the articles of Tang and Than [21], Mac Cormack, Harris, Hurwitz, and Spagon [11], or Ding [6].
The literature related to the processes limited by two tolerance limits is important, but it is not the same in the less frequent situation where only one limit is imposed. In this case two situations can be considered. The first one occurs when the variable of interest, because of its nature, cannot exceed a certain level which represents the target of the process. It is the case for example for concentricity or circularity where the observed measure is obviously positive or null, the target being equal to 0. The second one occurs when a drift of the mean in a direction appears much less serious to the user, so that he is induced to define only one single tolerance. It is this second situation which is the subject of this paper. We recall the rare indices which one finds in the literature, then propose a family of four indices having interpretations and properties, similar to the ones of the usual family p C , pk C , pm C , and pmk C . For normal data, the densities and the moments of the natural estimators are given. The expressions used in the case of a normal distribution can be easily generalized to the case of non-normal distributions by replacing the mean μ by the median M, and the natural variation The estimation of the extreme percentiles Up and Lp is however far from reliable on samples of reasonable size when one uses the traditional moment estimation, as Clements does [5]. The work of Shore [19,20] enables us to obtain estimations much more reliable which will be developed in an example.

Existing indices
Starting from the definition of ( ) in order to measure the capability of a process in the unilateral tolerance situation. Let us note that the CPU and CPL indices do not take into account the existence of a target which, for the index pk C , is supposedly, implicitedly located at the center of the tolerance interval. When the target is not centered, Kane [8] suggests an index referred to as * pk C from which he gets the indices ( ) , in the case of a unilateral tolerance. In the same way as with the usual indices, the value of these indices is equal to 0 if the previous calculus leads to a negative value. Chan, Cheng and Spiring [2], in the case where one-sided tolerance is required, have suggested generalizing pm C to ( ) . Moreover in order to generalize all the previous suggestions, Vänmann [23] puts forward two families of indices, , and . As Vänmann points out, the fact that there is an unilateral tolerance can be interpreted in such a way that a shift of μ away from T towards that tolerance is more serious than a shift of μ towards the

Suggestions of indices for a normal distribution
Our objective is to build four indices generalizing the usual properties and interpretations of the indices p C , pk C , pm C , and pmk C . Let us mention that for a normal distribution and a target centered in the tolerance interval, p C is linked to the proportion of non conforming items when the mean is on the target, and is equal to 1 for a proportion of 0,27% of non conforming. Moreover it is meant as the potential capability of the process, that is to say the maximum capability which one can obtain for a given dispersion, when the mean is on the target. The modern standard of quality deems that a process should not be considered capable if μ is far away from T, even if σ is small. p C which is unrelated to T does not satisfy this requirement. For this reason the indices pk C , pm C , and pmk C , which take into account the location of the process mean as well as the process variability, have been introduced. The 3 of them are maximum and equal to p C when the mean is on the target. The deviation of the mean is taken linearly into account by pk C , so that it is null when the mean reaches the tolerance limits. Thus, the ratio pk p C C allows to determine the position of the mean between the target and the tolerance towards which it deviates. pm C takes the deviation into account in a quadratic way. Thus, it is not null at the tolerance limits, but takes the same value. pmk C being the combination of pk C and pm C allows to take the deviation into account in a "quadratic" way and to obtain a null index when the mean reaches one of the tolerances. Thus, in the case of unilateral tolerance, in order to keep interpretations similar to the bilateral case, we will require a) that the potential capability indices pu C and pl C take value 1 for a 0,135% proportion of non conforming, when the mean is on the target b) that the indices of capability taking the position of the mean related to the target into account, are maximum and equal to potential capability when the mean is on the target. c) that the indices pku C and pkl C decrease linearly and take value 0 when the mean reaches the tolerance limit d) that the indices pmu C and pml C decrease in a quadratic way e) that the indices pmku C and pmkl C decrease in a "quadratic" way and take value 0 when the mean reaches the tolerance limit. The main difficulty in the building of indices lies in the fact that we have no knowledge of the risk taken when the mean moves away from the target in the opposite direction to the tolerance limit. However, even if there is no tolerance limit, obviously, a production manager cannot accept a too large deviation in a direction even if, a priori, it does not seem too serious to him. Thus, it appears fundamental to us that he should quantify the "not too serious". Is this twice, five times, ten times less serious? We consider thereafter that it has been decided that the risk is k time less serious. The choice of the constant k (>1) being rather approximate we require a last condition: f) when the means deviates towards the tolerance limit, the capability indices must be independent of the choice of k. Let us consider for the moment the case of an upper tolerance U, and note By using notations similar to Vänmann's [23], we can write the general formula This expression gives the four indices for the couples ( , ) u v = (0,0), (1,0), (0,1) and (1,1) again. Note that the letter u used in subscript is an abbreviation of the word upper and is independent of the first parameter of the indices family. The indices ( , ) It is obvious that when T μ = and pu , and a proportion of 0.135% of items Thus for h = ½ by example, the mean is halfway between the target and the tolerance limit. For T μ > ,  to the left is considered to be three times less important than to the right.
In a similar way, if the single tolerance is a lower limit L, we obtain the general formulation

Suggestions of indices for a non-normal distribution
As in the bilateral case, we replace the mean μ by the median M and the natural variation 3σ When ( , ) u v = (0,0), (1,0), (0,1) and (1,1), the previous formulas give generalizations of the indices p C , pk C , pm C , and pmk C , in the case of non-normal distributions. Assuming that , we find the same algebraical expressions and thus the same properties as those stated in the case of a normal distribution. Practically, to be able to use these results we need to estimate the percentiles M, p U , and p L . The determination of a percentile is easy when the distribution from which the observations ensue is known. To identify this distribution, the most usual method consists in fitting the distribution to a member of a family covering a great number of usual distributions. To achieve this adjustment the method of the moments requires the identification of 3 or 4 parameters according to the family being used. If the estimate of the mean and of the standard deviation is rather reliable, it is no longer the same for skewness and kurtosis, which are manifestly known for their great dispersion. Thence, the resulting estimates are not at all reliable for the extreme percentiles p U and p L . Thus for a non negative variable, Shore [19] suggests approaching the p percentile p Q by the relation By identification of the complete moments M Y the ith partial moment of Y. Shore's method leads to better estimates of the percentiles than those obtained by Clements's method [4], insofar as the expected values being similar, the mean squared error is much lower by Shore's method [19,20,21].

Distribution and moments of the estimators of the indices suggested for a normal distribution.
The studied characteristic of the process is supposed normally distributed with mean μ and variance σ² . Two natural estimators of In the case of a single upper tolerance U, considering that the risk of a deviation towards the left of the target is k time less serious than towards the upper limit U can be interpreted as the positioning of a lower limit L so that l u T L D kD − = = . In these conditions ( ) Assuming that

Estimation and distribution of pu C and pl C
As previously, if we consider that the choice of the constant k can be interpreted as the positioning of a lower limit L so that ( 1)/2 ( 1)/2 Moreover, since In a similar way we obtain the density and the moments of , pl n C and , 1 pl n C − replacing pu C by pl C .

Estimation and distribution of pku C and pkl C
If we consider that the choice of the constant k can be interpreted as the positioning of a lower limit L so that x given by Pearn, Lin and Chen [18].
Only the moments of order 1 and 2 of '' , 1 pk n C − are explicitly given by Pearn and Chen [13].
Thus, we use the equation (3) for the r-th moment which, in addition, leads to expressions of the first two moments simpler than those given by Pearn and Chen [13].

Estimation and distribution of pmu C and pml C
If we consider that the choice of the constant k can be interpreted as the positioning of a lower limit L so that (1).
From the appendix, The moments are obtained from the equation (3):  x The moments are obtained from the equation (3): 700 685 490 610 685 590 620 700 730 750 780 730 770 760 780 640 830 785  780 680 650 585 650 600 610 640 640 575 630 610 665 690 650 640 645 630 700  800 685 650 675 640 655 610 625 580 610 Figure 3 gives the histogram of the 86 observations studied. A khi square test carried out on the classes of this histogram gives a p-value of 0.003. This p-value leads us to reject the normality of the observations. Thus we use the indices defined for a non-normal distribution and estimate the quantiles by the method of Shore [18]. We obtain M = 685, 1  L T . The histogram shows however that no unit is below the tolerance, but that the median estimated at 685 ppm is largely above the target T = 480 ppm. However, even if the company assesses that this deviation is less serious than in the other direction, the resulting additional cost requires the use of indices of real capability taking into account the position of the median according to the target. As an example we will consider the 2 situations where the gravity was judged 3 times (k = 3), then 10 times (k = 10) less important. For k=3, we obtain ˆp when the deviation takes place on the right, the smaller the risk is, the closer to the potential capability the indices of real capability are.

Conclusion
In this paper we first review the existing generalizations of the basic capability indices p C , pk C , pm C , and pmk C which have been proposed to handle one-sided tolerances processes. In this situation, the risk of a deviation towards the side opposed to the limit of tolerance is not considered very serious and thus, in general, not quantified. However, for problems of cost among other things, this risk cannot be completely ignored. When a ratio is established, at least approximately, between the risks taken for deviations to the left and to the right of the target, we introduce indices of capability preserving the properties and interpretations of the usual family p C , pk C , pm C , and pmk C . For a process of normal distribution, the densities and moments of the estimators are given. They ensue from the properties established for the family " ( , ) p C u v , used in the case of asymmetrical tolerances. When the distribution of the process is not normal, we propose simple extensions using the median and the quantiles of order 0.00135 and 0.99865. Finally an example shows the way to calculate these quantities on a real example of non-normal data, using the method of Shore [18] in order to estimate the extreme quantiles.


extfieldsofstudy: Physics Mathematics

text:
Manifold and metric in numerical solution of the quasi-static electromagnetic boundary value problems

Classical vector analysis is the predominant formalism used by engineers of computational electromagnetism, despite the fact that manifold as a theoretical concept has existed for a century. This paper discusses the benefits of manifolds over the traditional approach in practical problems of modelling. With a structural approach, it outlines the role and interdependence of coordinate systems, metric, constitutive equations, and fields, and relates them to practical problems of quasi-static computational electromagnetics: mesh generation, open-boundary problems, and electromagnetic-mechanical coupled problems involving motion and deformation. The proposed procedures also imply improvements to the flexibility of the modelling software.

Introduction
The objective of solver software systems for electromagnetic boundary value problems (BVP) is to solve an arbitrary problem from a specific problem class. However, this is not possible in practice: One of the restricting factors is the difficulty of the mesh generation. This is related to the fact that to pose a BVP, the user needs to define a geometry of the BVP domain. The geometry is described with a coordinate system and certain choices of coordinate systems may cause too much round-off error with floating point arithmetic.
Often the problem geometries lend themselves to a particular coordinate system, and therefore the software systems usually offer a list of coordinate systems to choose from. To offer a user interface with a list of coordinate systems, many structures in the software are fixed: For example, the software often has one intrinsic coordinate system, usually Cartesian, which it uses in all calculations. Furthermore, the basic material parameters, such as the permittivity of a vacuum, as well as the metric (inner product) in the intrinsic coordinate system are pre-decided. When the user gives the model's geometry, for example, in the spherical coordinates, the software maps them to the intrinsic Cartesian system.
Frequently, the floating point problems leave the modeler without any mesh at all and prevent the solution of the BVP: Practical computations with the intrinsic coordinate system requires that the coordinate numbers are represented with floating point numbers. This makes mesh generation particularly challenging for BVPs with details of the domain varying significantly in size. For example, consider the power line on top of figure 1. The thickness of the cables and the structures in the supports are much smaller in scale than the distance between the two supports. Meshes are difficult to generate because the floating point representation gives us labels to only a finite, although large, number of points. For example, as a floating point number grows in magnitude, the distance between points separable by numeric labels increases, and, at some point, round-off errors become big enough to make labels of distinct points coincide unintentionally. Generally, mesh generation is based on geometric predicates [1], some of which are difficult to calculate with floating point numbers. Mesh generation becomes hard or even impossible with double precision floating point numbers when the scale variation, the ratio of the largest and the smallest dimensions, is in the order of 10 4 or higher. This estimate for scale variation is based on the experiences with the BVP presented in the section 7.
The mesh generation problems arise from the way the BVP was described in a particular coordinate system. Is it possible to alleviate these problems? With the help of modern mathematics, manifolds and differential geometry, the role of mathematical structures on the description of a BVP can be outlined. In this paper we strongly promote a structural approach. Particularly, we will show that for a given BVP, the elements of the triplet {coordinate system, metric, material parameters} are interdependent, and the triplet constitutes an entity. In fact, there is a whole class of these triplets for each BVP: Each triplet describes a geometry of the BVP domain and the constitutive equations, which are needed to define an electromagnetic BVP. Specifically, the numerical values of material parameters are respect to given metric. For example, the numerical value of the permittivity of vacuum is different with respect to meters than with respect to inches. Although the metric appears in the triplets, the electromagnetic theory essentially concerns the works related to displacements of particles and the total energy of systems. The role of metric is central only in the classical view-vector analysis-whereas the modern view of differential geometry reveals the independence of the electromagnetic theory from a particular metric. This has been known for long [2][3] [4][5] [6], but not that well recognized in engineering. Hence, the modeller may select the metric at will. This is much greater a freedom in the context of manifold than a redefinition of the inner product in the single codomain vector space of a vector field in the classical view. We use the modern view to show how to formulate electromagnetic problems with the classical vector formalism in different metrics. From the software design point of view, there is no need to fix the whole triplet in a software system. Just like the software system is not intended to solve a particular problem, but any problem from a class of problems, the solver software should provide the possibility to use any triplet from the triplet class. It should not be even restricted to a single triplet, because manifold promotes local coordinate systems that provide further flexibility. This will help overcome floating point problems in the calculations of the geometric predicates and thus will enlarge the class of practically solvable problems: Our approach suggests giving up the metric and the coordinate system strongly suggested by human intuition. For example, if the mesh generation fails with the coordinate system in the top of figure 1, one could use the coordinate system in the bottom of the same figure, when the calculations of the predicates are easier. This time it is more likely that one gets a mesh and thus a numerical solution. However, there is usually a price to pay: numerical error can be large because the mesh quality may be poor.
The problems of mesh generation are just one example of practical challenges in numerics that can be engaged with tools provided by manifold. The problems with motion and deformations, such as the deformation effects of magnetostriction, can be modelled with a single mesh: motions and deformations can be modelled as changes in the metric. This can speed up the solution process substantially. Finally, the solutions for open boundary problems, that are solved by finite domains with so-called infinite elements in the boundary region, are an example of an application where a convenient triplet is chosen for solution.

Structural Background
We need many mathematical structures to pose and describe an electromagnetic BVP. The starting point of this description is a point set that is the BVP domain. The points of this set correspond to observations on distance measurements. Then we will add one structure at a time to this point set, such that finally we have all that is needed for the description.

Metric Space
At first we need to specify a point set for the BVP domain. The points in the domain correspond to the points observed in reality, in the following sense: The observed points are considered distinct if and only if a nonzero distance can be measured between them. We denote this abstract point set by M . By the way the set M was constructed, it is naturally endowed with the structure of metric space [7]. The distance between a given pair of points refers to an observation of lengths in multiples of some reference object called the rigid body. (Observe here that numerical values of distances depend on the chosen rigid body but distances can be discussed without a specified reference body.) Furthermore, the metric induces the structure of topological space for M [7]. The dimension of M , which is two or three, must be the same at every point of space M .

Manifolds
For computational analysis, we need to characterize displacements arithmetically and thus must parameterize the point set. In the classical view this is done with coordinate systems, which assign unique coordinates-that is, a tuple of real numbers-to each point of M . Formally this means a mapping of type M → R n . Furthermore, R n is endowed with the norm topology and the mapping is homeomorphic, i.e. continuous bijection with continuous inverse. However, the abstract point set M is often ignored completely and the points of M are identified with the coordinate numbers. Furthermore, often a metric connotation is given to the coordinate numbers. Thus the domain of the BVP is taken to be some particular subset of R n . Yet the coordinates are just labels for the points of M , such that they allow arithmetic and need not carry metric information. The coordinates can be changed with diffeomorphic mappings of type R n → R n .
Instead of a particular coordinate system, the modern view uses manifolds as domains for BVPs. In manifolds, the abstract point set M is the primary object, which is not identified with any particular set of coordinate numbers. The coordinate systems, or charts, are considered secondary, which is natural because we can choose the charts in many ways, and the real numbers assigned to each point are not canonical. In particular, the physics-work and energy-does not depend at all on the chosen chart. Furthermore, M is not expected to be coverable by a single chart, as was the case for coordinate systems: Charts are homeomorphic mappings from the open connected subsets of M to the open subsets of R n . Moreover, every point of M is required to be in the domain of some chart. Any topological space, that can be covered with charts, is called a topological manifold [8]. Thus the structure of topological manifold emphasizes the existence of charts, rather than a particular chart.
Topological manifold is not enough for physics, because the differentiability of functions needs an extra structure. Let us require a collection of charts, that are bound together by their transition maps, or change-of-chart maps: The differentiability of a real-valued function h defined on M is defined indirectly through the differentiability of the functions h • f −1 , where f is a chart. The chart can be changed with the transition map f • g −1 as follows: We do not want the n-times differentiability of h to depend on the choice of chart, and therefore require that the transition map be (at least) n-times differentiable. If this holds, the charts are considered equivalent. This relation constitutes an equivalence class, a differentiable structure, of charts. The set M together with a differentiable structure is a differentiable manifold or here just manifold for short [8][9][10].

Standard Parameterization
The rigid body can be used as a pair of dividers to specify spheres in M . If the image of every sphere in M under a chart is a sphere also in R n when measured with the Euclidean 2-norm, we call the chart a standard parameterization. For example, the chart f in figure 2 is a standard parameterization. All standard parameterizations are identical up to scaling, rotation, and translation (plus reflection). The standard parameterizations are important because in literature the numerical values of the material parameters are given in a specific class of standard parameterizations. This is natural because the numerical values of the material parameters are meaningful only with respect to a given metric, and in practice, the given metric is always based on a length unit system in which some rigid body is chosen for reference.

Tangent Spaces
Charts enable us to talk about displacements in small neighborhoods, i.e. The electric field is an object that yields the virtual electromotive force dU (virtual work up to charge) for every virtual translation v. In other words, the electric field is a functional field over the tangent spaces. Correspondingly, the electric flux is an object that yields the virtual flux for every pair of virtual translations, which define a virtual surface. That is, the electric flux is also a functional field over the tangent spaces. The solutions of quasi-static BVPs are such functional fields.

Inner Product
In the classical view we use a vector field to express the electric field. This possibility is due to Riesz representation theorem [12] [13], that enables us to represent the electric field as a vector field, once an inner product is available. In the classical view the inner product is defined in the codomain vector space, but in the modern view each tangent space of M has its own inner product. The pointwise representation of the electric field functional by a vector E requires an inner product into each of them. At each point, with the inner product denoted by (·,·), the virtual emf dU (v) corresponding to the virtual translation vector v can be represented as dU (v) = (E, v). This representation makes the inner product an inherent part of the classical view : the vector E depends on the chosen inner product whereas dU does not. When the tangent spaces of a manifold are equipped with an inner product which varies smoothly from point to point, the manifold is called a Riemannian manifold [8][9] [10]. The inner products on the tangent spaces of Figure 4: Correspondences of tangent vectors under different charts. A tangent vector v at p ∈ M corresponding to vectors dr f at f (p) and dr g at g(p).
M are usually called metric tensors. They induce a metric 1 on M and also bring about the notion of angles between the tangent vectors [8].
For practical computations with tangent vectors, such as evaluations of metric tensors, we need the counterparts of the tangent vectors and the metric tensors on R n . For each tangent vector v at p, a unique counterpart vector (the so-called push forward [11]) dr f exists in R n , induced by the chart f , see figure 4. We can use these counterpart vectors dr f to define metric tensors on R n such that the tensors uniquely correspond to the ones on the tangent spaces of M : For each metric tensor (·, ·) M given for M , there is a unique equivalent metric tensor (·, ·) f for R n , that satisfies at all points and for all pairs dr f , dr f of R n and all pairs of corresponding tangent vectors v, v . Equivalent metric tensors on R n and M make f an isometry. However, the metric tensor chosen for R n need not be equivalent with the one on M , i.e. the charts need not be isometries.

Equivalent Descriptions of BVP
The theory of electromagnetics is independent of the chosen metric tensor and the chart, as stated in the introduction. In fact, Maxwell's equations do not need a metric tensor at all. They are meaningful on a differentiable manifold (not necessarily a Riemannian manifold) [14][15] [16]: the virtual works related to the displacements of a point charge are independent of the metric tensors and chosen charts. Naturally, the energy stored in a system does not depend on the chosen charts either. However, the constitutive equations need to be described to define the energy and there are at least two ways to proceed: First, a metric tensor is selected in order to represent the fields as vector fields, and constitutive equations are defined couple the two vector fields. Second way to proceed is to couple the functional fields directly,in which case the constitutive equations couple different types of fields. Although a metric tensor is not needed to construct the constitutive equations in this case, constitutive equations induce metric tensors on the manifold [16]. Yet in neither strategy does the energy itself depend on the selected or the induced metric tensor. Thus an electromagnetic BVP is independent of particular chart or metric tensor, and it is uniquely specified with manifold M , Maxwell's equations, the constitutive equations, and the boundary values. However, the material parameters, or rather the numeric values of permittivity, permeability, and conductivity depend on the chosen metric. Furthermore, to describe an abstract manifold concretely, we need to do that by charts, even though the manifold itself is independent of particular charts. Therefore, to pose the BVP, we need charts and metric tensors. When a chart and a metric tensor field are chosen, the numerical values of the material parameters can be specified. Then the material parameters are realor matrix-valued fields defined on the chart. Thus we have a triplet {chart, metric tensor, material parameters} that describes a geometry of the BVP domain and the constitutive equation. Because we can choose any chart from the differentiable structure, and define any metric tensor on the chart, there are infinitely many geometries that we can use to pose the BVP. In fact, there is an equivalence class of these triplets {chart, metric tensor, material parameters}, that describe the same unique BVP. We call the equivalence relation producing these classes the material equivalence.
We can now express a BVP on any chart: we need a triplet {chart, metric tensor, material parameters}, Maxwell's equations, and the boundary conditions. To express the BVP on any other chart, we only need to select another triplet, because Maxwell's equations are independent of the metric and of particular chart, as are the boundary conditions. Notice, that we can also keep the same chart and change the metric tensor and the material parameters. These descriptions of the BVP have equivalent solutions; that is, in any description, the electric vector field corresponds to the same virtual emf functional field. We obtain a new view of BVPs from the material equivalence: a BVP can be seen as a material equivalence class together with Maxwell's equations and boundary conditions. This view suggests an approach to construct solver software systems: Solver for BVPs should allow the user to select any triplet from the material equivalence class.
This way the user can use the triplet that is the most suitable for numerical solution.
The benefit of choosing any triplet is that we can solve numerically a larger class of problems. Let us derive the material parameters for given chart and metric tensor, such that the material equivalence holds with another triplet that is known.

Material Parameters for Given Chart And Metric Tensor
To determine the material parameters for given chart and metric tensor, we require that the energy stored in a system and the virtual works related to the displacements of a point charge be invariant of our choice of chart and metric tensor. For simplicity, let us now focus on electrostatics, because other quasi-static cases can be treated conceptually in a similar fashion. Following the standard convention, we use dr · dr to denote (dr, dr ) in R n . Also, for each metric tensor field · in R n , there is a unique symmetric positive definite matrix field S such that dr · dr = dr T S dr holds pointwise for all dr and dr . Let us assume we know a triplet {f i , · i , i }, consisting of chart f i of M , metric tensor · i on R n and the material parameters i . Next, we choose chart f j from the differentiable structure and use metric tensor · j on R n and then determine the material parameters j for the triplet {f j , · j , j } such that the two triplets are equivalent.

Invariance of Virtual Work
The counterparts or the push forwards of a tangent vector v at p ∈ M are dr i and dr j under charts f i and f j , respectively. The virtual emf dU (v) corresponding to the virtual translation v is E i · i dr i , or E j · j dr j , where · i and · j are metric tensors in R n . These expressions must equal pointwise: Let : ran(f i ) ⊂ R n → R n be the transition map from the range of f i to the codomain of f j . Denoting by J the Jacobian matrix [17] Substituting (3) for dr j in (2), and relying on the invariance of the virtual work, we get where S i and S j are the matrix presentations of the metric tensors · i and · j . Because (4) holds for every point and for every virtual translation dr i , we can express E i in terms of E j and matrices S i , S j , and J as follows: If we know either E i or E j , we know dU for every virtual translation, hence the solution of the BVP.

Invariance of Energy
The energy W stored in the electric field in f i (M ) can be expressed as where the matrix i contains the material parameters in the chart f i . Equation (6) induces an inner product for the fields. Here E k and E l are any two vector fields. Even though there is a metric tensor in the expression of the inner product of the fields, the energy is independent of chosen metric tensor: the product S i i describes the constitutive equation, such that the value of energy is the same regardless of the choice of S i . This is what is meant by the claim that electromagnetic theory is independent of the chosen metric. Note that the inner product defined in (7) allows us to formulate the Galerkin method. For any electric field, the energy stored in the field should be invariant under the choice of chart, or must hold for any equivalent fields E i and E j . Using the change of variables theorem [17], we can write the left hand side of (8) in the codomain of f j as where   J −1    is the determinant of the inverse of the Jacobian matrix J. Substituting the expression given in equation (5) for E i in (9), we get Now, because this holds for any vector field E j defined on f j (M ), the matrix j can be written in terms of matrices i , S i , S j , and J, as follows:

Equivalence of Numerical Solutions
The solutions for equivalent triplets are equivalent from a mathematical point of view, and so are the numerical solutions: Next we define equivalent meshes and show that under equivalent meshes, the system matrices in FEM are identical up to the round-off errors of the floating point representation. First, assume a mesh in M . The mesh points in M have corresponding points in every codomain of the chart of every triplet; consequently, we get an equivalent mesh for every chart (an example of equivalent meshes in figure 5) 2 3 . Assuming real number arithmetics, equation (8) implies that the system matrices related to them are identical by construction. Thus the FEM solutions are also equivalent. However, in practice the system matrices are assembled and the numerical solutions are calculated with floating point numbers, and this leads to round-off errors. Thus the system matrices are identical only up to the round-off errors.
Although the solutions are equivalent, in practice every BVP has tripletssometimes including the standard parameterizations-that do not allow mesh generation, and thus no numerical solution. In these cases, mesh generation may well be possible with some other triplets. However, if the mesh quality criteria are not changed according to the chosen triplet, then there is usually Figure 5: Example of equivalent meshes. Two equivalent meshes for the codomains of the charts f and g of figure 2. The mesh in the codomain of the standard parameterization f has long elements in the middle, which may imply a poor quality. However, it is not always possible to generate any reasonable mesh for the standard parameterizations. a price to pay: the mesh quality may be poor, which causes numerical error in the solution. For example some elements may be very long in the standard parameterization as seen in top of figure 5.

Applications
We now have all the necessary tools to express a BVP, attempt a solution, and interpret it with any triplet, and next suggest how to take advantage of this possibility in practice. First, we show how a chart can be changed when the metric tensor in R n is fixed, which is usually the case. Second, we suggest a scheme to use multiple charts in mesh generation and problem solution. Third application is open boundary problems which are solved by selecting a triplet in which the problem domain in the chart is bounded. Finally, we discuss how to model motion with only one mesh by changing the metric.

Reparameterization Under Fixed Chart Metric
Because the numerical values of the material parameters are always given in literature in some length unit system, it is easiest to describe the constitutive relations with isometric standard parameterizations. Consequently, the programmer of a finite element software system often assumes that the user applies only isometric standard parameterizations; that is, the user can control only the unit of length, and rotate and translate the model. Furthermore, due to the isometry assumption, the metric tensor is implicitly fixed such that it induces the Euclidean 2-norm. Therefore the matrix representing the metric tensor in the intrinsic coordinate system of the solver is implicitly hardwired to be the identity matrix.
Let us assume that we know how to express a BVP in some isometric standard parameterization f , that is, we know the triplet {f , S f , f }. However, instead of using f directly, we would like to reparameterize with chart g, which is not a standard parameterization, to avoid floating point problems in numerical modeling. Because of the same hardwired metric tensor in the codomains of f and g, the distances between any two coordinates (tuples of real numbers) do not change when changing chart from f to g. However, a change of chart makes the coordinates correspond with different points in the manifold M , and thus the change of chart in this case implies a change of metric, hence the chart g is not isometry.
To set up our BVP correctly under g, we must impose the material equivalence on g. For this we can use only the material parameters g , because g is chosen by the modeller, and the metric tensor is hardwired by the programmer. Equation (11) shows that the material parameters g require information on g relative to f . Yet everything related to f , including f itself, exists only in the modeller's mind. We are solving the BVP using a triplet {g, S g , g }, and we interpret this description of the BVP under g as being equivalent to that under f . However, the software programmer who assumes only standard parameterizations and has no information about f , might think of g as a standard parameterization and very likely interpret the problem geometry very differently. In this sense, we would be "misusing" the software.
To impose the material equivalence, we need the Jacobian between the charts and material parameters f . After substituting I for S f and S g in (11), we get the following equivalent material parameters: where the identity (Henrotte et al. [18] arrived at a similar expression but one of their motivations was open boundary problems, see section 6.3). It is easy to see from (12) that the material parameters in g may not be a scalar multiple of the identity matrix, even if they were in the chart f , a point expressed in [16] from the opposite point of view.
The solution of the problem are fields E g and D g . However, to present the fields in the standard parameterization f , where they may be more intuitive, we must use (5) and D f = f E f to obtain the equivalent fields.
Remark 1 D g cannot be transformed into D f with the same transformation as E g is transformed into E f in (5). That is, the electric field intensity E and the electric flux density D transform differently between charts. To get the D f , one can use (5) to transform E g to E f and then use the constitutive The above method is fully realizable with any mesh generator and solver software that allows for pointwise definition of material parameters as a matrix. However, if the solver does not allow for pointwise definition, only charts g with piecewise affine transition maps from the standard parameterization f can be used.

Modelling With Multiple Charts
Previously we had two charts that both cover the whole domain. However, manifold can be covered with multiple charts, such that only a part of the manifold is covered by each chart. All that is required is that every point of the manifold belongs to the domain of some chart and that all the charts belong to the same differentiable structure. The use of multiple charts can benefit numerical analysis, for example in the case of figure 6: First of all, the mesh generation gets easier, because the detailed regions could be covered separately with their own charts, with the origins located such that there are as many floating point numbers available as possible. Second, for the same reason, the calculations of the inner products for the system matrix are more accurate. Even though multiple charts are used, we still assume here that the manifold can be covered with one R n chart. This makes it possible to partition the manifold into separate domains such that the domains only overlap on their boundary. Furthermore, this guarantees that the manifold can be embedded into the n-dimensional Euclidean space. Next, based on these assumptions, we suggest an implementation of multiple charts into software systems.
First, the user implicitly decides on one isometric standard parameterization, which covers the whole domain. Let us call it the universal chart, because it works as the reference to the other charts, and the material parameters are given on it. However, the universal chart is not used in the calculation or in the mesh generation, but these are done in separate charts that cover only parts of the manifold: the user partitions the universal chart into multiple regions, such that the regions only overlap on their boundary (see figure 6). Then the user gives new charts for these regions and gives the transition maps to the universal chart. In practice the user only gives the separate charts and the transition maps to the software system, and then the universal chart is constructed from these by the software. It is also possible that the universal chart is not constructed at all in the software system but exists only implicitly. The metric tensors for these separate charts could be fixed or given by the user. The material parameters for the charts are given as if they are standard parameterizations, and the software calculates automatically the material parameters to be used in the calculations. This is possible because the metric tensors are known, as are the transition maps, whose Jacobians can be automatically calculated.
Finally, in order to attain compatibility of meshes in different charts, the meshes have to agree at the boundaries of the regions. The mesh generation could proceed as follows: First a mesh is generated on one of the regions with g f Figure 7: Open boundary problem. Top: The problem domain (manifold) consist of an device over the ground. Bottom: The codomains of two charts, f and g. The regions inside of the rectangles denoted with dashed line correspond to the same part of the domain, but in g the rest of the space is between the rectangles. the chart that covers it. Then the coordinates of the nodes on the common boundary with a second region are mapped with the transition map to the chart that covers the second region, and the rest of the mesh in the second region can be generated. The assembly of the system matrix can be done as usually.

Open Boundary Problems
The open boundary problems are often solved by truncating the domain and setting the fields zero at an artificial distant boundary. This usually gives a good enough approximation for the fields near the sources, because the fields tend to zero quite fast as the distance from the source increases. The distance of the artificial boundary from the sources and regions of importance is decided by the modeller. However, the error from this truncation is hard to estimate and the number of elements needed to cover all that empty space can be large.
To overcome these problems, many methods have been proposed over the years: The so called "ballooning method" [19], in which the true distance of the boundary is put very far away with a thin layer of special elements.
Another method uses infinite elements [20], which are special decaying basis functions used in FEM. Yet another method couples the FEM with analytical solutions as in [21]. Finally, the transformation (or shell transformation) methods presented for example in [18] [22] put the infinite boundary to a finite distance with the help of a suitable change of coordinates. However, many of these methods are hard to implement into most production solvers, because they require modification of the solver codes: For example, the modeller needs to define new basis functions or needs to give pointwise Jacobians of the change of coordinates.
The triplet approach of BVPs suggest a strategy to solve open boundary problems: Just select a suitable triplet, such that in the chart the infinity boundary is at a finite distance, such as the chart g in the figure 7. This is effectively identical to many of the above methods. For example, in the transformations method the Jacobian of the change of coordinates is implemented in the software code whereas in the triplet approach the Jacobian is included in the material parameters given by the user.
The triplet should be selected such that a small region including the source regions and other regions of interest is modelled as usually, but the complement of this region is mapped to a small bounded outer region as in the chart g in the figure 7. The elements nearest the infinity boundary of the domain are such that the points at the boundary are mapped to the infinity by the transition map form g to the corresponding standard parameterization. A suitable transition mapping could be the following: The region outside the unit disc is mapped to the ring with exterior radius two, such that the ring surrounds the unit disc and the boundary of the disc is the inner boundary of the ring, see figure 8. A point outside the unit disc at the distance r from the center of the unit disc is mapped to a point at the distance R from the center, such that the distances have the following relation: This method can be used with almost any solver: Only the material parameters have to be changed accordingly, which is a responsibility of the modeller. However, the procedure is easy to automate: the user only specifies the regions where the outer boundary is mapped to the infinity by the transition map from g to the standard parameterization, and the software system would take care of the rest. The system could have built-in such a mapping or mappings which adapt according to the coordinates of the infinity boundary and the boundary of the infinite regions. If necessary, the user gives the correspondence between these boundaries and then the system modifies the user-given material parameters according to the adapted mapping.

Modelling Motion With Single Chart
When modelling motion, we have to solve a sequence of BVPs whose only difference is that some object changes its position relative to the other objects. For an example, consider the force calculations of an electromagnet, see figure 9: The magnet attracts the load, which then moves upwards. One solution is to make whole model for each position of the moving object. For each position, we have to define a new geometry and generate a new mesh for it, and furthermore assembly the whole system matrix again. This is very time-consuming.
The triplet approach can speed up the process: Keep the chart and the material parameters fixed and change the metric tensor. This changes the geometry and takes into account the motion: When a rigid object moves in a fluid (air) such that the topology does not change, i.e. the moving object does not touch any other material than the fluid, then the underlying differentiable manifold of the domain does not change (this underlying manifold is often referred to as material manifold). However, the Riemannian manifold will change, because the metric tensor changes. This means that the distances between points of the manifold are changed. Thus we may select a single chart g from the differentiable structure, to cover all the different Riemannian manifolds.
The metric tensor for the selected chart g again depends on the material equivalence. This requires that some standard parameterization f for each Riemannian manifold is known, when the corresponding Jacobians are also I I Figure 9: Modelling motion. Electromagnet: coil with current I, ferromagnetic core and load, which is attracted upwards.
known. Then the matrix S g representing the metric tensor on g can be solved from equation (11): Assume the material parameter f to be a scalar and require the material parameters to be equal, that is g = f . Then S g in terms of the Jacobian can be solved: Now, for the example in the case of figure 9, the metric tensor only changes in the air between the magnet and the load according to the equation (13). Because the mesh remains the same, only the elements in the system matrix corresponding to the mesh elements in this air region will change. Thus only a partial re-assembly of system matrix is needed for each new step. Furthermore, in the case of iterative solvers, a preconditioner used for one step can be used effectively for many successive steps, because the changes in the system matrix can be very small. Also the previous solution may be a good initial guess. Thus, all these things-only one mesh, partial re-assembly of system matrix, the preconditioner, and the initial guess-can speed up the solution process substantially.
Many solvers do not offer the possibility to change the metric tensor, yet the modelling of motion with a single mesh is still possible: we can equivalently change the material parameters, see section 6.1. One needs to know a standard parameterization for each Riemannian manifold and then use the material parameters given by the equation (12). Finally, all things said about modelling motion apply also for modelling the geometric effects of deformations, such as the effects of magnetostriction.

Numerical Example
We demonstrate a combination of proposed applications: an open boundary problem (section 6.3) with reparameterization under fixed chart metric (section 6.1). The problem to solve is the electric potential of a three-phase high-voltage line, a 3D-Laplace's problem. Figure 10 shows a standard parameterization of the domain. The domain is the half-space above the ground apart from the line itself. The infinity boundary is mapped to a finite distance from the line, as explained in section 6.3. There is also great variation in the scale of the details of the domain: the length of the lines is order of hundreds of meters, whereas the smallest dimensions are order of centimeters. Based on experience, we can say that the mesh generation with the standard parameterization is difficult or even impossible. To avoid mesh generation problems, the lines are scaled down, as is the height of the pillars. The model used in the calculations is shown in figure 11. The result is shown in figure 12. The figures and the calculations are produced with GetDP [23] and Gmsh [24].

Conclusion
Most solver software systems for quasi-static electromagnetic boundary value problems are based on the formulation of electromagnetics with classical vector analysis. This has many restrictions for numerical solution, because many mathematical structures have limited freedom or they are completely fixed. For example, the fields are presented as mappings from some global coordinate system to a single vector space and the metrics for the coordinate systems and the vector space are fixed. The restrictions limit the possible coordinate systems that can be used. The formulation of electromagnetics with manifolds and differential geometry is less restrictive and thus can help in a numerical solution. This view exposes the independence of the electromagnetic theory on the chosen coordinates or metric. Based on this it is shown that there is an equivalence class of triplets {coordinate system, metric, material parameters} for each boundary value problem. A triplet describes a geometry of the problem domain and the constitutive equations which are needed to pose the problem in the coordinate system. We propose that the solvers should allow to use any triplet, that is, the modeller should be able to choose any coordinate system and metric at will. These possibilities come built-in in differential geometry, whereas they are not self-evident in classical vector analysis. We discuss how this can be done. It is not hard to apply some of the triplet strategies in most production solvers, because not all of them require modifications to the code. We show that despite the freedom in the choice of the triplet, the FEM system matrices of a given BVP, assembled for a mesh represented under different triplets, are identical.
We show how the triplet approach can be exploited in applications: Potentially prohibitive difficulties of mesh generation in problem domains that have large scale variation can be alleviated by a suitable choice of coordinate system. This extends the class of solvable problems, although possibly at the expense of mesh quality. Open boundary problems can be solved by choosing a convenient triplet for the solution. Solution of problems involving motion and deformations can be accelerated when a single mesh is used on every time step and only the metric is altered.