Gps-tsc.upc.es

NLOS MITIGATION BASED ON A TRELLIS SEARCH FOR WIRELESS LOCATION
Andreu Urruela, Helena Morata and Jaume Riba Signal Processing and Communications Group. Technical University of Catalonia (UPC) Campus Nord, Ed. D5, Jordi Girona 1-3, D5, 08034, Barcelona (Spain) Emails: {andreu,alhelena,jriba}@gps.tsc.upc.edu ABSTRACT
estimating the variance of the incoming observation, so this ap-proach is not very suitable for high dynamic scenario.
Wireless location using Time of Arrival (TOA) and Time-Difference- Other papers, like [4] are based on the idea that the NLOS Of-Arrival (TDOA) measurements has received considerable at- process could be intermittent. This means that the BSs can be in tention over the last years to be the best selection for cell phone LOS or in NLOS in an intermittent way. This model is intuitively location. The major problem related with these types of mea- correct when we are talking about in-car applications where the surements is the Non-Line-Of-Sight problem (NLOS) that happens mobile is supposed to be moving and seeing the BSs in LOS or when the direct path between the Base Stations (BSs) and the mo- NLOS intermittently. The idea proposed in [4] consists in exploit- bile is blocked. This paper presents a new technique to mitigate ing this intermittent effect of the NLOS to discard the NLOS mea- the NLOS effect in dynamic scenarios based on a Trellis search surements. The problem here is that in [4] it is assumed that the of the NLOS state. Numerical simulation shows that the proposed mobile is static during the observation period, which is contradic- technique outperforms the previous contributions because it is able tory with the hypothesis that the NLOS is intermittent.
to detect and reject the NLOS measurements.
Also based on the idea that the NLOS effect can be intermit- tent, other papers like [5] and [6] try to pre-filter the timing mea- 1. INTRODUCTION
surements to, somehow, eliminate the effect of the NLOS from themeasurements. The techniques used are a biased Kalman filter in Wireless location has received considerable attention in the last [5] and a filter based on a certain scattering model in [6]. Obvi- decades. This is partially motivated by the fact that location al- ously, the authors in these last papers made assumptions about the gorithms will be mandatory for public cell networks, and partially nature of the NLOS errors, so the performance of such algorithms motivated by the interesting market in position-based services as is logically based on how well the real measurements errors fit this car-routing or position-based publicity. Another interesting feature model. On the other side, this type of approaches does not use the is the possible optimization of the network planning based on the relationship between the measurements captured simultaneously from different BSs. This is, all the timing measurements coming To estimate the position of the mobile, several techniques have from a certain BS are pre-filtered independently from the rest of been proposed during the last years. The most successful tech- measurements, not exploiting the relationship between them.
nique seems to be based on timing measurements such as TOA In this direction, [7] tries to exploit the relationships between or TDOA. The advantages of using these observations are that the three TOA measurements captured simultaneously from different accuracy achieved is in the order of tens of meters, that is what BSs, trying to mitigate the NLOS errors. However, it is difficult to is required for emergency calls, and that the investments needed extent the proposed algorithm to more TOA measurements or even other types of measurements. Maybe the most significant contri- The main obstacle between the use of TOA or TDOA obser- butions in this direction are [8] and [9]. In both articles, the authors vations towards an accurate location system is the Non-Line-Of- formulate all the possible hypothesis of LOS/NLOS for a specific Sight effect. This effect mentioned as the killing effect in [1], is by instant of time, i.e. each hypothesis is characterized by which BSs far the most important cause of big errors in the positioning algo- are under LOS and which BSs are under NLOS. Latter on, the ML rithms based on timing measurements. Note that NLOS errors can position estimates are obtained for each hypothesis. Finally the be in the order of hundreds of meters.
solution presented in [8] consists in a linear weighted combination Good algorithms to mitigate the NLOS effect have been pre- of the partial position estimates associated to each hypothesis. The sented in the previous literature. One of the first approaches pre- algorithm presented in [9] tries to identify what is the most likely sented was based on the observation that the timing measurements hypothesis using the incoming timing measurements based on the captured in a NLOS channel are normally corrupted with a stronger application of the ML detection principle. The major problem with noise, i.e. bigger variance. Using this concept, [2] and [3] pro- these last approaches is that they exploit only the relationship be- posed to discard the timing measurements with bigger variances tween timing measurements captured at the same time (one snap- using a certain detection criteria. The main problem here is that shot), but they do not exploit the relationship between the timing the variance of a timing observation in a LOS scenario depends measurements along the time axis as in [5] and [6].
strongly on the type of scenario, i.e. urban, suburban or rural. An- This paper tries to extend the work presented in [9] to dynamic other problem is that the detection process could take some time scenarios, exploiting also the relationship between the timing ob- servations among different snapshots of measurements. In this This work has been partially supported by the European Comission under IST project EMILY IST-2000-26040, by the European Commis- sense the new algorithm exploits the relationship between mea- sion (FEDER) and by the Spanish/Catalan Government under projects surements in both directions: inside the snapshot and along the TIC2003-05482, TEC2004-04526 and 2001SGR-00268 time axis. The proposed algorithm assumes, like previous ap- proaches [4], that the NLOS effect is intermittent, but it is coherent A common assumption in wireless location is that the multi- also with assuming that the mobile is in a dynamic scenario. No path terms wk,n are independent in time among sources, and Gaus- specific assumption of the NLOS error is done, so the proposed al- gorithm exploits only the coherence between the incoming timingmeasurements using certain assumptions about the mobile trajec- E [wk,nwk ,n ] = σ2kδk−k δn−n tory. As it will be seen, the natural extension of [9] for dynamic scenarios leads to an unafordable complexity, so a Trellis-based p is the Kronecker Delta. For a generalization with un- known variance, we refer the reader to [10].
search algorithm is proposed to track the state (LOS or NLOS) of previous one but three common assumptions are done in previous A model for the intermittent effect of the NLOS is also pre- sented to simulate the proposed algorithm. Simulation results willshow that the proposed algorithm outperforms the previous ap- (a) The NLOS term for each source (subindex k) is indepen-
proaches and it is near to the optimum that is ideally defined as dent from that of the other sources.
the ML estimator that only uses the timing measurements under (b) In NLOS conditions, a significant bias will be added to the
measurement and this bias will be positive for TOA mea-surements. The distribution for this bias is usually assumed 2. SIGNAL MODEL
as uniform but no standard models have been presented sothey will be assumed to be completely unknown.
Let us assume an scenario where L BSs are under the visibility (c) The NLOS effect in dynamic scenarios is intermittent, so
of the mobile. In this scenario, the mobile can periodically obtain the NLOS term will be zero when the measurement is un- snapshots of L TOA measurements associated to each one of the der LOS conditions and will be more or less constant in BSs or L − 1 independent TDOA measurements associated to all the independent pairs between the reference BSs and the neigh-bor BSs. To unify the notation of the development for the caseof TOA and TDOA observations, let us define K as the number 3. SINGLE SNAPSHOT ESTIMATION
of measurements we have at each instant of time (each snapshot),i.e. K=L for TOA measurements and K=L − 1 for TDOA mea- This section presents the ML estimation of the position at the n-th surements. These K measurements will be referred to as sources instant of time, xn, using only a single snapshot of K measure-
(because they produce a new measurement at each instant of time) ments as shown in (1). The goal here is to obtain the ML estima- indicating single BSs for the TOA case and pairs of BSs for TDOA tion of xn being robust against the NLOS effect. This represents
an introduction to the multiple snapshot case and a review of pre- Let us define the K observations (k=[1,K]) taken at the n-th If we have K measurements that can be in LOS or in NLOS, there are Q = 2K possible hypothesis regarding the NLOS state.
tk,n = gk (xn) + wk,n + lk,n
For the q-th hypothesis, we can define γLOS vectors containing the subindex of the sources that are supposed to where wk,n is the noise term caused by the additive noise and the multipath effect in the TOA/TDOA estimation, lk,n is the error Under the q-th hypothesis, we can divide the K original sam- caused by the NLOS effect and gk (xn) is the non-linear function
ples tk,n in two subvectors tLOS
and tNLOS
relating the observations and the parameter of interest, i.e. xn that
measurements that are in LOS and NLOS respectively as follows: is the position of the mobile at the n-th instant of time.
The non-linear function gk (xn) in (1) does not depend on the
tk,n ∈ tLOS
temporal subindex n and it is defined specifically for the type of tk,n ∈ tNLOS
In a coherent way we can define also the non-linear functions as- (xn)
xn − xBS
(xn)
xn − xBS
xn − xBS
gk (xn) gLOS
n,q (xn)
where xBS
is the known position of the k-th BS.
gk (xn) gNLOS
(xn)
If the mobile captures N snapshots of measurements at a con- stant rate of r snapshots per second, we would like to model the In the same direction, we can divide the original K variances σ2k evolution of the mobile position at each instant of time as a func- in two diagonal matrices containing the variance associated with tion of a finite number of parameters as: the measurements that are supposed to be in LOS and NLOS re-spectivelly under the q-th hypothesis as follows: xn = hn (x, u)
k ∈ diag RLOS
n (x, u) does not depend on the source subindex k and
shows the mobile position at the n-th instant of time based on the position x at the end of the observed window and the rest of move-
k ∈ diag RNLOS
ment parameters included in u such as speed, acceleration,etc. For
instance, for a constant speed linear trajectory, we have: and finally, we can define the vector containing all the unknown hn (x, u) = x + (n − N)s
where u = s is the constant unknown speed vector of the mobile.
lk,n ∈ ln,q ⇔ k ∈ γNLOS
Using these definitions, we can define the ML estimation of 4. MULTIPLE SNAPSHOT ESTIMATION
the position xn and the unknown NLOS biases ln,q under the q-th
hypothesis (x
The goal now is to extend the work presented in [9] to the multi- n,q and ln,q respectively) as follows:
ple snapshot case using the movement model presented in (4). The xn,q, ln,q = arg min ΨLOS
major drawback in [9] was that the decisions about the LOS/NLOS n,q (xn) + ΨNLOS
(xn, ln,q)
xn,ln,q
hypothesis were taken instantaneously at each instant of time as shown in (22), non exploiting the high correlation between the po- sition in two consecutive instants of time shown in (4) or (5) northe high correlation in time of the NLOS phenomena.
n,q (xn) = φ tLOS
gLOS
n,q (xn) , RLOS
The goal now is not estimating the mobile position at each in- stant of time in an independent way, but also estimating the mobile (xn, ln,q) = φ tNLOS
gNLOS
(xn) ln,q, RNLOS
position at the end of the observed window x and the movement
φ (m, R) = ln |R| + mT R1m
parameters u shown in (4).
Taking into account that we have Q possible simple hypothe- As shown in [11], under the assumtion that xn is known, the
sis regarding the NLOS state for each instant of time, now we have ML estimation of ln,q can be obtained from (15) as:
P = QN global hypothesis regarding the NLOS state for the entireobservation window of N consecutive snapshots of measurements.
ln,q (xn) = tNLOS
gNLOS
(xn)
Let us define qp as the vector containing the simple hypothesis for
and now using (17) in (15), the compressed ML estimation of x
each instant of time under the p-th global hypothesis. This is, its -th element, qp (n) [1, Q], is the simple hypothesis selected
for the n-th instant of time under the p-th (p ∈ [1, P ]) global hy- xn,q = arg min ΨLOS
n,q (xn) + ΨNLOS
xn, ln,q (xn) (18)
Under this p-th global hypothesis we have as unknowns, the xn,q = arg min ΨLOS
position at the end of the observation window x, the movement
n,q (xn)
parameters u and the set of unknown NLOS biases for each instant
that can be solved using a closed form algorithm like in [12] or [13]. Note that the second term in (18) becames constant with respect to x
n and can be removed from the minimization. This p = lT
1,qp(1), . . . , lT
N,qp(N)
last expression (19) basically demonstrates that, if we do not have any knowledge of the NLOS biases, the ML position estimator Analogously to the previous section, the ML formulation of is obtained only using the measurements that are supposed to be all the unknown parameters under the p-th hypothesis can be for- Another conclusion from (19) is that the covariance matrix of n,q, Cn,q = E (xn,q − xn) (xn,q − xn)T
p, up, lp = arg min
p (x, u, lp)
x,u,lp
χp (x, u, lp) =
Cn,q ≈ GTn,q (xn,q) RLOS
Gn,q (xn,q)
n,qp(n) (hn (x, u)) + ΨNLOS
n,qp(n) hn (x, u) , ln,qp(n)
gLOS
gLOS
where matrix Gn,q (xn)=
n,q (xn)
n,q (xn)
[xn]1
[xn]2
n,q (xn) and ΨNLOS
(xn, ln,q) are defined in (16).
Until now, we have Q possible position estimates xn,q , q =
Again, as shown in [11] and assuming that all the unknown 1, · · · , Q corresponding to the ML position estimates under the Q parameters included in lp are independent, we can obtain the ML
possible hypothesis regarding the NLOS state. In order to obtain a estimation lp and compress the ML function χp (x, u, lp). Con-
final position estimate considering the Q possible hypothesis, two cretely, the ML estimation of l
approaches have been proposed to solve this last step. Pi-Chun p, lp = lT
1,qp(1), . . . , lT
N,qp(N) is
Chen proposed in [8] to obtain the estimations as a combination of the Q partial estimates as follows: ln,qp(n) = tNLOS
n,qp(n) gNLOS
n,qp(n)(hn (x, u))
and the compressed ML estimations of x and u are obtained as:
xn=
n,q (xn,q)
n,q (xn,q) xn,q
xp, up = arg min χp x, u, lp
x,u
In [9], the authors proposed just to extend the ML search into xp, up = arg min
n,q
the hypothesis domain applying the ML detection principle as fol- x,u
p(n) (hn (x, u))
Here, we obtain a similar result: the mobile position x and
xn , q = arg min ΨLOS
the movement parameters u are obtained as if we had only the
n,q (xn) + ln Γ1
available LOS measurements (under the p-th hypothesis).
where Γq are some tuning parameters to modify the a-priori proba- Finally, if we apply the same idea shown in (22), we can select bility of each hypothesis and to compensate for the different num- ber of unknown parameters in each hypothesis.
As exposed in [9], this second alternative obtains better per- x, u, p = arg min
n,qp(n) (hn (x, u)) + ln Γ1
formance, specially in severe NLOS scenarios since it is able to x,u,p
5. PROPOSED ALGORITHM
of simple hypothesis associated to all the instants of time until thecurrent iteration. These Q partial hypothesis will be the most likely It is clear that the complexity associated with the last expression hypothesis that end with a particular NLOS state.
makes it useless for a real implementation since it is not possible At the M -th instant of time, these Q partial hypothesis are de- to test P = 2K·N global hypothesis to select the best one. This fined by the vectors βM
q , q = 1, · · · , Q. These vectors contain section presents the proposed simplification for a real implemen- the simple hypothesis for each instant of time until the current in- tation. The proposed idea is inspired in the coding theory based on stant of time indicated by M , this is βM
q (n) [1, Q] n ≤ M , the Trellis Coded Modulation (TCM), where it is not possible to eval- n-th element of βM
q , is the simple hypothesis considered for the uate the likelihood of all possible paths, so only the most likely paths survive at each iteration. While in TCM schemes we try to Coherently with the definition, the last element inside the vec- extract the codified source information, here we try to extract the tor βM
NLOS state at each iteration. While in TCM schemes there are q has to be q, this is βM
q (M ) = q. Note that βM
length vector since it contains only the simple hypothesis selected forbiden transitions, here we have low probability transitions (i.e.
for the first M instants of time. Note also, that each partial hypoth- all NLOS to all LOS). Finally, while in TCM schemes we extract esis corresponds to a family of the P global hypothesis since the information from the received signal to, somehow, estimate the partial hypothesis shows only the simple hypothesis selected until probability of each state based on the previous probability of all the states, here we do exactly the same.
Let us assume then that we are at the M -th instant of time and First of all, note that (28) shows the ML estimation of the posi- that we have Q partial hypothesis defined by the vectors βM
tion x and model parameters u using only the LOS measurements
clear that initially, we can form Q ·Q partial hypothesis as follows: we have at different instants of time defined by the p-th globalhypothesis. This estimation can be approximated in a two-stepsestimation process. In the first step, partial position estimates are βM
q,q = βM
q = 1, · · · , Q q = 1, · · · , Q obtained for each instant of time using only the LOS we have. Inthe second step the partial position estimates are fused using the Generically, βM
q,q contains all the simple hypothesis considered movement model shown in (4). This two steps estimation process until the M -th instant of time by βM
can be mathematically expressed as follows: about the NLOS state for the M + 1-th instant of time.
Since we have Q new vectors that end with a certain hypoth- xp, up = arg min
φ xn,q
x,u
p(n) hn (x, u) , Cn,qp(n)
esis q , these are βM
q,q ∀q, we have to select the most likely one using a simplified version of (32). Concretelly we implement this where xn,q is the partial estimate under the q-th hypothesis for the
n-th instant of time as defined in (19) and Cn,q is its approximate
βM+1
βM
Under this approximation, we can redefine the hypothesis se- qsel = arg min Ξ βM
x, u, p =
Note that the search for each state (q ) is only over the Q pos- sibilities consisting in the transitions between all the previous state φ xn,q
and the new one. It is clear that the unafordable complexity of (32) x,u,p
p(n) hn (x, u) , Cn,qp(n) +ln Γ1
qp(n) (31)
has been dramatically reduced since here we implement a searchof Q candidates instead of the QM candidates that we had in the where it can be seen that all the P hypothesis are based on the partial position estimates and their covariance matrices. This Finally, it can be seen that the minimization over the posi- is, the P = QN global hypothesis are all the ways to combine the tion x and movement parameters u in (33) can be implemented
partial position estimates that we have per each one of the N iteratively with a linear Kalman filter reducing the complexity of evaluating the minimization at each iteration.
For clarity reasons, let us isolate the hypothesis selection in p = arg min Ξ (qp, N)
6. NUMERICAL SIMULATIONS
Ξ (q, M ) =
Numerical simulations have been performed in a scenario with fiveBSs uniformly distributed in a circle of radius 2 Km centered in the origin of coordinates. The mobile is moving with a constant φ xn,q(n) hn (x, u) , Cn,q(n) + ln Γ1
x,u
qp(n)
speed of s =
of coordinates in the middle of the time observed window. The The cost function Ξ (q, M ) shows the likelihood of a certain
measurement rate is r = 1 snapshot of TOAs per second and the hypothesis expressed by the vector q. This includes the simple
length of the observed window is N = 100 snapshots.
hypothesis for the M instants of time.
Joining the assumptions (a), (b) and (c) shown in section 2
Note that (32) is just showing that we decide the hypothesis p about the NLOS effect, we have simulated the NLOS effect of based on the likelihood of the possible P global hypothesis using each BS as an independent process. The NLOS phenomena for all the N available snapshots of measurements.
each BS is modelled with a Markov Chain with two states: LOS The idea we propose is to limit the search over the entire P and NLOS. There is a transition at each instant of time with the possible hypothesis, iteration by iteration as in TCM.
The proposed technique consists in only considering Q partial hypothesis at each iteration. A partial hypothesis is a collection performs the previous contributions since it is able to exploit a movement model over the intermittent effect of the NLOS.
7. REFERENCES
[1] M. I. Silventoinen and T. Rantalainen, “Mobile station emer- gency locating in gsm,” IEEE International Conference on Personal WirelessCommunications, pp. Page(s): 232 –238, [2] M. P. Wylie and J. Holtzman, “Non-line of sight problem in mobile location estimation, the,” IEEE InternationalConfer- ence on Universal Personal Communications, vol. Volume: Optimum per iteration Simple ML−Detection (max 2 NLOS) 2, pp. Page(s): 827 –831, September 1996.
Simple ML−Detection (max 1 NLOS) Pi−Chun Chen approach (max 2 NLOS) [3] J. Borras, P. Hatrack, and N. B. Mandayam, “Decision the- oretic framework for nlos identification,” IEEE Vehicular Technology Conference VTC, vol. Volume: 2, pp. Page(s):1583 –1587, May 1998.
Fig. 1. C.D.F. of the position error with σ
[4] M. P. Wylie-Green and S. S. Wang, “Robust range estimation in the presence of the non-line-of-sight error,” IEEE Vehicu-lar Technology Conference VTC, vol. VTC 2001 Fall. IEEE VTS 54th, Volume: 1, pp. Page(s): 101 –105, 2001.
[5] N. J. Thomas, D. G. M. Cruickshank, and D. I. Laurenson, “A robust location estimator architecture with biased kalman filtering of TOA data for wireless systems,” IEEEInterna- tional Symposium on Spread Spectrum Techniques and Ap- plications, vol. Volume: 1, pp. Page(s): 296 –300, September2002.
[6] S. Al-Jazzar, Caffery J., Jr., and H. R. You, “A scattering model based approach to NLOS mitigation in TOA locationsystems,” IEEE Vehicular Technology Conference VTC, vol.
Spring 2002. IEEE 55th, Volume: 2, pp. Page(s): 861 –865, Optimum per iteration Simple ML−Detection (max 2 NLOS) [7] S. Venkatraman, Caffery J., Jr., and H. R. You, “Location Pi−Chun Chen approach (max 2 NLOS)Pi−Chun Chen approach (max 1 NLOS) using los range estimation in nlos environments,” IEEE Ve- hicular Technology Conference VTC, vol. Volume: 2, Spring- 2002, pp. Page(s): 856 –860, 2002.
[8] Pi-Chun Chen, “A non-line-of-sight error mitigation algo- Fig. 2. C.D.F. of the position error with σk=60 meters
rithm in location estimation,” Wireless Communications andNetworking Conference, vol. vol.1, pp. Pages:316 – 320,September 1999.
Each time there is a transition from LOS to NLOS, a new bias is [9] J. Riba and A. Urruela, “A non-line-of-sight mitigation tech- selected for the samples lk,n using a uniform distribution in the nique based on ml-detection,” IEEE International Confer- range [200,400] meters. This bias is used in all the consecutive ence on Acoustics, Speech, and Signal Processing (ICASSP), iterations until the next transition from NLOS to LOS.
vol. Volume: 2, pp. Pages:ii – 153–6 vol.2, May 2004.
Figures (1) and (2) show the sample cumulative density func- [10] Andreu Urruela and Jaume Riba, “Efficient mobile loca- tion (CDF) of the position error of the proposed technique when tion from time measuremenets with unknown variances in the standard deviation of the TOA measurements is σk = 30meters dynamic scenarios,” Signal Processing advances in wireless and σk = 60meters repectivelly.
In the same figure, two theoretical limits have been plotted for [11] Yihong Qi and H. Kobayashi, “Cramer-rao lower bound for comparison. We can see the Optimum per iteration limit consisting geolocation in non-line-of-sight environment,” IEEE Inter- in selecting always the LOS measurements at each iteration and national Conference on Acoustics, Speech, and Signal Pro- compute the position xn iteration by iteration. On the other hand,
cessing (ICASSP), vol. vol 3, pp. Pages:III–2473 – III–2476, we can also see the Optimum limit consisting in selecting always the LOS measurements and compute the mobile position x and the
[12] A. Urruela and J. Riba, “Novel closed-form ML position es- movement parameters u exploiting the movement model presented
timator for hyperbolic location,” IEEE International Confer- in (5). Logically, both limits can not be implemented since it is not ence on Acoustics, Speech, and Signal Processing (ICASSP), vol. Volume: 2, pp. Pages:ii – 149–52, May 2004.
For comparison, we also present the performance of the al- [13] Y. T. Chan and K. C. Ho, “A simple and efficient estimator gorithm presented in (22) [12] as simple-ML-detection and the for hyperbolic location,” IEEE Transactions on Signal Pro- algorithm presented in (21) [8] as Pi-Chun Chen algorithm. In cessing, vol. Volume: 42 Issue: 8, pp. Page(s): 1905 –1915, both cases, we have limited the NLOS hypothesis search to 1 or 2 NLOS. It can be seen how the proposed algorithm clearly out-

Source: http://gps-tsc.upc.es/comm/jriba/NLOS%20mitigation%20based%20on%20a%20Trellis%20Search%20for%20Wireless%20Location.pdf

Microsoft word - ateliers def gb.doc

8TH LASAIRE BIENNIAL 13 AND 14 JANUARY 2005 WORKSHOPS PROGRAMS Workshop 1 : How Europe can work efficiently with 25 member countries? Workshops leaders: Mario Dehove, Joël Maurice, Antonio Lettieri Thursday 13 January 2005 from 2:30 p.m to 4:15 p.m ROUND TABLE A: How should we assess Europe's economic performance? Speakers: Eva Belabed, European Economic

Pharmacy

Pharmaceuticals Industry Profile From: Yahoo Finance, September 13, 2005 A The US leads the world in pill popping, claiming both the largest market share and five of the ten largest druggernauts (Bristol-Myers Squibb, Johnson & Johnson, Merck & Co., Pfizer, and Abbott Laboratories). Pfizer, already the top drug company, pushed itself even farther out in front with its 2003 buy of

Copyright © 2010 Medicament Inoculation Pdf