1

Sequential Channel Estimation in the Presence of
Random Phase Noise in NB-IoT Systems

arXiv:1706.04350v3 [] 29 Sep 2017

Fredrik Rusek and Sha Hu
Department of Electrical and Information Technology, Lund University, Lund, Sweden
{firstname.lastname}@eit.lth.se

Abstract—We consider channel estimation (CE) in narrowband
Internet-of-Things (NB-IoT) systems. Due to the fluctuations in
phase within receiver and transmitter oscillators, and also the
residual frequency offset (FO) caused by discontinuous receiving
of repetition coded transmit data-blocks, random phase noises
are presented in received signals. Although the coherent-time
of fading channel can be assumed fairly long due to the lowmobility of NB-IoT user-equipments (UEs), such phase noises
have to be considered before combining the the channel estimates
over repetition copies to improve their accuracies. In this paper,
we derive a sequential minimum-mean-square-error (MMSE)
channel estimator in the presence of random phase noise that
refines the CE sequentially with each received repetition copy,
which has a low-complexity and a small data storage. Further, we
show through simulations that, the proposed sequential MMSE
estimator improves the mean-square-error (MSE) of CE by 1
dB in the low signal-to-noise ratio (SNR) regime, compared to a
traditional sequential MMSE estimator that does not thoroughly
consider the impact of random phase noises.

I. I NTRODUCTION
Narrowband Internet-of-Things (NB-IoT) is a new radioaccess technology designed for ultra-low-end IoT applications
[1]–[4]. In the latest Long-Term Evolution (LTE) Rel-13, the
3rd Generation Partnership Project (3GPP) [3] has introduced
a number of key features for NB-IoT. A user-equipment (UE)
operates in the downlink (DL) of NB-IoT system using 12
subcarriers with a subcarrier bandwidth 15 kHz, and in the
uplink (UL) using a single subcarrier with bandwidth either
3.75 or 15 kHz, or alternatively 3, 6 or 12 subcarriers with
a subcarrier bandwidth 15 kHz. A legacy global system for
mobile communications (GSM) [5] operator can have a “standalone” NB-IoT deployment by replacing one GSM carrier
(200 kHz) with NB-IoT, while an LTE operator can allocate
one physical resource blocks (PRBs) of 180 kHz to NB-IoT
inside an LTE carrier, or deploying NB-IoT in the guard-band,
which are know as “inband” and “guard-band” deployments,
respectively. NB-IoT reuses the LTE design extensively to
achieve a harmonious co-existence of current networks.
To offer a wide-area coverage and improve the robustness
of data-transmissions, repetition coding [4] is used for both
DL and UL in NB-IoT systems. The UE has to start decoding
attempts for multiple possible repetitions which are limited
by a maximum number of repetitions, which is 2048 for NB
physical downlink shared channel (NPDSCH) and 128 for NB
physical uplink shared channel (NPUSCH), respectively. There
are three DL control channel indicator (CCI) formats specified
which have different bit positions for the DCI information [4],
and the number of repetitions is written into a field of the DCI.

The UE can therefore combine channel estimates over received
repetition copies in the channel coherent-time to improve the
CE accuracies and the data-transmission performance.
Due to the fluctuations in phase within receiver and transmitter oscillators and residual frequency offset (FO), there
are random phase noises presented in received signals, which
can severely degrade the mean-square-error (MSE) of channel
estimation (CE). The performance losses in relation to signalto-noise ratio (SNR) and bit error rate (BER) in OFDM
systems have been well-studied in e.g., [6]–[11]. With NBIoT systems considered, the problem is somewhat different
from what has been considered before. Firstly, as NB-IoT
systems only use 12 subcarriers for data-transmission, the
phase noises can be approximately assumed to be identical
through one PRB and yields a common phase error (CPE)
estimation problem [10]. Secondly, for low-end NB-IoT UEs,
the CE algorithm needs to have a low computational cost. For
instance, the maximum a posteriori CE algorithm in [12] and
a majorization-minimization based joint phase noise and CE
algorithm in [13] are over complex and consumes too much
power. Thirdly, the data-storage of the CE algorithm needs to
be kept small in order to save the memory.
In this work, we consider CE in the presence of phase noise
with repetition coded1 transmit data-blocks in NB-IoT systems, and derive the sequential minimum-mean-square-error
(MMSE) based CE using NB reference symbols (NRSs). The
proposed sequential MMSE estimator has a low-complexity
and only the correlation matrix needs to be stored and updated
over each iteration. Further, the proposed estimator can be
reviewed as an extension of a traditional sequential MMSE
estimator with taking the impact of phase noise into consideration, and has a structure of a Kalman-filter [14]. Numerical
results show that in the low the SNR regime, the proposed
estimator performs around 1 dB better than a traditional
sequential MMSE estimator that compensates the phase noise,
but do not take their estimation accuracies into consideration.
II. R ECEIVED S IGNAL M ODEL IN NB-I OT S YSTEMS
We consider an NB-IoT system with a single transmit
and receive antennas, and the NRS pattern is depicted in
Fig. 1. For stand-alone and guard-band deployments, no LTE
resource needs to be protected, and NPDCCH, NPDSCH or
NRS can utilize all the resource elements (REs) in one PRB
1 Note that, as we only consider the reference symbols based CE, the
proposed estimator also works for different transmit data-blocks as long as
the channels are in the coherent-time.

2

N

2 −1
j2πnk
1 X
s̃` [n]e− N
S̃` [k]= √
N
N

n=−
N
2

=

N
2

H` [k] X
N
N
n=−

2

−1

X

S` [m]ejφ` [n] e

j2πn(fe +m−k)
N

+ w̃[k]

m=− N
2
N
2

−1
X

=H` [k]S` [k]I` (0) +

m=− N
2

S` [m]I` [m]+ w̃[k],

(2)

,m6=k

|

{z

}

v[k]

where H` [k] is the fading channels on frequency domain, I` [k]
is a function of phase noises and FO, v[k] comprises both
the interference and noise. For each received repetition copy,
although the channel taps H` [k] can be assumed to be identical
as they are in the coherent-time, however, the term
N
2

I` (0) =

−1
X

ej (φ` [n]+

2πnfe
N

)

(3)

n=− N
2

changes due to the randomness of φ` [n]. As in NB-IoT, only
the middle 12 subcarriers are used for data-transmission, I` (0)
can be assumed to be identical for one received repetition copy.
Further, as from [10], [11], the module of I` (0) is very close
to one, and can be seen as a symbol rotation in the complex
plane. That is, we can approximate I` (0) by some ejφ` .
From (2), stacking the received signal corresponding to all
NPS symbols for a mth repetition copy in a vector form yields
S̃ m ≈ ejφm h

NPS Port 1
LTE CRS Port 0-3
LTE PDCCH

…...

…...

even-indexed slots

odd-indexed slots

One subframe

Fig. 1. The NRS pattern in one PRB (180 KHz) for stand-alone, guard-band
or in-band deployments of NB-IoT systems in LTE Rel. 13.

the term v m comprising both the interference and noise is
modeled as AWGN with a variance γ for all copies. The
random phase rotation φm is a real value, whose probability
density function (PDF), for analytical tractability, can be
assumed to a uniform distribution over [0, 2π). For generality,
we assume that there are K NRS in one repetition copy.
By removing S m from both sides in (4), we obtain a signal
model corresponding to the least-square (LS) estimates of the
channel vector as
r m ≈ ejφm h + ṽ m ,

2

−1

NPS Port 0

One PRB

pair. However, for in-band deployment, NPDCCH, NPDSCH
or NRS cannot be mapped to the REs taken by LTE cellspecific reference symbols (CRS) and control channels. NBIoT is designed to allow a UE to learn the deployment mode
and cell identity through initial acquisition.
With phase noise φ` [n] caused by oscillator fluctuations and
a residual FO fe normalized by the subcarrier frequency, the
time-domain baseband received signal at the nth sampling
time of the `th orthogonal-frequency-division-multiplexing
(OFDM) symbol can be expressed as


N/2−1
X
2πn(fe +k)
1
? h` [n] + w(n),
N
S` [k]ej
s̃` [n] = ejφ` [n] √
N k=−N/2
(1)
where “?” denotes the linear convolution, S` [k] is the transmit
symbol on the `th OFDM symbol and the kth subcarrier,
h[n] is discrete fading channel taps, w(n) is additive-whiteGaussian-noise (AWGN), and N is the Fast-Fourier-Transform
(FFT) size. An FFT operation is implemented at the UE to
obtain the received frequency domain signals,

S m + vm ,

(4)

where “ ” is the Hadamard product, S̃ m is the received signal
vector, S m is the transmitted NPS vector at the mth repetition
copy with zero-means and satisfying E[S m S †m ] = I. The
elements in the channel vector h are independent complex
Gaussian variables with zero-means and unit-variances which
can be assumed to be the same for a whole repetition period,

(5)

where the noise ṽ m has the same statistical properties as v m .
Given the received signal mode (5), we consider estimating
h in a sequential manner based on each received repetition
copy, which saves the required storage and also simplifies the
operations for low-end NB-IoT UEs.
III. T HE P ROPOSED S EQUENTIAL MMSE BASED CE
Without loss of generality, we let φ0 = 0 for the first received
copy, and denote ĥ as the CE obtained in a previous step. With
a newly received mth copy, the vector ejφm h can be estimated.
However, this estimate cannot be directly combined with ĥ
due to the phase rotation ejφm . An intuitive way is to use ĥ
to estimate φm first, and then obtain an update estimate ĥ after
compensating the phase rotation. This can be implemented by
a Kalman-filter based approach, but the optimal filtering needs
to be solved. Instead, we directly derive a sequential MMSE
estimator for h in he following.
The sequential MMSE estimator has the form [14]
Z
E(h|r) =
hp(h|r)dh
ZC Z
hp(r|h, φ)p(h)p(φ)
dhdφ
=
p(r)
C RZ Z
1
=
hp(r|h, φ)p(h)p(φ)dhdφ, (6)
p(r) C R
where the conditional PDF p(r|h, φ) with a given pair (h, φ)
and under SNR 1/γ is complex Gaussian and equal to


1
kr − h exp(jφ)k2
p(r|h, φ) =
exp
−
.
(7)
(πγ)K
γ
Assuming that ĥ is the estimate of h obtained with the previous received copy, the PDF p(h) is also complex Gaussian
and can be estimated by


1
† −1
p(h|ĥ) =
exp
−(
ĥ
−
h)
R
(
ĥ
−
h)
, (8)
(π)K det(R)

3

where R = E[hh† ] denotes the correlation matrix of h.
With (7) and (8), the integrals in (6) and p(r) are expressed
in (9) and (10), respectively. Therefore, the estimator (6) is
ξ1
E(h|r) = ,
ξ2

a newly received repetition copy. The MSE matrix

R̃ =

I+

R
γ

−1
,

(15)

(11)
and the factor ζ equals

where ξ1 and ξ2 are further derived in (12) and (13), respectively. Inserting them back into (11) yields a final form of the
proposed MMSE estimator, which is stated in Theorem 1.
Theorem 1. The sequential MMSE estimator of channel
vector h in the presence of random phase noise is


ζ
E(h|r) = R̃ ĥ + Rr ,
(14)
γ

ζ=

β1 r † R̃ĥ
,
β2 r † R̃ĥ

(16)

with β1 and β2 calculated as

where ĥ is the estimate of h obtained in the previous step
with R being its correlation matrix, 1/γ is the SNR, and r is




exp(jφ) †
β1 =
exp −jφ+2R
r R̃ĥ p(φ)dφ,
γ



ZR
exp(jφ) †
β2 =
exp 2R
r R̃ĥ p(φ)dφ.
γ
R
Z

(17)
(18)

Z Z
hp(r|h, φ)p(h)p(φ)dhdφ


Z Z


1
kr − h exp(jφ)k2
1
† −1
h
exp
−
exp
−(
ĥ
−
h)
R
(
ĥ
−
h)
dhdφ
=
(πγ)K (π)K det(R) C R
γ


†
1
1
1
=
exp − krk2 − ĥ R−1 ĥ
(πγ)K (π)K det(R)
γ




 
Z Z
†
1 †
1
×
h exp −h† R−1 + I h + 2R
r exp(jφ) + ĥ R−1 h
p(φ)dhdφ
γ
γ
C R
|
{z
}
C

R

(9)

ξ1

Z Z
p(r)=
p(r|h, φ)p(h)p(φ)dhdφ
C R


Z Z


1
1
kr − h exp(jφ)k2
exp
−
exp −(ĥ − h)† R−1 (ĥ − h) dhdφ
=
K
K
(πγ) (π) det(R) C R
γ


†
1
1
1
−1
2
exp
−
krk
−
ĥ
R
=
ĥ
(πγ)K (π)K det(R)
γ


 


Z Z
† −1
1 †
1
†
−1
r exp(jφ) + ĥ R
h
p(φ)dhdφ (10)
×
exp −h R + I h + 2R
γ
γ
C R

−1
|
{z
}

Z 
R−1 + γ1 I
ξ2
1


ξ1 =
r exp(−jφ) + R−1 ĥ
γ
R
det R−1 + γ1 I


−1 
!
† −1
1 †
1
1
−1
−1
× exp
r exp(jφ) + ĥ R
R + I
r exp(−jφ) + R ĥ
p(φ)dφ
γ
γ
γ

−1
!



−1
R−1 + γ1 I
−1
†
1 †
1
−1
−1

 exp
=
r γR + I
r exp ĥ I + R
R ĥ
γ
γ
det R−1 + γ1 I
(


−1 )!
Z 
1
exp(jφ) †
R
−1
×
r exp(−jφ) + R ĥ exp 2R
r I+
ĥ
p(φ)dφ.
(12)
γ
γ
γ
R

−1 
!
1
1
−1
−1


ξ2 =
exp
R + I
r exp(−jφ) + R ĥ
dφ
γ
γ
R
det R−1 + γ1 I
!



−1
−1
†
1
1 †
1
−1
−1

 exp
=
r γR + I
r exp ĥ I + R
R ĥ
γ
γ
det R−1 + γ1 I
(

−1 )!
Z
exp(jφ) †
R
r I+
ĥ
p(φ)dφ.
× exp 2R
γ
γ
R
1

Z



†
1 †
r exp(jφ) + ĥ R−1
γ

(13)

4

To obtain closed-from expressions for β1 and β2 , we assume
φ to be uniformly distributed in [0, 2π) which yields


2 †
β1 = I 1
r R̃ĥ ,
(19)
γ


2 †
(20)
r R̃ĥ ,
β2 = I 0
γ
using [15, Formulas 3.937.1 and 3.937.2], where I1 and I0 are
the modified Bessel functions of the first kind.
The variable ζ in (16) comprises two parts: the first part
is to estimate the phase noise with r † R̃ĥ |r † R̃ĥ|, and
the second part
the estimate by a factor
. isto scale down

I0 γ2 r † R̃ĥ , which is smaller than 1 and
I1 γ2 r † R̃ĥ
asymptotically equals 1 as SNR increases.
Now let’s consider two simple cases to illustrate the sequential MMSE estimator (14). The first case is R = I, that is, the
elements in h are identical and independently distributed (IID).
In this case, (14) becomes

−1 

1
ζ1
E(h|r) = 1 +
ĥ + r
γ
γ
γ
ζ1
=
ĥ +
r,
(21)
γ+1
γ+1
where
 † 
2|r ĥ|
I1 γ+1
r † ĥ

ζ1 = 
.
(22)
2|r † ĥ|
† ĥ
r
I0 γ+1
The estimator (21) is the traditional sequential MMSE estimator [14], but with an extra factor ζ1 which is introduced by
the random phase.
The second case is R = 1K×K , that is, the elements in h are
fully dependent and only one scalar h needs to be estimated.
Following the same derivations leading to (14), it can be shown
that the sequential MMSE estimator in this case is
γ
ζ2
E(h|r) =
ĥ +
r̃,
(23)
γ+1
γ+1
where r̃ is the mean value of r and


2K |r̃ ĥ|
I1
γ+1
r̃† ĥ

ζ2 = 
.
(24)
2K |r̃ ĥ|
† h̃
r̃
I0
γ+1
This is equivalent to the first case with a single observation r̃,
but both the variances of h and the noise are decreased by a
factor of K, which make an intuitive sense.
Using Theorem 1, after E(h|r) is computed, the correlation
matrix R needs to be updated at each step. By assuming that
ζ results in a perfect estimation of the phase and as β1 /β2 ≈ 1
at high SNR, R can be updated according to
−1

Rm
.
(25)
Rm+1 = Rm I +
γ
The initialization of R can be set to a heuristic value or estimated through a long-term averaging. With ĥm , the estimate
of phase noise φm of the mth repetition copy is obtained as
 †

φ̂m = ∠ ĥm r m ,
(26)

where “∠” takes the angle of the input variable.
Algorithm 1 The proposed sequential MMSE estimator.
Initialize m = 0 and ĥm = 0, and assume the channel
correlation matrix R0 and γ are known.
2: With the (m+1)th received repetition copy r m+1 , update
the parameters according to
−1

Rm
,
R̃m+1 = I +
γ


ζ
ĥm+1 = R̃m+1 ĥm + Rm r m+1 ,
γ
−1

Rm
.
Rm+1 = Rm I +
γ

1:

3:

where ζ is calculated according to (16).
Repeat Step 2 until the transmit data-block is successively
decoded or all repetition codes have been received.

Based on Theorem 1, the proposed sequential MMSE estimator is summarized in Algorithm 1. Note that, under the
case that there is no phase noise φ presented, we can set ζ = 1
and then Algorithm 1 is identical to the traditional sequential
MMSE estimator [14] for h.
IV. N UMERICAL R ESULTS
In this section, we evaluate MSE of the proposed sequential
MMSE estimator in comparison to the traditional sequential
MMSE estimator, which refers to an estimator that only
compensates the phase noise by setting ζ to
ζ=

r † R̃ĥ

,

(27)

r † R̃ĥ
in Algorithm 1, that is, setting β1 /β2 = 1 in (16) and ignores
the impact of β1 and β2 , i.e., the quality of the phase estimates.
In addition, we also evaluate the MSE under ideal cases when
there is no phase noise presented, in which cases, the proposed
sequential MMSE and traditional estimators are identical. The
MSE at the mth repetition copy is measured2 according to

N −1
2
1 X
MSE(m) =
ĥn,m exp(j φ̂n,m )−hn exp(jφn,m )
,
N K n=0
where we simulate N = 20000 channel realizations, and ĥn,m ,
φ̂n,m are estimates of hn and φn,m at the mth repetition copy.
In Fig. 2, we evaluate the MSE of the proposed and traditional sequential MMSE estimators, and the MSE obtained
when there is no phase noise present. We test under AWGN
channels with SNR equals to -4, -2 and 0 dB respectively.
As can be seen, the proposed sequential MMSE estimator
has around 1 dB gain over the traditional estimator without
considering β1 /β2 in (16) at SNR -4 dB. In Fig. 3, we repeat
the tests in Fig. 2 under ETU-3Hz channels with SNR equal
to -3, 0 and 3 dB, respectively. As can be seen, the proposed
estimator also outperforms the traditional estimator around 1
dB at SNR -3 dB.
2 The

MSE can also be equivalently calculated as the trace of Rm .

5

0

0
From top to down SNR=-4, -2 and 0 dB

-2

From top to down SNR=-3, 0 and 3 dB
-5

-6

MSE [dB]

MSE [dB]

-4

-8
-10

-10

-15

-12
-20
-14
-16

0

5

10

15

20

25

30

35

40

0

5

10

15

20

25

Numbe of received repetition copies

Numbe of received repetition copies

Fig. 2. The evaluation of MSE under AWGN channel.

Fig. 3. The evaluation of MSE under ETU-3Hz channel.

V. S UMMARY
We have considered channel estimation (CE) in narrowband Internet-of-Things (NB-IoT) systems in the presence of
random phase noise caused by the fluctuations of oscillators
and the residual frequency offset (FO). We have derived a
sequential minimum-mean-square-error (MMSE) estimator for
the CE, which updates the estimates in a sequential way with a
low computational cost and a small data storage. Moreover, we
show through simulations that, the proposed sequential MMSE
estimator performs 1 dB better in the low signal-to-noise ratio
(SNR) regime in terms of mean-square-error (MSE) than a
traditional sequential MMSE estimator that did not thoroughly
consider the impact of the phase noises.
R EFERENCES
[1] Y. P. E. Wang, X. Lin, A. Adhikary, A. Grövlen, Y. Sui, Y. Blankenship,
J. Bergman, and H. S. Razaghi, “A primer on 3GPP narrowband internet
of things,” IEEE Commun. Mag., vol. 53, no. 3, pp. 117-123, Mar. 2017.
[2] 3GPP TS 45.820, “Technical Specification Group GSM/EDGE Radio
Access Network; Cellular system support for ultra-low complexity and
low throughput Internet of Things (CIoT),” v13.1.0, Dec. 2015.
[3] 3GPP TS 36.201, “Evolved Universal Terrestrial Radio Access (EUTRA); Physical Layer General Description,” v13.3.0, Mar. 2017.

35

40

2

4

6

8
4
0

MSE [dB]

In both Fig. 2 and 3, when SNR increases, the proposed
sequential MMSE estimator starts to perform close to the traditional estimator, which is due to the fact that β1 /β2 quickly
converges to 1 as SNR increases. Further, as the number of
received repetition copies increases, the MSE renders errorfloors under both cases compared to the ideal cases, as a
consequence of the presence of phase noises.
As shown in Fig. 3, the MSE converges within 10∼20
repetition copies. In Fig. 4, we plot the MSE obtained at
the 20th received copy under ETU-3Hz channels and at
different SNR values. As can be seen, when SNR is less
than 0 dB, the proposed sequential MMSE estimator renders
considerable gains over the traditional estimator, with only
a minor complexity increment to compute the values β1 /β2
through e.g., look-up-table operations. We also investigate the
impact of different initializations of R0 . Setting R0 to the
ideal correlation matrix yields around 2 dB gains over the
case initializing R0 with an identity matrix. When there is no
phase noise present, the MSE gains are larger with ideal R0 .

30

-4
-8
-12
-16
-20
-24
-10

-8

-6

-4

-2

0

SNR [dB]

Fig. 4. The MSE measured at the 20th received copy under ETU-3Hz channel.
[4] 3GPP TS 36.212, “Evolved Universal Terrestrial Radio Access (EUTRA); Multiplexing and channel coding”, v14.2.0, Mar. 2017.
[5] S. Hu, H. Kröll, Q. Huang, and F. Rusek, “A Low-complexity channel
shortening receiver with diversity support for evolved 2G device,” IEEE
Int. Conf. on Commun. (ICC), Kuala Lumpur, May, 2016, pp. 1-7.
[6] T. Pollet, M. Bladel, and M. Moeneclaey, “BER sensitivity of OFDM
systems to carrier frequency offset and Wiener phase noise,” IEEE Trans.
on Commun., vol. 43, no.2, pp. 191-193, Feb. 1995.
[7] A. G. Armada and M. Calvo, “Phase noise and sub-carrier spacing
effects on the performance of an OFDM communication system,” IEEE
Commun. Lett., vol. 2, pp. 11-13, Jan. 1998.
[8] L. Tomba, “On the effect of Wiener phase noise in OFDM systems,”
IEEE Trans. on Commun., vol. 46, pp. 580-583, May 1998.
[9] M. S. El-Tanany, Y. Wu, and L. Hazy, “Analytical modeling and
simulation of phase noise interference in OFDM-based digital television
terrestrial broadcasting systems,” IEEE Trans. on Broadcast., vol. 47, pp.
20-31, Mar. 2001.
[10] D. Petrovic, W. Rave, and G. Fettweis, “Common phase error due to
phase noise in OFDM-estimation and suppression, ” in Proc. IEEE Int.
Symp. on Personal, Indoor and Mobile Radio Commun. (PIMRC), Sep.
2004, pp. 1901-1905.
[11] S. Wu and Y. Bar-Ness, “OFDM Systems in the presence of phase noise:
Consequences and solutions,” IEEE Trans. on Commun., vol. 52, no. 11,
pp. 1988-1996, Nov. 2004.
[12] D. D. Lin, R. A. Pacheco, T. J. Lim, and D. Hatzinakos, “Joint estimation
of channel response, frequency offset, and phase noise in OFDM,” IEEE
Trans. on Signal Process., vol. 54, no. 9, pp. 3542-3554, Sep. 2006.
[13] Z. Wang, P. Babu, and D. P. Palomar, “Effective low-complexity
optimization methods for joint phase noise and channel estimation in
OFDM,” IEEE Trans. on Signal Process., vol. 65, no. 12, pp. 32473260, Jun. 2017.
[14] S. M. Kay, Fundamentals of statistical signal processing, volume I:
Estimation theory, Prentice Hall signal processing series, 1993.
[15] A. Abramowitz and I. A. Stegun, Handbook of mathematical functions:
With formulas, graphs, and mathematical tables, the ninth reprint, Dover
Publications, INC., New York, Dec. 1972.

