# Nonlinear response of the irreversible work via generalized relaxation
functions

Pierre Nazé<EMAIL_ADDRESS>Universidade Estadual Paulista, 14800-090,
Araraquara, São Paulo, Brazil

###### Abstract

The nonlinear response of the excess work, when made via series expansion in
the parameter perturbation of the average thermodynamic work, requires
adjustments to agree with the Second Law of Thermodynamics. In this work, I
present a well-behaved nonlinear response of the irreversible work, based on
its well-known cumulant series expansion. From the generalization of the
fluctuation-dissipation relation derived from it, I define the terms of the
series expansion in the parameter perturbation of the irreversible work by the
terms of the cumulants. Since every cumulant depends on raw moments, I define
from them the generalized relaxation functions, whose arbitrary constants were
chosen guaranteeing the accomplishment of the Second Law of Thermodynamics. A
procedure to calculate the nonlinear response of the irreversible work is then
provided.

##  I Introduction

The theoretical calculation of the irreversible work is a hard task to
accomplish since the probabilistic distribution of the phase space is not
completely known for arbitrary driven processes. A partial solution for this
problem is considering approximations, like linear response theory, where the
parameter perturbation is slightly changed. In this case, the only requirement
is knowing the solution of the non-perturbed system. Higher orders of such
theory undoubtedly would be desired to have a better range of performance in
the parameter perturbation.

In Ref. [1] the development of the nonlinear response theory was made
considering thermally isolated systems. Unfortunately, the final result was
disappointing since an ad hoc procedure was necessary for the series to agree
with the Second Law of Thermodynamics. Indeed, without such intervention, the
excess work would diverge for large switching times. This is a sign that this
nonlinear response theory captures this usual behavior that occurs in time-
independent perturbation theory in classical mechanics.

In this work, I present a nonlinear response theory for systems performing
isothermal processes. In order to avoid any possible divergence for long
switching times, I adopt an alternative procedure to the calculation of the
irreversible work: the well-known expansion in cumulants [2]. In this manner,
the irreversible work can be expressed in terms involving raw moments of the
distribution, which are written in their turn in generalized relaxation
functions. To agree with the Second Law of Thermodynamics, the arbitrary
constants on which the generalized relaxation functions depend are chosen for
the moments to nullify for large switching times. This immediately implies
that the cumulants will have the same behavior, and so will the terms of the
irreversible work. This will avoid the diverging behavior observed for
thermally isolated systems. A procedure for the calculation of the
irreversible work using this nonlinear response theory is provided at the end.

##  II Preliminaries

Consider a classical system weakly coupled with a heat bath of temperature
$\beta^{-1}$. Their Hamiltonians are $\mathcal{H}_{\mathcal{S}}$ and
$\mathcal{H}_{\mathcal{B}}$, such that their sum are represented by
$\mathcal{H}_{0}$. A linear perturbation $\lambda(t)\mathcal{H}_{1}$ is added
to the composed system, where $\lambda(t)=\lambda_{0}+g(t/\tau)\delta\lambda$.
Here $\tau$ is the switching time of the process, $\delta\lambda$ the
parameter perturbation, and $g(0)=0$ and $g(\tau)=1$. Our interest is to
calculate the irreversible work

$W_{\rm irr}(\tau)=\langle W\rangle-\Delta F,$ (1)

where $\langle W\rangle$ is the average thermodynamic work

$\langle
W\rangle=\int_{0}^{\tau}\langle\mathcal{H}_{1}\rangle(t)\dot{\lambda}(t)dt,$
(2)

and $\Delta F$ is the difference of Helmholtz’s free energy between the final
and initial equilibrium state of the composed system. The average is
calculated in the probability distribution $\rho$, which satisfies Liouville’s
equation

$\frac{\partial\rho}{\partial t}=\mathcal{L}\rho,$ (3)

where
$\mathcal{L}(\bullet)=-\\{\bullet,\mathcal{H}_{0}+\lambda(t)\mathcal{H}_{1}\\}$
is called Liouville’s operator.

##  III Cumulant series expansion

It is well known that the irreversible work can be defined by the series
expansion [2]

$\beta W_{\rm irr}=\sum_{n=2}^{\infty}\beta\omega_{n},$ (4)

with

$\beta\omega_{n}=\frac{(-\beta)^{n}}{n!}\kappa_{n},$ (5)

where $\kappa_{n}$ are the $n$-th cumulant of the probability distribution of
the thermodynamic work. This cumulant series expansion agrees with the Second
Law of Thermodynamics. Indeed

$\beta W_{\rm
irr}=\sum_{n=2}^{\infty}\frac{[i(i\beta)]^{n}}{n!}\kappa_{n}=\ln{\phi(i\beta)}+\beta\langle
W\rangle,$ (6)

which leads to Jarzynski’s equality [2]

$\langle e^{-\beta W}\rangle=e^{-\beta\Delta F}.$ (7)

Observe that the expression (5) can be called a generalization of the
fluctuation-dissipation relation. The idea is that by using such an expansion
the diverging problem observed before for thermally isolated systems [1] will
disappear.

##  IV Nonlinear response of the irreversible work

To generate the nonlinear response of $W_{\rm irr}$, the cumulants need to be
expanded in the parameter perturbation of the process, and equal orders have
to be summed up. Therefore

$\omega_{n}=\sum_{m=0}^{\infty}\omega_{nm}=\sum_{m=0}^{\infty}\frac{(-\beta)^{n}}{n!}\kappa_{nm},$
(8)

where $\omega_{nm}$ and $\kappa_{nm}$ are the $m$-th term of the series
expansion in the parameter perturbation of $\omega_{n}$ and $\kappa_{n}$. The
$m$-th term of the series expansion of $W_{\rm irr}$ in the parameter
perturbation is

$\beta W_{\rm
irr}^{(m)}=\sum_{n=2}^{\infty}\omega_{nm}=\sum_{n=2}^{m}\frac{(-\beta)^{n}}{n!}\kappa_{nm},$
(9)

which can be called a generalization of the fluctuation-dissipation relation
in the context of nonlinear response theory. To illustrate the result with a
simple example, the linear response of $W_{\rm irr}$ is

$W_{\rm irr}^{(2)}=\frac{\beta}{2}\kappa_{22}=\frac{\beta}{2}\sigma^{2}_{2},$
(10)

where the term $\sigma^{2}_{2}$ is the variance of the work calculated until
its second order in the parameter perturbation. This result has already been
predicted [3].

##  V Generalized relaxation functions

The calculation of cumulants is not simple, and considering expansions in the
parameter perturbation requires much more attention. In the next paragraphs, I
will translate this procedure into the calculation of generalized relaxation
functions, as is done in linear response theory [4].

To start, cumulants can be defined in terms of the $i$-th raw moments
$\mu_{i}$ of the probability distribution

$\kappa_{n}=\sum_{i=1}^{n}(-1)^{i-1}(i-1)!\mathcal{B}_{n,i}(\mu_{1},...,\mu_{n-i+1}),$
(11)

where $\mathcal{B}_{n,i}$ are the partial Bell polynomials. In particular, the
$m$-th expansion term $\kappa_{nm}$ will be expressed as

$\kappa_{nm}=\sum_{i=1}^{n}\sum_{m_{1}+...+m_{(n-i+1)}=m}(-1)^{i-1}(i-1)!\mathcal{B}_{n,i}\left(\mu_{1m_{1}},...,\mu_{(n-i+1)m_{n-1+1}}\right).$
(12)

Now each raw moment has to be expanded in a series of the parameter
perturbation

${(-\beta)}^{n}\mu_{n}={(-\beta)}^{n}\langle W^{n}\rangle,$ (13)

which furnishes

${(-\beta)}^{n}\mu_{n}=\int_{0}^{\tau}...\int_{0}^{\tau}{(-\beta)}^{n}\left\langle\prod_{j=1}^{n}\mathcal{H}_{1}(t_{j})\right\rangle\dot{g}(t_{1})...\dot{g}(t_{n})dt_{1}...dt_{n}.$
(14)

The term on the average is calculated at different times. In this case, I
consider that the probability distribution has the form

$\rho=\rho(z_{t_{1}},...,z_{t_{n}},t_{1},...,t_{n}),$ (15)

such that

$\int_{\Gamma_{1}}...\int_{\Gamma_{n}}\rho(z_{t_{1}},...,z_{t_{n}},t_{1},...,t_{n})dz_{t_{1}}...dz_{t_{n}}=1.$
(16)

Demanding that the probability distribution should be invariant over time, we
have

$d\rho=\sum_{i=1}^{n}\left(\frac{d\rho}{dt_{i}}\right)dt_{i}=0,$ (17)

meaning

$\frac{d\rho}{dt_{i}}=0\quad\rightarrow\quad\frac{\partial\rho}{\partial
t_{i}}=-\mathcal{L}_{i}\rho$ (18)

The solution must satisfy each of Liouville’s equations, requiring that the
probability distribution starts from the canonical ensemble. Let us try the
most trivial solution

$\rho(z_{t_{1}},...,z_{t_{n}},t_{1},...,t_{n})=\rho_{1}(z_{t_{1}},t_{1})...\rho_{n}(z_{t_{n}},t_{n}),$
(19)

with

$\rho(z_{0},...,z_{0},0,...,0)=\rho_{c}(z_{0}).$ (20)

Under such circumstances, the moment will be

${(-\beta)}^{n}\mu_{n}=\int_{0}^{\tau}...\int_{0}^{\tau}{(-\beta)}^{n}\prod_{j=1}^{n}\left\langle\mathcal{H}_{1}(t_{j})\right\rangle_{j}\dot{g}(t_{1})...\dot{g}(t_{n})dt_{1}...dt_{n},$
(21)

where the index $j$ in the average indicates the $j$-th probability
distribution of $\rho$. Using Liouville’s theorem, one has

${(-\beta)}^{n}\mu_{n}=\int_{0}^{\tau}...\int_{0}^{\tau}{(-\beta)}^{n}\prod_{j=1}^{n}\left\langle\mathcal{H}_{1}(0)\right\rangle_{j}\dot{g}(t_{1})...\dot{g}(t_{n})dt_{1}...dt_{n}.$
(22)

Applying now the series expansion in the parameter perturbation

${(-\beta)}^{n}\mu_{n}=\int_{0}^{\tau}...\int_{0}^{\tau}{(-\beta)}^{n}\prod_{j=1}^{n}\sum_{i=0}^{\infty}\left\langle\mathcal{H}_{1}(0)\right\rangle_{ji}\dot{g}(t_{1})...\dot{g}(t_{n})dt_{1}...dt_{n}.$
(23)

one can rewrite in the following form

${(-\beta)}^{n}\mu_{n}=\sum_{i_{1}=0}^{\infty}...\sum_{i_{n}=0}^{\infty}\int_{0}^{\tau}...\int_{0}^{\tau}{(-\beta)}^{n}\prod_{j=1}^{n}\left\langle\mathcal{H}_{1}(0)\right\rangle_{ji_{j}}\Bigg{|}_{i_{1},...,i_{n}}\dot{g}(t_{1})...\dot{g}(t_{n})dt_{1}...dt_{n}.$
(24)

The $m$-th term of the $n$-th moment will be

${(-\beta)}^{n}\mu_{nm}=\sum_{n+i_{1}+...+i_{n}=m}\int_{0}^{\tau}...\int_{0}^{\tau}{(-\beta)}^{n}\prod_{j=1}^{n}\left\langle\mathcal{H}_{1}(0)\right\rangle_{ji_{j}}\Bigg{|}_{i_{1},...,i_{n}}\dot{g}(t_{1})...\dot{g}(t_{n})dt_{1}...dt_{n}.$
(25)

To each one of those terms can be put in terms of generalized relaxation
functions. If all $i_{j}=0$, one has

${(-\beta)}^{n}\prod_{j=1}^{n}\left\langle\mathcal{H}_{1}(0)\right\rangle_{j0}=\beta[\Psi_{n0}(t_{1},...,t_{n})+\mathcal{C}_{n0}],$
(26)

where

$-\Psi_{n0}(t_{1},...,t_{n})={(-\beta)}^{n-1}\left[\left\langle\prod_{j=1}^{n}e^{\mathcal{L}_{0}t_{j}}\mathcal{H}_{1}(0)\right\rangle_{0}-\mathcal{C}_{n0}\right],$
(27) $\mathcal{C}_{n0}=\langle\mathcal{H}_{1}(0)\rangle_{0}^{n}.$ (28)

If some $i_{j}\geq 1$, the calculation becomes more involved. Indeed

${(-\beta)}^{n}\prod_{j=1}^{n}\left\langle\mathcal{H}_{1}(0)\right\rangle_{ji_{j}}\Bigg{|}_{i_{1},...,i_{n}}=\left\langle(-\beta)^{n}\prod_{j=1}^{n}\mathcal{K}_{ji_{j}}(\Delta_{i_{j}},...,\Delta_{1})\mathcal{H}_{1}(0)\right\rangle_{0}.$
(29)

where the operator acts as

$\mathcal{K}_{ji_{j}}(\Delta_{i_{j}},...,\Delta_{1})(\bullet)=(-1)^{i_{j}}\prod_{k=1}^{i_{j}}\int^{s_{k}}_{0}\mathcal{L}_{1}(0)e^{-\mathcal{L}_{0}\Delta_{k}}(\bullet)g(s_{k-1})ds_{k-1},$
(30)

with $\Delta_{k}=s_{k}-s_{k-1}$ and $\Delta_{i_{j}}=t_{j}-s_{i_{j}-1}$.
Observe that

$\mathcal{L}_{1}(0)=\frac{d}{ds_{k-1}},$ (31)

from which, performing an integration by parts, one has

$\mathcal{K}_{ji_{j}}(\Delta_{i_{j}},...,\Delta_{1})(\bullet)=(-1)^{i_{j}}\prod_{k=1}^{i_{j}}\left[(\bullet)g(s_{k})-\int^{s_{k}}_{0}[e^{-\mathcal{L}_{0}\Delta_{k}}(\bullet)+\mathcal{C}_{jk}]\dot{g}(s_{k-1})ds_{k-1}\right],$
(32)

or also

$\mathcal{K}_{ji_{j}}(\Delta_{i_{j}},...,\Delta_{1})(\bullet)=\prod_{k=1}^{i_{j}}\left[\int^{s_{k}}_{0}[(e^{-\mathcal{L}_{0}\Delta_{k}}-1)(\bullet)+\mathcal{C}_{jk}]\dot{g}(s_{k-1})ds_{k-1}\right].$
(33)

Returning to Eq. (29)

${(-\beta)}^{n}\prod_{j=1}^{n}\left\langle\mathcal{H}_{1}(0)\right\rangle_{ji_{j}}\Bigg{|}_{i_{1},...,i_{n}}=\left\langle(-\beta)^{n}\prod_{j=1}^{n}\prod_{k=1}^{i_{j}}\left[\int^{s_{k}}_{0}[(e^{-\mathcal{L}_{0}\Delta_{k}}-1)\mathcal{H}_{1}(0)+\mathcal{C}_{jk}]\dot{g}(s_{k-1})ds_{k-1}\right]\right\rangle_{0},$
(34)

which is the same as

${(-\beta)}^{n}\prod_{j=1}^{n}\left\langle\mathcal{H}_{1}(0)\right\rangle_{ji_{j}}\Bigg{|}_{i_{1},...,i_{n}}=\left\langle(-\beta)^{n}\prod_{j=1}^{n}\int^{t_{j}}_{0}...\int^{s_{1}}_{0}\prod_{k=1}^{i_{j}}\left[[(e^{-\mathcal{L}_{0}\Delta_{k}}-1)\mathcal{H}_{1}(0)+\mathcal{C}_{jk}]\right]\dot{g}(s_{i_{j}-1})...\dot{g}(s_{0})ds_{i_{j}-1}...ds_{0}\right\rangle_{0}.$
(35)

Observe that the property of switching the product and integrals can be
performed one more time, leading to

${(-\beta)}^{n}\prod_{j=1}^{n}\left\langle\mathcal{H}_{1}(0)\right\rangle_{ji_{j}}\Bigg{|}_{i_{1},...,i_{n}}=\left(\prod_{j=1}^{n}\int^{t_{j}}_{0}...\int^{s_{1}}_{0}\right)\beta\Psi_{ni_{1}...i_{n}}(\\{\Delta_{1}\\}_{j},...,\\{\Delta_{n}\\}_{j})\left(\prod_{j=1}^{n}\dot{g}(s_{i_{j}-1})...\dot{g}(s_{0})ds_{i_{j}-1}...ds_{0}\right),$
(36)

where

$\beta\Psi_{ni_{1}...i_{n}}(\\{\Delta_{1}\\}_{j},...,\\{\Delta_{n}\\}_{j})=(-\beta)^{n}\left\langle\prod_{j=1}^{n}\prod_{k=1}^{i_{j}}[(e^{-\mathcal{L}_{0}\Delta_{k}}-1)\mathcal{H}_{1}(0)+\mathcal{C}_{jk}]\right\rangle_{0},$
(37)

The relaxation function
$\Psi_{n\\{i_{1}...i_{n}\\}}(\\{\Delta_{1}\\}_{j},...,\\{\Delta_{n}\\}_{j})$
is the sum over the permutations of the set $\\{i_{1},...,i_{n}\\}$ of
$\Psi_{ni_{1}...i_{n}}(\\{\Delta_{1}\\},...,\\{\Delta_{n}\\})$

$\Psi_{n\\{i_{1}...i_{n}\\}}(\\{\Delta_{1}\\}_{j},...,\\{\Delta_{n}\\}_{j})=\sum_{n+i_{1}+...+i_{n}=m}\Psi_{ni_{1}...i_{n}}(\\{\Delta_{1}\\}_{j},...,\\{\Delta_{n}\\}_{j}).$
(38)

To determine the terms $\mathcal{C}_{jk}$, I use the criterion that the
generalized relaxation functions should nullify when
$\Delta_{j}\rightarrow\infty$. This will guarantee that at each expansion the
moment will nullify at infinity, and also each cumulant and term of the
irreversible work. In that limit, the canonical ensemble decorrelates at each
random variable. Therefore

$\beta\Psi_{ni_{1}...i_{n}}(\\{\infty\\}_{j},...,\\{\infty\\}_{j})=(-\beta)^{n}\prod_{j=1}^{n}\prod_{k=1}^{i_{j}}\left\langle[(e^{-\mathcal{L}_{0}\Delta_{k}}-1)\mathcal{H}_{1}(0)+\mathcal{C}_{jk}]\right\rangle_{0}=0\quad\rightarrow\quad\mathcal{C}_{jk}=0.$
(39)

Therefore, the relaxation function is the sum over the permutations of the set
$\\{i_{1},...,i_{n}\\}$ of the following function

$\beta\Psi_{ni_{1}...i_{n}}(\\{\Delta_{1}\\}_{j},...,\\{\Delta_{n}\\}_{j})=(-\beta)^{n}\left\langle\prod_{j=1}^{n}\prod_{k=1}^{i_{j}}[(e^{-\mathcal{L}_{0}\Delta_{k}}-1)\mathcal{H}_{1}(0)]\right\rangle_{0}.$
(40)

##  VI Examples

The result seems to be consistent. Indeed, we recover linear response case,
where $n=2$ and $i_{1},i_{2}=0$,

$\Psi_{2\\{00\\}}(t_{1},t_{2})=\beta\left\langle\mathcal{H}_{1}(t_{1})\mathcal{H}_{1}(t_{2})\right\rangle_{0}-\beta\left\langle\mathcal{H}_{1}(0)\right\rangle_{0}^{2}.$
(41)

Considering $n=3$ and $i_{1},i_{2},i_{3}=0$ one has

$\Psi_{3\\{000\\}}(t_{1},t_{2},t_{3})=-\beta^{2}\left\langle\mathcal{H}_{1}(t_{1})\mathcal{H}_{1}(t_{2})\mathcal{H}_{1}(t_{3})\right\rangle_{0}+\beta^{2}\left\langle\mathcal{H}_{1}(0)\right\rangle_{0}^{3}.$
(42)

For its turn, when $n=2$ and $i_{1},i_{2}=1$, one has

$\Psi_{2\\{11\\}}(t_{1},t_{2})=\beta\left\langle(\mathcal{H}_{1}(t_{1})-\mathcal{H}_{1}(0))(\mathcal{H}_{1}(t_{2})-\mathcal{H}_{1}(0))\right\rangle_{0},$
(43)

Also, when $n=3$ and $i_{1},i_{2},i_{3}=1$, one has

$\Psi_{3\\{111\\}}(t_{1},t_{2},t_{3})=\beta^{2}\left\langle(\mathcal{H}_{1}(t_{1})-\mathcal{H}_{1}(0))(\mathcal{H}_{1}(t_{2})-\mathcal{H}_{1}(0))(\mathcal{H}_{1}(t_{3})-\mathcal{H}_{1}(0))\right\rangle_{0},$
(44)

Observe now for $n=2$ with $i_{1}=0$ and $i_{2}=1$

$\Psi_{2\\{01\\}}(t_{1},t_{2})=\beta\left\langle\mathcal{H}_{1}(t_{1})(\mathcal{H}_{1}(t_{2})-\mathcal{H}_{1}(0))\right\rangle_{0}+\beta\left\langle\mathcal{H}_{1}(t_{2})(\mathcal{H}_{1}(t_{1})-\mathcal{H}_{1}(0))\right\rangle_{0},$
(45)

Finally, consider $n=2$ with $i_{1}=0$ and $i_{2}=2$. One has

$\Psi_{2\\{01\\}}(t_{1},\\{s_{1},s_{2}\\})=\beta\left\langle\mathcal{H}_{1}(t_{1})(\mathcal{H}_{1}(s_{1})-\mathcal{H}_{1}(0))(\mathcal{H}_{1}(s_{2})-\mathcal{H}_{1}(0))\right\rangle_{0}+\beta\left\langle\mathcal{H}_{1}(s_{1})(\mathcal{H}_{1}(t_{1})-\mathcal{H}_{1}(0))(\mathcal{H}_{1}(t_{2})-\mathcal{H}_{1}(0))\right\rangle_{0},$
(46)

Observe now the role of each index in the generalized relaxation function: the
index $n$ indicates the number of terms that are multiplied in the average,
while the index $i_{j}$ indicates the number of differences
$(\mathcal{H}_{1}(t)-\mathcal{H}_{1}(0))$ involved in each one of those terms,
except when $i_{j}=0$, where it becomes $\mathcal{H}_{1}(t)$. The sum over the
permutations of $\\{i_{1},...,i_{n}\\}$ indicates that the generalized
relaxation functions are symmetric under the change of the vectors
$\\{\Delta_{j}\\}$. Also, when there is $i_{j}\neq 0$, the constant equals
zero. The sum $n+i_{1}+...+i_{n}$ indicates the order in the perturbation in
the $n$-th moment.

Observe that the calculation of the relaxation function is nothing more than a
linear combination of the terms

$\mathcal{R}_{uv}=\left\langle\left(\prod_{i=1}^{u}\mathcal{H}_{1}(t_{i})\right)\left(\prod_{i=1}^{v}\mathcal{H}_{1}(0)\right)\right\rangle_{0},$
(47)

with $u,v\geq 0$. These terms produce a relaxation matrix $\mathcal{R}(m)$ of
order $m\times m$ for an approximation until the $m$-th order in the parameter
perturbation. Because of properties from the average, the number of effective
terms to be calculated are $m^{2}-m$.

##  VII Procedure to calculate the approximation of the irreversible work

The procedure to calculate the approximation of the irreversible work until
its $n$-th order in the parameter perturbation is

  * •

Choose $n\geq 2$;

  * •

Consider $m=2$;

  * •

Identify the cumulants $\kappa_{l}$, with $2\leq l\leq m$;

  * •

Identify the moments $\mu_{p}$, with $1\leq p\leq l$, that compose the
cumulants $\kappa_{l}$;

  * •

Calculate the remaining terms of relaxation matrix $\mathcal{R}(m)$;

  * •

Calculate from $\mathcal{R}_{uv}$ all $\Psi_{p\\{i_{1},...,i_{p}\\}}$, with
$p+i_{1}+...+i_{p}=m$;

  * •

Use Eq. (25) to calculate $\mu_{pm}$;

  * •

Use Eq. (12) to calculate $\kappa_{lm}$;

  * •

Use Eq. (9) to calculate to $m$-th order of the irreversible work;

  * •

If $m=n$ is true, go to the next step, otherwise repeat the process from the
third step considering $m+1$;

  * •

Sum up all the $n-1$ terms calculated before to find the approximation until
$n$-th order.

##  VIII Final remarks

On this work, I developed a nonlinear response theory for the irreversible
work. The use of cumulants series expansion is an attempt to avoid diverging
behavior in this case, since it establishes generalized relaxation functions
that nullify at infinity by the choice of constants that provide such
behavior. Providing however an example of a system performing an isothermal
process until higher orders to verify such behavior is a herculeous work. A
new version of this work will be shared when such an example is calculated.

## References

  * Nazé and Bonança [2022] P. Nazé and M. V. S. Bonança, Journal of Statistical Physics 186, 23 (2022).
  * Jarzynski [1997] C. Jarzynski, Phys. Rev. Lett. 78, 2690 (1997).
  * Nazé [2023] P. Nazé, arXiv preprint arXiv:2304.11965 (2023).
  * Kubo  _et al._ [2012] R. Kubo, M. Toda, and N. Hashitsume,  _Statistical physics II: nonequilibrium statistical mechanics_, Vol. 31 (Springer Science & Business Media, 2012).


# Machine-Learning-based High-resolution DOA Measurement and Robust DM for
Hybrid Analog-Digital Massive MIMO Transceiver

Feng Shu, , Ling Xu, Jinyong Lin, Jingsong Hu,
Linqing Gui, Jun Li, and Jiangzhou Wang This work was supported in part by the
National Natural Science Foundation of China (Nos. 61472190, and
61501238).Feng Shu, Ling Xu, Jinsong Hu, Yijin Zhang, Linqing Gui, and Jun Li
are with School of Electronic and Optical Engineering, Nanjing University of
Science and Technology, 210094, CHINA. E-mail:{shufeng,
<EMAIL_ADDRESS>Shu,and Jinsong Hu are the College of Computer and
Information Sciences, Fujian Agriculture and Forestry University, Fuzhou
350002, China (e-mail: shufeng0101@163.com).Jinyong Lin is with Beijing
Aerospace Automatic Control Institute, Beijing 100854, China (e-mail:
ljiny3771@sina.com)Jiangzhou Wang is with the School of Engineering and
Digital Arts, University of Kent, Canterbury CT2 7NT, U.K. Email:
<EMAIL_ADDRESS>
###### Abstract

At hybrid analog-digital (HAD) transceiver, an improved HAD rotational
invariance techniques (ESPRIT), called I-HAD-ESPRIT, is proposed to measure
the direction of arrival (DOA) of desired user, where the phase ambiguity due
to HAD structure is addressed successfully. Subsequently, a machine-learning
(ML) framework is proposed to improve the precision of measuring DOA. In the
training stage, the HAD transceiver works as a receiver and repeatedly
measures the values of DOA via I-HAD-ESPRIT to form a slightly large training
data set (TDS) of DOA. From TDS, we find that the probability density function
(PDF) of DOA measurement error (DOAME) is approximated as a Gaussian
distribution by the histogram method in ML. This TDS is used to learn the mean
of DOA and the variance of DOAME, which are utilized to infer their values and
improve their precisions in the real-time stage. The HAD transceiver rapidly
measures the real-time value of DOA some times to form a relatively small
real-time set (RTS), which is used to learn the real-time mean and variance of
DOA/ DOAME. Then, three weight combiners are proposed to combine the-maximum-
likelihood-learning outputs of TDS and RTS. Their weight factors depend
intimately on the numbers of elements in TDS and RTS, and signal-to-noise
ratios during the two stages. Using the mean and variance of DOA/DOAME, their
PDFs can be given directly, and we propose a robust beamformer for directional
modulation (DM) transmitter with HAD by fully exploiting the PDF of DOA/DOAME,
especially a robust analog beamformer on RF chain. Simulation results show
that: 1) The proposed I-HAD-ESPRIT can achieve the HAD CRLB; 2) The proposed
ML framework performs much better than the corresponding real-time one without
training stage, 3) The proposed robust DM transmitter can perform better than
the corresponding non-robust ones in terms of bit error rate and secrecy rate.

###### Index Terms:

hybrid analog and digital, ESPRIT, statistical learning, DM, robust precoder

##  I Introduction

Direction of arrival (DOA) measurement, also known as direction finding,
belongs to one of the most important aspects in array signal processing [1, 2,
3, 4, 5]. It is also a very old research field and has a long and successful
history[6, 7, 8, 9]. In the recent years, massive MIMO promises to provide a
ultra-high-resolution DOA estimation [10]. Therefore, the research field of
estimating DOA regains its new life. The need for high resolution DOA
estimation always exists because of its wide applications in many practical
scenarios such as radar, sonar, seismic exploration, electronic
countermeasures and radio astronomy[11]. More importantly, with the rapid
development of some new technologies and appearance of new application
scenarios such as directional modulation(DM), device to device(D2D), satellite
communications, internet of things (IoT), uninhabited aerial vehicle(UAV), 5G
and beyond mobile works like millimeter wave communications, in our view, the
demand for DOA estimation will emerge as an explosive growth in the coming
future.

Classic and representative DOA estimation methods are categorized into two
main kinds: conventional and subspace methods. Conventional methods, also
called classical methods, first compute a spatial spectrum and then find DOAs
by numerical search of local maxima of the spectrum. Typical Methods falling
in this class are Capon and Bartlett methods. But these methods suffer from
angular resolution loss. Alternatively, two classic high-angular-resolution
subspace methods such as MUSIC (Multiple signal classification) and ESPRIT
(Estimation of signal parameters via rotation invariance techniques)
algorithms are more popular and achieve wide applications[12, 13].

Recently, compressed sensing(CS) was also applied to improve the precision of
DOA estimation. In [14], the authors developed an effective sparse
representation $\mathrm{L1-SVD}$ algorithm, which first used singular value
decomposition(SVD) of data matrix to make a reduction of dimensionality and
then the problem of DOA estimation was casted a second-order cone(SOC)
optimization. In [15], a mixed $l_{2,0}$ approximation DOA algorithm was
proposed to solve the joint-sparse recovery problem, which could distinguish
closely spaced and highly correlated sources without knowing the number of
emitters in advance.

The authors in [10] applied deep learning into massive MIMO system to conduct
DOA estimate. Specifically, training data set (TDS) is formed by using
randomly generated input and output data, and the weight parameters of deep
neural network(DNN) is trained to learn the statistic property of angle
domain. Once DNN is trained well, DNN can output the real-time value of DOA
rapidly. In [16], an adaptive DOA estimation method, based on convolutional
neural network(CNN) with long short term memory, was developed. In [17], the
authors proposed a CNN-based DOA estimation method for broad band-width signal
with the phase of short-time Fourier transform coefficients as inputs of CNN
to implement training process.

However, all these DOA related works consider the fully digital (FD)
structure, in which each antenna is equipped with one dedicated RF chain. If
the number of antennas goes to large-scale, it will lead to a significant
increasement in circuit cost and energy consumption. Therefore, a hybrid
analog-to-digital (HAD) MIMO structure is an excellent alternative solution to
FD due to its low circuit cost and high energy efficiency. The authors in [18]
developed a low-complexity HAD Root-MUSIC algorithm in MIMO receive array with
phase alignment, which could achieve the CRLB of HAD.

Compared to ESPRIT-based algorithms, the MUSIC-based methods are more
accurate, and stable especially in the high SNR scenarios[19]. But, the need
for gratuitous subarray calibration information brings more inconveniences for
MUSIC-based ones owing to that the phase drifts will reduce the reliability of
received data[20]. More important, for Root-MUSIC-based methods, to resolve
multiple emitters, the receiver should know the number of incident waves in
advance while ESPRIT doesn’t need such a number and is even capable of
estimating this number [13].

In this paper, to avoid the disadvantages of MUSIC-based algorithms, an
ESPRIT-based DOA estimator is proposed for a HAD receiver. Here, the phase
ambiguity is addressed by maximizing the output receive power, which is
similar to the Root-MUSIC in [18]. However, the kind of methods will face a
performance loss. With the help of machine learning framework, we further
improve the performance of this kind of method. Subsequently, to obtain more
precise variance, we make a further investigation of the statistical feature
of DOAMEs and conduct a weight combination of the DOA means and variances of
TDS and new on-line data set or real-time set (RTS). Since we complete the
density estimation of DOA, we propose a robust beamforming scheme for DM using
HAD structure. Before we summarize our main contributions, let us review the
literature concerning DM as follows.

As a secure physical-layer transmission technology, DM need to know the angle
direction information of desired user in advance, which can be obtained by the
preceding methods. The existing research works about DM synthesis can be
divided into two categories. The first one was based on RF front-end[21, 22,
23]. In [21], the authors developed a near-field direct antenna modulation
technique using a large amount of reflectors and switches to adjust the
amplitude and phase of signal, thus distorting signal constellation severely.
In [22], the authors proposed a phased-array-based DM, which produced phase
and amplitude of each symbol in desired direction by shifting phase of each
array element. Similarly, in [23], an improved enhancing direct antenna
modulation was presented by increasing element spacing to improve the
direction error rate in desired direction. The second way of DM synthesis is
implemented in baseband with the aid of artificial noise (AN) to disturb the
constellation pattern along undesired directions. In [24], the orthogonal
vector method with the help of AN was proposed to synthesize the waveform of
DM. In [25], a closed-form expression of null-space-projection(NSP) beamformer
was derived and a robust DM synthesis was proposed using the rule of
conditional minimum mean square error for single-desired-user scenario.
Furthermore, in [26, 27], the concept of robust DM is extended to two new
scenarios: multi-beaming broadcasting and MU-MIMO, respectively. In [28, 29],
the relay-aided DM is exploited to enhance the physical layer security. In
[28, 30, 31], the authors proposed a new secure concept, called secure precise
wireless transmission (SPWT). The former achieved the SPWT with the help of
multi-relay while the latter adopted a random frequency diverse array plus DM
and random subcarrier selection plus DM with the help of AN, respectively. In
[32], the optimal power allocation strategy between confidential message and
AN was optimized to maximize the secrecy rate (SR) of DM system. An
alternative iteration method to maximize SR by resorting to GPI was proposed
in [33]. Additionally, DM is utilized to harvest energy by its directional
property [29] due to its excellent directive property. In particular, in a
massive MIMO scenario, an extremely high energy-efficiency and security for
harvesting energy and transmitting confidential messages can be achieved. From
this aspect, we will predict the fact that DM will become a dramatically
important secure wireless transmission way in the future wireless networks.

However, all the above research works concerning DM focus on fully-digital
transmitter, if a massive transmit antenna array is adopted, a hybrid analog
and digital DM transmitter is a perfect choice of striking a good balance by
taking circuit cost, energy consumption, and computation complexity into
account. Traditional research works pertaining HAD architectures focus on the
design of two beamforming matrices including digital precoder and analog
precoder at transmitter. In [34], the authors first introduced HAD and
verified the number of RF chains needed to realize the FD precoding in single-
data-stream scenario. In [35, 36], the HAD is generalized to multi-user MIMO
scenarios. In [37], an effective alternating minimization algorithms was
proposed for both fully-connected and sub-connected HAD architectures. In
[38], a energy-efficient precoder based on successive interference cancelation
was conceived for partially-connected architecture. The authors in [39]
considered security in HAD systems by maximizing the projection between FD
precoder and the hybrid precoder. In [40], an overview of hybrid precoding
techniques especially in mmWave MIMO systems was presented. However, all the
preceding works optimized the two beamformers: analog and digital, without AN
projection.

To the best of our knowledge, there is still no scheme of combining HAD and DM
to provide a robust secure transmission with low-complexity and low cost in
line-of-sight (LoS) channels. Therefore, in this paper, a robust HAD plus DM
transmitter is presented with the aid of AN to achieve a robust and secure
physical-layer transmission by exploiting the PDF of measured DOA. Our main
contributions are summarized as follows:

  1. 1.

A ESPRIT DOA measurement method for a sub-connected receiver with HAD
structure is proposed to estimate directions of desired and undesired users,
respectively. By utilizing the rotational invariance among each subarray in
sub-connected hybrid architecture, the ESPRIT is well applicable to hybrid
analog and digital sub-connected receiver. And the phase ambiguity due to
antenna subarray structure is addressed by maximizing the output power over a
set of multiple potential solutions induced by periodicity of phase.
Subsequently, after applying the proposed ESPRIT repeatedly to form a set of
measured values of DOA, via histogram method in machine learning (ML), we find
that the density distribution of DOA or its measurement error (DOAME) can be
approximated as a Gaussian distribution, which significantly simplifies the
following ML-based improved method of DOA estimation.

  2. 2.

To improve the estimate performance of ESPRIT, a ML-based framework is
proposed. In this framework, in the first stage, called training or offline
stage, ESPRIT is repeatedly used to form a TDS, which will provide a high-
precision estimate for the mean of DOA and variance of DOAME. Here, the number
of elements in TDS can be chosen to sufficiently large in order to achieve a
high-performance predictor for variance of DOAME and mean of DOA in the next
stage. In the second stage, called real-time or new-data, the value of DOA are
rapidly estimated by several snapshots, and yields the mean of real-time DOA
and a coarse estimate of DOAME variances by using maximum likelihood learning.
Finally, the outcomes from the above stages are combined to form new outputs
by a weight method, where the weighted factors rely heavily on the receive
SNR, and the numbers of TDS and RTS. Now, we can directly have the probability
density function (PDF) of DOAME. The above process actually is the density
estimation. From simulation results, it follows that, compared to the
corresponding method without training stage, the proposed ML method can
dramatically enhance the performance of DOA and DOAME variance in terms of
variance. This paves a way to design a robust beamforming method of DM for
HAD-MIMO structure.

  3. 3.

After learning the mean values of DOAs and PDFs of DOAMEs of desired user and
eavesdropper, we can start to design a robust beamformer for DM in HAD MIMO
system, which is made up of three beamforming vectors/matrices required to be
devised: digital precoding vector of confidential message in baseband, digital
AN projection matrix in baseband, and analog beamforming matrix on RF chain.
The first two beamformers are optimized via minimizing the Euclidean distance
between the HAD beamforming matrix and corresponding FD precoder,
respectively. Our emphasis is mainly on the last analog beamforming matrix. By
making use of conditional expectation, a robust analog beamforming matrix is
proposed. Simulation results reveal that, compared to non-robust beamforming
methods, our proposed robust hybrid method performs much better in terms of
both SR and BER.

The remainder of this paper is organized as follows. Section II describes the
system models of HAD receiver to learn DOA and HAD DM transmitter to achieve a
robust secure transmission. In Section III, the ESPRIT DOA estimation using
HAD structure is proposed, the density of DOAMEs is approximately attained by
histogram method, and the TDSs of DOA measurements are constructed. Finally, a
weight combination between the outputs of TDS and RTS is to yield a more
precise mean and variance of DOA. In Section IV, based on conditional
expectation, a robust secure beamformer is proposed for DM HAD transmitter
with emphasis on analog part. Simulation results will be presented in Section
V. Finally, we draw our conclusions in Section VI.

Notations: throughout the paper, matrices, vectors, and scalars are denoted by
letters of bold upper case, bold lower case, and lower case, respectively.
$(\cdot)^{T}$, $(\cdot)^{*}$ and $(\cdot)^{H}$ denote transpose, conjugate,
and conjugate transpose, respectively. $\parallel\cdot\parallel_{2}$ and
$\parallel\cdot\parallel_{F}$ denote the $l_{2}$ norm of a vector and
Frobenius norm of a matrix, respectively. $\text{tr}(\cdot)$ and $vec(\cdot)$
are matrix trace and matrix vectorization. $\otimes$ denotes the Kronecker
products between two matrices.

##  II System Model

(a) Hybrid DOA estimation (b) Hybrid DM precoder

Figure 1: Block diagram of DOA measurement and DM for HAD transceiver

In Fig. 1, the schematic diagram of a transceiver using HAD is plotted. In the
figure, there are two subfigures (a) and (b). They show the fact that the
working process of the practical DM system consists of two steps. Before
performing the beamforming operation of DM, we should measure the directional
angle of desired user or its statistical property as shown in subfigure (a).
Subsequently, in subfigure (b), using the measured value of DOA or its
statistical property, the beamforming scheme can be designed. In both
subfigures, the sub-connected hybrid architecture is adopted, where each array
is made up of $K$ sub-arrays of antennas and each sub-array has $M$ antennas.
In other words, $N=KM$, where $N$ is the total number of antennas at Alice.
Below, we establish the system models of subfigures (a) and (b).

Fig. 1 (a) is a HAD receiver at Alice. In this subfigure, a narrow-band signal
$s(t)e^{jw_{c}t}$ from a far-field emitter Bob impinges on the HAD receive
array at Alice, the receive signal is

$\displaystyle\mathbf{y}_{a}(t)=\mathbf{F}_{RF}^{H}\mathbf{a}(\theta)s(t)e^{jw_{c}t}+\mathbf{n}(t)$
(1)

where $s(t)$ is the baseband signal, $w_{c}=2\pi f_{c}$ with $f_{c}$ being
carrier frequency, and $\mathbf{a}(\theta)$ is the array manifold defined as

$\displaystyle\mathbf{a}(\theta)=[1,e^{j\frac{2\pi}{\lambda}d\cos\theta},\cdots,e^{j\frac{2\pi}{\lambda}(N-1)d\cos\theta}]^{T}$
(2)

which is used to estimate the angle range from $0^{\circ}$ to $180^{\circ}$
and convenient to design beamformer for DM in the following sections. In (1),
$\mathbf{F}_{RF}$ is the $N\times K$ analog receive beamforming matrix
consisting of $K$ analog weight vectors $\mathbf{f}_{k}$, with each element
having the same amplitude $\frac{1}{\sqrt{M}}$ but different phases, as
follows

$\displaystyle\mathbf{F}_{RF}=\left[{\begin{array}[]{*{20}{c}}{\mathbf{f}_{1}}&{\mathbf{0}}&{\cdots}&{\mathbf{0}}\\\
{\mathbf{0}}&{\mathbf{f}_{2}}&{\cdots}&{\mathbf{0}}\\\
{\vdots}&{\vdots}&{\ddots}&{\vdots}\\\
{\mathbf{0}}&{\mathbf{0}}&{\cdots}&{\mathbf{f}_{K}}\end{array}}\right]$ (7)

where
$\mathbf{f}_{k}=\frac{1}{\sqrt{M}}\left[e^{j\alpha_{k,1}},\cdots,e^{j\alpha_{k,M}}\right]^{T}$.
Then, after performing the down-conversion and analog-to-digital conversion
(ADC) operations, the receive signal at Alice can be expressed as

$\displaystyle\mathbf{y}_{a}(n,\theta)=\mathbf{F}_{RF}^{H}\mathbf{a}(\theta)s(n)+\mathbf{n}$
(8)

Finally, after the digital beamforming operation, the receive signal can be
written as

$\displaystyle
r(n)=\mathbf{F}_{BB}^{H}\mathbf{F}_{RF}^{H}\mathbf{a}(\theta)s(n)+\mathbf{F}_{BB}^{H}\mathbf{n}$
(9)

By making use of the above model, we will develop an improved ESPRIT DOA
estimator for HAD architecture to learn the direction angle of desired user.
The specific DOA estimation process will be presented in Section III.

After obtaining the value of DOA of Bob, now we show how to implement a DM
transmitter using HAD beamforming as indicated in Fig. 1 (b). Here, Alice
works on transmit model. The DM transmitter is equipped with $N$ uniformly
spaced linear antennas array. The transmit signal is

$\displaystyle\mathbf{s}=\sqrt{\beta
P_{s}}\mathbf{V}_{RF}\mathbf{v}_{BB}x+\sqrt{(1-\beta)P_{s}}\mathbf{V}_{RF}\mathbf{T}_{BB}\mathbf{z}$
(10)

where $P_{s}$ is the total transmit power constraint, $\beta$ denotes the
power allocation (PA) factor of confidential messages, and $1-\beta$ stands
for the PA factor of AN. Additionally, $\mathbf{V}_{RF}\in\mathbb{C}^{N\times
K}$ and $\mathbf{v}_{D}\in\mathbb{C}^{K\times 1}$ represent the analog and
digital beamforming martices of carrying confidential messages to desired
user, respectively. $\mathbf{V}_{RF}$ has the following form

$\displaystyle\mathbf{V}_{RF}=\left[{\begin{array}[]{*{20}{c}}{\mathbf{v}_{1}}&{\mathbf{0}}&{\cdots}&{\mathbf{0}}\\\
{\mathbf{0}}&{\mathbf{v}_{2}}&{\cdots}&{\mathbf{0}}\\\
{\vdots}&{\vdots}&{\ddots}&{\vdots}\\\
{\mathbf{0}}&{\mathbf{0}}&{\cdots}&{\mathbf{v}_{K}}\end{array}}\right]$ (15)

where $\mathbf{v}_{k}$ is the analog beamforming vector of subarray $k$,
$\mathbf{T}_{BB}\in\mathbb{C}^{K\times N_{s}}$ denotes the digital AN
projection matrices with constraints
$\parallel\mathbf{V}_{RF}\mathbf{T}_{BB}\parallel^{2}_{F}=\parallel\mathbf{V}_{RF}\mathbf{v}_{BB}\parallel^{2}=1$.
$x\sim\mathcal{CN}(0,1)$ and $\mathbf{z}\sim\mathcal{CN}(0,\mathbf{I})$ are
the confidential message and AN interfering the eavesdropper, respectively.
Let us define $\Omega({\theta})=\frac{2\pi d}{\lambda}\cos({\theta})$ and the
total steering vector (SV) in free space as

$\displaystyle\mathbf{h}^{H}(\theta)=[$ $\displaystyle
e^{j(1-\frac{N+1}{2})\frac{2\pi
d}{\lambda}\cos(\theta)},e^{j(2-\frac{N+1}{2})\frac{2\pi
d}{\lambda}\cos(\theta)},\cdots,$ $\displaystyle
e^{j(N-\frac{N+1}{2})\frac{2\pi
d}{\lambda}\cos(\theta)}]=[\mathbf{h}_{1}^{H}(\theta),\cdots,\mathbf{h}_{K}^{H}(\theta)].$
(16)

where $\mathbf{h}_{k}(\theta)^{H}$ corresponds to the steering vector of
subarray $k$ given by

$\displaystyle\mathbf{h}_{k}(\theta)^{H}=\frac{1}{\sqrt{M}}\left[e^{j\alpha_{k,1}},\cdots,e^{j\alpha_{k,M}}\right]$
(17)

where

$\displaystyle\alpha_{k,m}=\left[(k-1)M+m-\frac{N+1}{2}\right]\Omega(\theta).$
(18)

Interestingly, the total steering vector $\mathbf{h}^{H}(\theta)$ is also
written in the following form

$\displaystyle\mathbf{h}^{H}(\theta)=\mathbf{h}^{H}_{\Theta}(\theta)\otimes\mathbf{h}^{H}_{1}(\theta)$
(19)

where $\mathbf{h}^{H}_{\Theta}(\theta)$ is defined as

$\displaystyle\mathbf{h}^{H}_{\Theta}(\theta)=\left[e^{j0},~{}e^{jM\Omega(\theta)},~{}\cdots,e^{-j(K-1)M\Omega(\theta)}\right].$
(20)

The receive signal at the desired user Bob can be expressed as

$\displaystyle y(\theta_{d})$
$\displaystyle=\mathbf{h}^{H}(\theta_{d})\mathbf{s}+n_{d}$
$\displaystyle=\sqrt{\beta
P_{s}}\mathbf{h}^{H}(\theta_{d})\mathbf{V}_{RF}\mathbf{v}_{BB}x$
$\displaystyle+\sqrt{(1-\beta)P_{s}}\mathbf{h}^{H}(\theta_{d})\mathbf{V}_{RF}\mathbf{T}_{BB}\mathbf{z}+n_{d}$
(21)

where$n_{d}$ is the complex additive white Gaussian noise(AWGN) with
$n_{d}\sim\mathcal{C}\mathcal{N}(0,\sigma_{d}^{2})$. Similarly, the receive
signal at eavesdropper is

$\displaystyle y(\theta_{e})$
$\displaystyle=\mathbf{h}^{H}(\theta_{e})\mathbf{s}+n_{e}$
$\displaystyle=\sqrt{\beta
P_{s}}\mathbf{h}^{H}(\theta_{e})\mathbf{V}_{RF}\mathbf{v}_{BB}x$
$\displaystyle+\sqrt{(1-\beta)P_{s}}\mathbf{h}^{H}(\theta_{e})\mathbf{V}_{RF}\mathbf{T}_{BB}\mathbf{z}+n_{e}$
(22)

where $n_{e}$ is the complex AWGN following
$n_{e}\sim\mathcal{C}\mathcal{N}(0,\sigma_{e}^{2})$. In particular, it is
noted that Eve is assumed to be a passive eavesdropper in our paper. According
to (7), the achievable rate of the desired user Bob is

$\displaystyle R(\theta_{d})=log_{2}(1+\frac{\beta
P_{s}\parallel\mathbf{h}^{H}(\theta_{d})\mathbf{V}_{RF}\mathbf{v}_{BB}\parallel^{2}}{(1-\beta)P_{s}\parallel\mathbf{h}^{H}(\theta_{d})\mathbf{V}_{RF}\mathbf{T}_{BB}\parallel^{2}+\sigma^{2}})$
(23)

Similarly, the achievable rate of eavesdropper is

$\displaystyle R(\theta_{e})=log_{2}(1+\frac{\beta
P_{s}\parallel\mathbf{h}^{H}(\theta_{e})\mathbf{V}_{RF}\mathbf{v}_{BB}\parallel^{2}}{(1-\beta)P_{s}\parallel\mathbf{h}^{H}(\theta_{e})\mathbf{V}_{RF}\mathbf{T}_{BB}\parallel^{2}+\sigma^{2}})$
(24)

Then the SR can be defined as the rate difference between $R(\theta_{b})$ and
$R(\theta_{e})$

$\displaystyle R_{S}=\max\\{0,R(\theta_{d})-R(\theta_{e})\\}.$ (25)

##  III Proposed ML-based framework for DOA estimation using HAD structure

In this section, when classic ESPRIT is applied to HAD receiver, there exists
a serious problem of phase ambiguity. By maximizing the receive power over a
set of finite potential angles, we completely solve this problem to form an
improved ESPRIT. Using this improved ESPRIT and histogram method in ML, we
find an important fact: the DOAME obeys a Gaussian distribution. Then, TDS is
constructed. Via ML framework, the value of DOA and its density is inferred
from this TDS and the following RTS.

###  III-A Improved ESPRIT-based HAD DOA Estimator

First, let us assume the phases of all elements of analog beamforming matrix
are equal to zeros, that is, each analog subarray satisfies

$\displaystyle\mathbf{f}_{k}=\frac{1}{\sqrt{M}}\left[1,\cdots,1\right]^{T}.$
(26)

Then, in accordance with (4), we can directly obtain the output vector after
analog beamforming as follows

$\displaystyle\mathbf{y}_{a}(n,\theta)=\mathbf{F}_{RF}^{H}\mathbf{a}(\theta)s(n)+\mathbf{n}$
$\displaystyle=\frac{1}{\sqrt{M}}\left[{\begin{array}[]{*{20}{c}}{1\cdots
1}&{\mathbf{0}}&{\cdots}&{\mathbf{0}}\\\ {\mathbf{0}}&{1\cdots
1}&{\cdots}&{\mathbf{0}}\\\ {\vdots}&{\vdots}&{\ddots}&{\vdots}\\\
{\mathbf{0}}&{\mathbf{0}}&{\cdots}&{1\cdots
1}\end{array}}\right]\mathbf{a}(\theta)s(n)+\mathbf{n}$ (31)
$\displaystyle=\frac{1}{\sqrt{M}}f(\theta)\mathbf{a}_{K}(\theta)s(n)+\mathbf{n}$
(32)

where

$\displaystyle
f(\theta)=\sum_{m=0}^{M-1}e^{j\frac{2\pi}{\lambda}md\cos\theta}=\frac{1-e^{j2\pi/\lambda
Md\cos\theta}}{1-e^{j2\pi/\lambda d\cos\theta}},$ (33)

and

$\displaystyle\mathbf{a}_{K}(\theta)=\left[1,e^{j\frac{2\pi}{\lambda}Md\cos\theta},\cdots,e^{j\frac{2\pi}{\lambda}(K-1)Md\cos\theta}\right]^{T}.$
(34)

Let
$\mathbf{a}_{D}(\theta)=\frac{1}{\sqrt{M}}f(\theta)\mathbf{a}_{K}(\theta)$,
then $\mathbf{y}_{a}(n)$ can be recasted as

$\displaystyle\mathbf{y}_{a}(n)=\mathbf{a}_{D}(\theta)s(n)+\mathbf{n}.$ (35)

Next, we adopt ESPRIT to estimate the directional angle of emitter. We choose
the first $K-1$ sub-array of antennas to be Array 1, and choose the last $K-1$
sub-array of antennas to be Array 2, then the output can be expressed as

$\displaystyle\mathbf{y}_{a}^{1}(n,~{}\theta)=\mathbf{a}_{D}^{1}(\theta)s(n)+\mathbf{n}^{1},$
(36)

and

$\displaystyle\mathbf{y}_{a}^{2}(n)$
$\displaystyle=\mathbf{a}_{D}^{2}(\theta)s(n)+\mathbf{n}^{2}$
$\displaystyle=\mathbf{a}_{D}^{2}(\theta)\mathbf{\Phi}s(n)+\mathbf{n}^{2},$
(37)

respectively, where

$\displaystyle\mathbf{a}_{D1}(\theta)=\frac{f(\theta)}{\sqrt{M}}[1,e^{j\frac{2\pi}{\lambda}Md\cos\theta},\cdots,e^{j\frac{2\pi}{\lambda}(K-2)Md\cos\theta}]^{T}$
(38)

and $\mathbf{\Phi}=\text{diag}\\{e^{j\frac{2\pi}{\lambda}Md\cos\theta}\\}$. By
computing the co-variance matrix of $\mathbf{y}_{a}^{1}(n)$ and
$\mathbf{y}_{a}^{2}(n)$, separately, we can obtain the signal spaces
$\mathbf{E}_{1}$ and $\mathbf{E}_{2}$. Because of the translational
displacement, there must exit a non-singular transformation matrix
$\mathbf{\Psi}$ such that

$\displaystyle\mathbf{E}_{1}\mathbf{\Psi}=\mathbf{E}_{2},$ (39)

and there must exit a non-singular transformation matrix $\mathbf{T}$ such
that

$\displaystyle\mathbf{E}_{1}=\mathbf{a}_{D1}\mathbf{T},$ (40)

and

$\displaystyle\mathbf{E}_{2}=\mathbf{a}_{D1}\mathbf{\Phi}\mathbf{T}.$ (41)

According to Eqs. (13), (14) and (15), we can derive a crucial relationship

$\displaystyle\mathbf{T}\mathbf{\Psi}\mathbf{T}^{-1}=\mathbf{\Phi}.$ (42)

Therefore, it is evident that the eigenvalues of $\mathbf{\Psi}$ are equal to
the diagonal elements of $\mathbf{\Phi}$, i.e.

$\displaystyle\lambda_{i}=\exp\left\\{j\frac{2\pi}{\lambda}Md\cos\theta\right\\}.$
(43)

Now the problem is how to calculate the rotation operator $\mathbf{\Psi}$.
Next, we use the total least-squares in [13] method to estimate
$\mathbf{\Psi}$. The correlation matrix of the entire antenna array is

$\displaystyle\mathbf{R}_{yy}$
$\displaystyle=E[\mathbf{y}_{a}(n)\mathbf{y}_{a}^{H}(n)]\approx\frac{1}{L}\\{\mathbf{Y}_{a}\mathbf{Y}_{a}^{H}\\}$
$\displaystyle\approx\mathbf{a}_{D}\sigma_{s}^{2}P_{b}^{2}\mathbf{a}_{D}^{H}+\sigma_{n}^{2}\mathbf{I}_{K}$
(44)

where
$\mathbf{Y}_{a}=[\mathbf{y}_{a}(1),\mathbf{y}_{a}(2),\cdots,\mathbf{y}_{a}(L)]$
is the matrix of spatial-time sampling points. The correlation matrix of the
entire antenna array has $K$ eigenvalues, and the associated $K$ eigenvectors
$\mathbf{E}=[\mathbf{e}_{1}\mathbf{e}_{2}\cdots\mathbf{e}_{K}]$. Signal
subspace $\mathbf{E}_{s}$ consists of the eigenvectors corresponding to the
largest eigenvalues of the matrix $\mathbf{R}_{yy}$. Since signal spaces
$\mathbf{E}_{1}$ and $\mathbf{E}_{2}$ can be constructed from the entire
antennas array signal space $\mathbf{E}_{s}$, we choose the first $K-1$ rows
of $\mathbf{E}_{s}$ to construct $\mathbf{E}_{1}$, and similarly, the last
$K-1$ rows to compose $\mathbf{E}_{2}$. Next, we utilize signal subspaces to
construct a $2\times 2$ matrix.

$\displaystyle\mathbf{C}=\left[{\begin{array}[]{*{20}{c}}{\mathbf{E}_{1}^{H}}\\\
{\mathbf{E}_{2}^{H}}\\\
\end{array}}\right]\left[\mathbf{E}_{1}\mathbf{E}_{2}\right]=\mathbf{E}_{C}\Lambda\mathbf{E}_{C}^{H}$
(47)

where $\mathbf{E}_{c}$ can be derived from the EVD (Eigen value decomposition)
of matrix $\mathbf{C}$ and $\Lambda=\text{diag}\\{\lambda_{1},\lambda_{2}\\}$.
Furthermore, $\mathbf{E}_{c}$ can be decomposed into four subarrays as follows

$\displaystyle\mathbf{E}_{c}=\left[{\begin{array}[]{*{20}{c}}{\mathbf{E}_{11}}&{\mathbf{E}_{12}}\\\
{\mathbf{E}_{21}}&{\mathbf{E}_{22}}\\\ \end{array}}\right]$ (50)

Then, we can obtain the rotation operator $\mathbf{\Psi}$

$\displaystyle\mathbf{\Psi}=-\mathbf{E}_{12}\mathbf{E}_{22}^{-1}$ (51)

By calculating the eigenvalue $\lambda_{1}$ of $\mathbf{\Psi}$, the angle of
emitter estimated by traditional ESPRIT method is

$\displaystyle\hat{\theta}=\arccos(\frac{\lambda\arg(\lambda_{1})}{2\pi Md})$
(52)

However, the periodicity in equation (43) generates an effect of phase
blurring, which makes estimated angle no longer the expression in
equation(52), thus, we need to extend the estimated angle to a set of
potential solutions

$\displaystyle\hat{\Theta}=\\{\hat{\theta_{i}},i\in{0,1,\cdots,M-1}\\}$ (53)

where

$\displaystyle\hat{\theta_{i}}=\arccos(\frac{\lambda(arg(\lambda_{1})+2\pi
i)}{2\pi Md})$ (54)

Then, computing the receive signals power corresponding to the $M$ angles
individually, we choose the angle of maximizing the output power as the
optimal value casted as

$\hat{\theta}=\mathop{\arg\min}_{\hat{\theta_{i}}\in\hat{\Theta}}\ \
P_{r}(\hat{\theta_{i}})=\frac{1}{LN^{2}}\sum_{n=1}^{L}(r(n)r(n)^{H})$ (55)

where

$\displaystyle r(n)$ $\displaystyle=\mathbf{F}_{BB}^{H}y_{a}(n)$
$\displaystyle=\sum_{k=1}^{K}e^{-j\alpha_{k}}y_{a}(n)$
$\displaystyle=\frac{1}{\sqrt{M}}s(n)\sum_{k=1}^{K}e^{j\frac{2\pi}{\lambda}(k-1)Md(\cos\theta_{0}-\cos\hat{\theta}_{i})}$
$\displaystyle\times\sum_{m=1}^{M}e^{j\frac{2\pi}{\lambda}(m-1)d(\cos\theta_{0}-\cos\hat{\theta}_{i})}+\sum_{k=1}^{K}e^{-j\frac{2\pi}{\lambda}(k-1)Md\cos\hat{\theta_{i}}}n$
(56)

Since we have obtained the value of DOA of incoming signal, we can further
estimate the SNR of incoming signal. Given the estimated angle $\hat{\theta}$,
the analog beamforming vector (ABV) of each subarray can be designed as

$\displaystyle\hat{\mathbf{f}}_{k}=\frac{1}{\sqrt{M}}[e^{j\frac{2\pi}{\lambda}((k-1)M+1)d\cos\hat{\theta}},\cdots,e^{j\frac{2\pi}{\lambda}((k-1)M+M)d\cos\hat{\theta}}]^{T}$
(57)

Then, after the ABV $\hat{\mathbf{f}}_{k}$ is performed in RF chain, the
output signal of subarray $k$ in the next time slot can be given by

$\displaystyle\hat{y}_{a,k}(n)=\hat{\mathbf{f}}_{k}^{H}\mathbf{a}_{k}(\theta_{0})s(n)+n.$
(58)

If the estimated angle $\hat{\theta}$ is equal to the exact angle
$\theta_{0}$, then $\hat{\mathbf{f}}_{k}^{H}\mathbf{a}_{k}(\theta_{0})=1$, and
(58) can be recasted as

$\displaystyle\hat{y}_{a,k}(n)=s(n)+n$ (59)

which yields the output vector of all subarrays

$\displaystyle\hat{\mathbf{y}}_{a}(n)=\left[\hat{y}_{a,1}(n),\hat{y}_{a,2}(n),\cdots,\hat{y}_{a,K}(n)\right]^{T}.$
(60)

Then we have the overall estimated SNR $\hat{\rho}$ in HAD architecture as
follows

$\displaystyle\hat{\rho}=\frac{\hat{\mathbf{y}}_{a}(n)^{T}\hat{\mathbf{y}}_{a}(n)^{*}-K\sigma_{n}^{2}}{K\sigma_{n}^{2}}.$
(61)

where $\sigma_{n}^{2}$ is taken to be the least eigen-value of matrix (III-A).
Until now, we complete the DOA and SNR estimate by the improved ESPRIT, called
improved HAD-ESPRIT (I-HAD-ESPRIT) in what follows.

###  III-B Proposed ML-based DOA and its density estimation

Due to the effect of channel noise or co-channel signals, there exists DOA
measurement error in the estimated DOA in the preceding subsection.
Considering DOAME, the estimated angle can be modeled as

$\displaystyle\tilde{\theta}(n)=\theta_{0}+\Delta\theta(n),~{}~{}n=1,\cdots,N_{DOA}$
(62)

where $\tilde{\theta}(n)$ represents the practical observation angle, $\theta$
is the ideal angle we want to obtain, and $\Delta\theta(n)$ is the measurement
error between the ideal angle and measured one, $n=1,2,\cdots,~{}N$ denotes
the number of time slots. Here, the proposed I-HAD-ESPRIT is used one time to
measure DOA per time slot.

Figure 2: DOA measurement error
Figure 3: Flow diagram of ML-based DOA estimator

Using the proposed I-HAD-ESPRIT in the preceding subsection repeatedly for
several hundreds or thousands, Fig. 2 depicts the histogram of DOAME
$\Delta\theta(n)$ provided that $N=32$, $M=4$, and SNR=0dB. From this figure,
we can find the important fact that the measured value of DOA approximately
obeys a Gaussian distribution with zero mean.

$\displaystyle
f\left(\tilde{\theta}(n);\theta,\sigma^{2}\right)=\frac{1}{\sqrt{2\pi\sigma^{2}}}\exp\\{-\frac{1}{2\sigma^{2}}(\tilde{\theta}(n)-\theta)^{2}\\}$
(63)

If we can learn its variance, then we can give its density function. This is
just the density estimation in ML field. Below, we will show how to infer the
density of DOAME and improve the accuracy of DOA measurement by making use of
ML-based framework.

In practice, a rapid real-time or on-line DOA measurement is required in order
to quickly locate the position of emitter, especially in military field. Thus,
it is not allowable to take too many time slots. This means that the offline
training is very important to infer the density of DOAME and improve the
accuracy of DOA measurement. First, the TDS of DOA measurements is offline
constructed by using the proposed I-HAD-ESPRIT in the previous subsection. Let
us define the TDS as
$S_{TDS}=\left\\{\tilde{\theta}_{1},\tilde{\theta}_{2},\cdots,\tilde{\theta}_{N_{TDS}}\right\\}$,
where $N_{NDS}$ stands for the number of elements in TDS.

Considering the DOA measurement is approximated as a Gaussian distribution and
in terms of the learning rule of maximum likelihood in [41], we have its mean

$\displaystyle\hat{\theta}_{TDS}=\frac{1}{N_{TDS}}\sum_{\tilde{\theta}_{i}\in
S_{TDS}}\tilde{\theta}_{i},$ (64)

and sampling variance

$\displaystyle\hat{\sigma}^{2}_{TDS}=\frac{1}{N_{TDS}-1}\sum_{i=1}^{N_{TDS}}(\tilde{\theta}_{i}-\hat{\theta}_{TDS})^{2}$
(65)

in the training stage. Similarly, in the real-time measuring stage, we have

$\displaystyle\hat{\theta}_{RTS}=\frac{1}{N_{RTS}}\sum_{\tilde{\theta}_{i}\in
S_{RTS}}\tilde{\theta}_{i},$ (66)

and real-time variance

$\displaystyle\hat{\sigma}^{2}_{RTS}=\frac{1}{N_{RTS}-1}\sum_{i=1}^{N_{RTS}}(\tilde{\theta}_{i}-{\hat{\theta}}_{RTS})^{2}.$
(67)

If Bob keeps fixed or moves away from Alice along the fixed desired direction
during both offline and online periods, then we call this scenario a stable
state. Clearly, the above two estimated values mean and variance can be
exploited to improve the accuracy of DOA measurement mean and variance.
Conversely, if Bob moves fast away from the desired direction, we call this
situation as a moving state. In such a situation, at least the estimated
variance during the training stage can be used to improve the accuracy of DOA
measurement variance by predicting the variance of real-time measuring stage.
In what follows, we discuss the two situations. In the stable state, Alice
doesn’t need to adjust its beamforming direction. However, in the moving
state, Alice is required to adjust its beam direction to follow the real-time
varying angle of the desired user.

Due to the two important facts: the SNR of training stage may be different
from the SNR of real-time measuring stage, and the number $N_{TDS}$ of
elements in TDS is also possibly different from that $N_{RTS}$ in RTS, we
should take the two factors into account when we combine the learning output
of the two sets. In general, $N_{TDS}>>N_{RTS}$. In the stable state, the
weight output of means and variances due to TDS and RTS are given by

$\displaystyle\hat{\theta}_{ML}=\alpha_{1}{\hat{\theta}}_{TDS}+\alpha_{2}\hat{\theta}_{RTS}=\theta+\delta\theta_{ML}$
(68)

with the constraint $\alpha_{1}+\alpha_{2}=1$ with $\alpha_{1}\geq 0$ and
$\alpha_{2}\geq 0$. In the above equation, $\alpha_{1}$ and $\alpha_{2}$
depends mainly on the following factors:
$\hat{\rho}_{TDS},\hat{\rho}_{RTS},N_{TDS}$, and $N_{RTS}$, where
$\hat{\rho}_{TDS}$, and $\hat{\rho}_{RTS}$ are the receive SNRs at Alice in
the training and real-time measuring stages, respectively. Based on this, we
propose three convex weight combinators with weight factors as follows

$\displaystyle\alpha_{1}=\frac{\hat{\rho}_{TDS}}{\hat{\rho}_{TDS}+\hat{\rho}_{RTS}},$
$\displaystyle\alpha_{2}=\frac{\hat{\rho}_{RTS}}{\hat{\rho}_{TDS}+\hat{\rho}_{RTS}},$
(69) $\displaystyle\alpha_{1}=\frac{N_{TDS}}{N_{TDS}+N_{RTS}},$
$\displaystyle\alpha_{2}=\frac{N_{RTS}}{N_{TDS}+N_{RTS}},$ (70)

and

$\displaystyle\alpha_{1}=\frac{\hat{\rho}_{TDS}\cdot
N_{TDS}}{\hat{\rho}_{TDS}\cdot N_{TDS}+\hat{\rho}_{RTS}\cdot N_{RTS}},$
$\displaystyle\alpha_{2}=\frac{\hat{\rho}_{RTS}\cdot
N_{RTS}}{\hat{\rho}_{TDS}\cdot N_{TDS}+\hat{\rho}_{RTS}\cdot N_{RTS}},$ (71)

respectively. Here, for the sake of convenience, in what follows, (III-B),
(III-B) and (LABEL:Both-Weighting) are called the proposed weight methods i@,
ii@, and iii@, respectively. Similarly, we have the corresponding weight
combiners for DOA variance estimation

$\displaystyle\hat{\sigma^{2}}_{ML}$
$\displaystyle=\alpha_{1}\tilde{\sigma}^{2}_{TDS}+\alpha_{2}\hat{\sigma}^{2}_{NDS}$
(72)

where the weight parameters $\alpha_{1}$ and $\alpha_{2}$ are given by
(III-B), (III-B), and (LABEL:Both-Weighting), respectively. And

$\tilde{\sigma}^{2}_{TDS}=\frac{\hat{\rho}_{TDS}}{\hat{\rho}_{RTS}}\hat{\sigma}^{2}_{TDS}$
(73)

In the moving scenario, it is particularly noted that the learned variance in
training stage has an intimately direct relationship with real-time measuring
one. At least it is used to predict the variance of the real-time measuring
case in terms of their SNRs, and numbers of RTS and TDS. However, due to
moving, the direction angles of Bob and Eve will change with time, the learned
DOA mean in RTS will be independent of the learned DOA mean from TDS. In other
words, the learned DOA mean from RTS will be used as the final learned DOA. In
such a situation, the weight factors $\alpha_{1}$ and $\alpha_{2}$ in (68) are
set to 1 and 0, respectively.

##  IV Proposed robust HAD beamformer for DM

In the above section, we have attained the measured value $\hat{\theta}_{ML}$
and variance $\hat{\sigma^{2}}_{ML}$ of DOA, then we have the ML-based DOA
output modeled as

$\displaystyle\theta=\hat{\theta}_{ML}+\Delta\hat{\theta}$ (74)

with its PDF as follows

$\displaystyle
f\left(\theta;\hat{\theta}_{ML},\hat{\sigma}^{2}_{ML}\right)=\frac{1}{\sqrt{2\pi\hat{\sigma}^{2}_{ML}}}\exp\left\\{-\frac{1}{2\hat{\sigma}^{2}_{ML}}\left(\theta-\hat{\theta}_{ML}\right)^{2}\right\\},$
(75)

which yields the PDF of DOAME

$\displaystyle
f\left(\Delta\hat{\theta};\hat{\sigma}^{2}_{ML}\right)=\frac{1}{\sqrt{2\pi\hat{\sigma}^{2}_{ML}}}\exp\left\\{-\frac{1}{2\hat{\sigma}^{2}_{ML}}\left(\Delta\hat{\theta}\right)^{2}\right\\}.$
(76)

Considering an important fact that the DOAME falls into the interval
$[\pi,\pi]$ , similar to [26, 29], the PDF of DOAME can be approximated by the
truncated Gaussian density function

$\displaystyle
f(\Delta\hat{\theta})=\begin{cases}\frac{1}{K_{d}\sqrt{2\pi\hat{\sigma}^{2}_{ML}}}\exp\left(-\frac{(\Delta\hat{\theta})^{2}}{2\hat{\sigma}^{2}_{ML}}\right),&{\frac{\Delta\hat{\theta}}{\Delta\theta_{\max}}\in[-1,1]}\\\
0,\ &\mbox{others}\end{cases}$ (77)

where the truncated factor $K_{d}$ is given by

$\displaystyle
K_{d}=\int_{-\Delta\theta_{\max}}^{\Delta\theta_{\max}}\frac{1}{\sqrt{2\pi\hat{\sigma}^{2}_{ML}}}\exp\left(-\frac{({\Delta\hat{\theta}})^{2}}{2\hat{\sigma}^{2}_{ML}}\right)d(\Delta\hat{\theta}),$
(78)

and $\Delta\theta_{\max}$ may be chosen to be a value less than or equal to
$\pi$ but larger than a multiple of $\sigma$. For example,
$\Delta\theta_{\max}$ can be taken to be half of the main-beam width.

Figure 4: Block diagram of the proposed robust DM synthesis method

In accordance with the above approximation and in terms of conditional
expectation, three robust beammforming matrices $\mathbf{V}_{RF}$,
$\mathbf{v}_{BB}$, and $\mathbf{T}_{BB}$, are in order presented for HAD-based
DM systems in the following.

###  IV-A Proposed robust Analog beamforming matrix $\mathbf{V}_{RF}$

Since we have already obtained the direction of desired user in Section III, a
direct method to design the analog precoder is to directly implement phase
alignment (PA), which aligns the phases of phase shifters of subarray $k$ in
analog part to the direction of desired user. The block submatrix $k$
$\mathbf{v}_{k}$ in $\mathbf{V}_{RF}$ is given by

$\displaystyle\mathbf{v}_{NRAB,k}\left(\hat{\theta}_{ML}\right)=\mathbf{h}_{k}\left(\hat{\theta}_{ML}\right)$
(79)

where the SV $\mathbf{h}_{k}(\bullet)$ has been defined in (17). The PA method
to design the analog precoder in (46) doesn’t take the DOA estimation error
into consideration, which completely depends on the precision of DOA
estimation. However, there always exists measurement error in the measured
angle as shown in (68). To combat this type of errors, we will exploit the
statistical property of DOAME, i.e. its PDF to propose a robust analog
beamformer (RAB). In this precoder, given the PDF of DOAME and measured DOA,
minimizing the mean square error criterion directly yields the following
robust beamforming vector for subarray $k$

$\displaystyle\mathbf{v}_{RAB,k}$
$\displaystyle=\text{E}\left\\{\mathbf{h}_{k}\left(\theta\right)|\hat{\theta}_{ML},f\left(\Delta\theta_{ML}\right)\right\\}$
$\displaystyle=\text{E}_{\Delta\theta_{ML}}\left\\{\mathbf{h}_{k}\left(\hat{\theta}_{ML}+\Delta\theta_{ML}\right)\right\\}$
(80)

where the $m$ th element of $\mathbf{v}_{RAB,k}$ is derived in Appendix A.
Placing all $\mathbf{v}_{RAB,k}$ into (15), we obtain the robust analog
beamforming matrix $\mathbf{{V}}_{RF,RAB}$.

###  IV-B Proposed robust digital beamforming vector $\mathbf{v}_{BB}$

Now, we compute the robust digital beamforming (RDB) vector $\mathbf{v}_{BB}$.
Similar to [37, 40], the problem of solving $\mathbf{v}_{BB}$ can be casted as

$\displaystyle\min\limits_{\mathbf{v}_{D}}{\kern
1.0pt}~{}~{}~{}~{}~{}~{}\parallel\mathbf{v}_{FD}-\mathbf{V}_{RF,RAB}\mathbf{v}_{BB}\parallel_{2}^{2}$
$\displaystyle\text{subject~{}to}~{}~{}~{}~{}~{}\parallel\mathbf{V}_{RF,RAB}\mathbf{v}_{BB}\parallel^{2}=1.$
(81)

where $\mathbf{v}_{FD}$ stands for the optimal FD beamforming vector of
confidential message. Observing the above the above optimization problem, we
find an obvious fact: if we design $\mathbf{v}_{FD}$ well, then we easily
address the problem of optimzing $\mathbf{v}_{BB}$ with $\mathbf{V}_{RF,RAB}$
available in the precedeing subsection. Here, $\mathbf{v}_{FD}$ is chosen to
be the beamforming vector of confidential message of robust null-space
projection method in [42]. Due to the special structure of
$\mathbf{V}_{RF,RAB}$, we have the identity
$\mathbf{V}^{H}_{RF,RAB}\mathbf{V}_{RF,RAB}=\mathbf{I}$, and further
$\parallel\mathbf{V}_{RF,RAB}\mathbf{v}_{BB}\parallel^{2}=\parallel\mathbf{v}_{BB}\parallel^{2}=1$.
So the optimization problem in (IV-B) can be simplified as

$\displaystyle\min\limits_{\mathbf{v}_{BB}}{\kern
1.0pt}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\parallel\mathbf{v}_{FD}-\mathbf{V}_{RF,RAB}\mathbf{v}_{BB}\parallel^{2}_{2}$
$\displaystyle\text{subject~{}to}~{}~{}~{}~{}~{}\parallel\mathbf{v}_{BB}\parallel^{2}=1$
(82)

whose objective function can be expanded as

$\displaystyle
g\left(\mathbf{v}_{BB}\right)=\parallel\mathbf{v}_{FD}-\mathbf{V}_{RF,RAB}\mathbf{v}_{BB}\parallel^{2}_{2}$
$\displaystyle=(\mathbf{v}_{FD}-\mathbf{V}_{RF,RAB}\mathbf{v}_{BB})^{H}(\mathbf{v}_{FD}-\mathbf{V}_{RF,RAB}\mathbf{v}_{BB})$
$\displaystyle=\mathbf{v}_{FD}^{H}\mathbf{v}_{FD}-\mathbf{v}_{FD}^{H}\mathbf{V}_{RF,RAB}\mathbf{v}_{BB}-\mathbf{v}_{BB}^{H}\mathbf{V}_{RF,RAB}^{H}\mathbf{v}_{FD}$
$\displaystyle~{}~{}~{}~{}+\mathbf{v}_{BB}^{H}\mathbf{V}_{RF,RAB}^{H}\mathbf{V}_{RF,RAB}\mathbf{v}_{BB}$
$\displaystyle=\mathbf{v}_{FD}^{H}\mathbf{v}_{FD}-\mathbf{v}_{FD}^{H}\mathbf{V}_{RF,RAB}\mathbf{v}_{BB}$
$\displaystyle-\mathbf{v}_{BB}^{H}\mathbf{V}_{RF,RAB}^{H}\mathbf{v}_{FD}+\mathbf{v}_{BB}^{H}\mathbf{v}_{BB}.$
(83)

Differentiating the objective function $g\left(\mathbf{v}_{BB}\right)$ with
respect to $\mathbf{v}_{BB}$ and setting it equal to zero yields

$\displaystyle\frac{\partial
g\left(\mathbf{v}_{BB}\right)}{\partial\mathbf{v}_{BB}^{H}}=\mathbf{v}_{BB}-\mathbf{V}_{RF,RAB}^{H}\mathbf{v}_{FD}=0,$
(84)

which gives

$\displaystyle\mathbf{v}_{BB}=\mathbf{V}_{RF,RAB}^{H}\mathbf{v}_{FD},$ (85)

which is normalized to a unit vector

$\displaystyle\mathbf{v}_{RBB,~{}BB}=\frac{\mathbf{V}_{RF,RAB}^{H}\mathbf{v}_{FD}}{\parallel\mathbf{V}_{RF,RAB}^{H}\mathbf{v}_{FD}\parallel_{2}}.$
(86)

with the constraint $\mathbf{v}_{BB}^{H}\mathbf{v}_{BB}$=1 in (IV-B). This
completes the construction of $\mathbf{v}_{BB}$.

###  IV-C Proposed robust AN projection matrix $\mathbf{T}_{BB}$

In this subsection, we will optimize the design of the digital AN projection
matrix $\mathbf{T}_{BB}$. Actually, its designing is very similar that of
$\mathbf{v}_{BB}$ in the previous subsection. In the same fashion, after the
FD AN projection matrix $\mathbf{T}_{BB}$ is set to the robust NSP matrix in
[42], the optimization problem of finding $\mathbf{T}_{BB}$ can be casted as

$\displaystyle\min\limits_{\mathbf{T}_{BB}}{\kern
1.0pt}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\parallel\mathbf{T}_{FD}-\mathbf{V}_{RF,RAB}\mathbf{T}_{BB}\parallel^{2}_{F}$
$\displaystyle\text{subject~{}to}~{}~{}~{}~{}~{}\parallel\mathbf{V}_{RF,RAB}\mathbf{T}_{BB}\parallel^{2}_{F}=1$

Because of
$\parallel\mathbf{V}_{RF,RAB}\mathbf{T}_{BB}\parallel_{F}^{2}=\parallel\mathbf{T}_{BB}\parallel_{F}^{2}=1$,
we can further simplify the above optimization problem as

$\displaystyle\min\limits_{\mathbf{T}_{BB}}{\kern
1.0pt}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\parallel\mathbf{T}_{FD}-\mathbf{V}_{RF,RAB}\mathbf{T}_{BB}\parallel^{2}_{F}$
$\displaystyle\text{subject~{}to}~{}~{}~{}~{}~{}\parallel\mathbf{T}_{BB}\parallel^{2}_{F}=1$

The objective function
$\parallel\mathbf{T}_{opt}-\mathbf{V}_{RF,RAB}\mathbf{T}_{BB}\parallel^{2}_{F}$
can be represented as

$\displaystyle\parallel\mathbf{T}_{FD}-\mathbf{V}_{RF,RAB}\mathbf{T}_{BB}\parallel^{2}_{F}$
$\displaystyle=\parallel\text{v}ec(\mathbf{T}_{FD}-\mathbf{V}_{RF,RAB}\mathbf{T}_{BB})\parallel_{2}^{2}$
$\displaystyle=\parallel\text{v}ec(\mathbf{T}_{FD})-\text{v}ec(\mathbf{V}_{RF,RAB}\mathbf{T}_{BB})\parallel_{2}^{2}$
$\displaystyle=\parallel\text{v}ec(\mathbf{T}_{FD})-(\mathbf{I}_{K}\otimes\mathbf{V}_{RF,RAB})\text{v}ec(\mathbf{T}_{BB})\parallel_{2}^{2}$
$\displaystyle=\parallel\mathbf{g}-\mathbf{Q}\mathbf{t}\parallel_{2}^{2}$ (89)

where $\mathbf{g}=\text{v}ec(\mathbf{T}_{FD})$,
$\mathbf{Q}=\mathbf{I}_{K}\otimes\mathbf{V}_{RF,RAB}$ and
$\mathbf{t}=\text{vec}(\mathbf{T}_{BB})$. And also the constraint is
reformulated as $\parallel\mathbf{t}\parallel_{2}^{2}=1$. Hence, the original
optimization problem in (42) reduces to

$\displaystyle\min\limits_{\mathbf{T}_{D}}{\kern
1.0pt}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\parallel\mathbf{g}-\mathbf{Q}\mathbf{t}\parallel^{2}_{2}$
$\displaystyle\text{subject~{}to}~{}~{}~{}~{}~{}\parallel\mathbf{t}\parallel^{2}_{2}=1,$
(90)

whose the objective function can be expanded as

$\displaystyle\parallel\mathbf{g}-\mathbf{Q}\mathbf{t}\parallel^{2}_{2}$
$\displaystyle=(\mathbf{g}-\mathbf{Q}\mathbf{t})^{H}(\mathbf{g}-\mathbf{Q}\mathbf{t})$
$\displaystyle=\mathbf{g}^{H}\mathbf{g}-\mathbf{g}^{H}\mathbf{Q}\mathbf{t}-\mathbf{t}^{H}\mathbf{Q}^{H}\mathbf{g}+\mathbf{t}^{H}\mathbf{Q}^{H}\mathbf{Q}\mathbf{t}$
$\displaystyle\mathop{=}\limits^{(a)}\mathbf{g}^{H}\mathbf{g}-\mathbf{g}^{H}\mathbf{Q}\mathbf{t}-\mathbf{t}^{H}\mathbf{Q}^{H}\mathbf{g}+\mathbf{t}^{H}\mathbf{Q}^{H}\mathbf{Q}\mathbf{t}.$
(91)

Because

$\displaystyle\mathbf{t}^{H}\mathbf{Q}^{H}\mathbf{Q}\mathbf{t}$
$\displaystyle=\mathbf{t}^{H}(\mathbf{I}_{K}\otimes\mathbf{V}_{RF})^{H}(\mathbf{I}_{K}\otimes\mathbf{V}_{RF})\mathbf{t}$
$\displaystyle=\mathbf{t}^{H}(\mathbf{I}_{K}^{H}\otimes\mathbf{V}_{RF}^{H})(\mathbf{I}_{K}\otimes\mathbf{V}_{RF})\mathbf{t}$
$\displaystyle=\mathbf{t}^{H}(\mathbf{I}_{K}^{H}\mathbf{I}_{K}\otimes\mathbf{V}_{RF}^{H}\mathbf{V}_{RF})\mathbf{t}$
$\displaystyle=\mathbf{t}^{H}(\mathbf{I}_{K}\otimes\mathbf{I}_{K})\mathbf{t}$
$\displaystyle=\mathbf{t}^{H}\mathbf{t},$ (92)

the objective function reduces to

$\displaystyle\parallel\mathbf{g}-\mathbf{Q}\mathbf{t}\parallel^{2}_{2}=\mathbf{g}^{H}\mathbf{g}-\mathbf{g}^{H}\mathbf{Q}\mathbf{t}-\mathbf{t}^{H}\mathbf{Q}^{H}\mathbf{g}+\mathbf{t}^{H}\mathbf{t}.$
(93)

Taking the derivative of
$\parallel\mathbf{g}-\mathbf{Q}\mathbf{t}\parallel^{2}_{2}$ with respect to
$\mathbf{t}^{H}$ and setting it to zero, we have

$\displaystyle\frac{\partial\parallel\mathbf{g}-\mathbf{Q}\mathbf{t}\parallel^{2}_{2}}{\partial\mathbf{t}^{H}}=\mathbf{t}-\mathbf{Q}^{H}\mathbf{g}=0,$
(94)

Similarly, the optimal solution is

$\displaystyle\mathbf{t}_{opt}=\frac{\mathbf{Q}^{H}\mathbf{g}}{\parallel\mathbf{Q}^{H}\mathbf{g}\parallel}.$
(95)

where $\mathbf{Q}^{H}\mathbf{g}$ can be rewritten in the following form

$\displaystyle\mathbf{Q}^{H}\mathbf{g}=$
$\displaystyle\left[{\begin{array}[]{*{20}{c}}{\mathbf{V}_{RF,RAB}^{H}}&{\mathbf{0}}&{\cdots}&{\mathbf{0}}\\\
{\mathbf{0}}&{\mathbf{V}_{RF,RAB}^{H}}&{\cdots}&{\mathbf{0}}\\\
{\vdots}&{\vdots}&{\ddots}&{\vdots}\\\
{\mathbf{0}}&{\mathbf{0}}&{\cdots}&{\mathbf{V}_{RF,RAB}^{H}}\end{array}}\right]\left[{\begin{array}[]{*{20}{c}}{\mathbf{g}_{1}}\\\
{\mathbf{g}_{2}}\\\ {\vdots}\\\ {\mathbf{g}_{N}}\end{array}}\right]$ (104)
$\displaystyle=$
$\displaystyle\left[{\begin{array}[]{*{20}{c}}{\mathbf{V}_{RF,RAB}^{H}\mathbf{g}_{1}}\\\
{\mathbf{V}_{RF,RAB}^{H}\mathbf{g}_{2}}\\\ {\vdots}\\\
{\mathbf{V}_{RF,RAB}^{H}\mathbf{g}_{N}}\end{array}}\right].$ (109)

Using (95) and $\mathbf{t}=\text{vec}(\mathbf{T}_{BB})$, the robust digital AN
projection matrix $\mathbf{T}_{opt}$ is readily gotten by just the reverse
operation of matrix-to-vector operator $\text{vec}(\bullet)$.

##  V Simulation Results and discussions

In this section, simulation results are presented to evaluate the performance
of the proposed I-HAD-ESPRIT DOA estimator,its enhanced version ML-based
scheme, and robust HAD-based beamforming schemes for DM. Simulation parameters
are chosen as follows: $\theta_{d}=50^{\circ}$ and $\theta_{e}=70^{\circ}$.
The number of snapshots per DOA measurement is set to be $L=64$, and the total
number of antennas at Alice is $N=64$ with $M\in\\{4,~{}8,~{}16$}. The PA
factor $\beta=0.9$ and the antenna spacing is $d=\frac{\lambda}{2}$.
$\sigma_{d}^{2}=\sigma_{e}^{2}=\sigma^{2}$ and QPSK modulation is used to
assess bit error rate (BER) performance.

Figure 5: RMSE versus SNR of the proposed ESPRIT with
$N=32withM=\\{4,8,16\\}$.
Figure 6: RMSE versus SNR of the proposed SNR estimator.
Figure 7: RMSE versus of SNR of the proposed three weight combiners of DOA.
Figure 8: RMSE versus SNR of three weight methods of estimating variance
($\mathrm{SNR_{TDS}=10dB}$). Figure 9: Curves of SR versus SNR of the proposed
robust synthesis method for HAD-based DM with $SNR_{TDS}=-5$dB. Figure 10:
Curves of BER versus SNR of the proposed robust synthesis method for HAD-based
DM with $SNR_{TDS}=-5$dB.

Fig. 5 demonstrates the performance of root mean squared error(RMSE) versus
SNR of the proposed I-HAD-ESPRIT method with $N=64,L=64$ and
$M\in\\{2,8,16\\}$. As can be seen from Fig. 5, we find an interesting
phenomenon that the performance gap between I-HAD-ESPRIT and hybrid CRLB
becomes smaller and smaller as the value of $M$ increases from 2 to 16.
Eventally, this gap tends to zero. And when $M$ is up to 16 (large enough),
the proposed I-HAD-ESPRIT can achieve has the same performance as the hybrid
Root-MUSIC in [18] and achieves hybrid CRLB in the medium and high SNR
regions. However, the proposed I-HAD-ESPRIT performs slightly worse than
hybrid Root-MUSIC in [18] for most cases. The main reason is as follows: when
$M$ increases, the number of RF chains will accordingly reduce such that a
smaller total number of spatial-time sampling points are harvested to
calculate covariance matrix of receive signal.

Fig. 6 plots the RMSE performance of the estimated SNR versus exact SNR.
Similar to Fig. 5, the performance curves of SNR measurement shows the same
trend as that in Fig. 5. As the number of subarray antennas increases, the
precision of SNR estimation decreases gradually. In other words, reducing the
number of RF chains will result in a RMSE performance loss over the estimated
SNR.

Fig. 7 demonstrates the curves of variances of the proposed three different
weight DOA combiners versus SNR for SNRTDS=5dB. The learning output of RTS is
used as a performance upper bound. Due to the fact that the RTS has fewer
number of values of DOA measurement compared to TDS, it performs worst. It can
be seen from Fig. 7: the three proposed weight methods performs better than
the learning output of RTS without use of TDS in the low SNR region because
they all fully exploit the DTS to improve the estimate precision. The proposed
method III is the best one among the three proposed methods due to the fact it
not only utilizes the element numbers of TDS and RTS but also exploits the SNR
difference of TDS and RTS.

Fig. 8 illustrates the curves of variances of the proposed three weight
estimators of angle variance versus SNR with SNRTDS=10dB. The learning output
of RTS is adopted as a performance benchmark. It can be seen from Fig. 8 that
the proposed three weight methods exceeds the learning output of RTS without
using TDS in terms of variance performance. Due to the fact that the learning
output of TDS is normalized to the same magnitude of the estimation output of
RTS, the element numbers of the two sets TDS and RTS will play a dominant role
in the weighting process. The result can be seen from Fig. 8. That is, the
proposed method II is the best one among the proposed three methods.

Fig. 9 plots the curves SR versus SNR of our proposed robust synthesis method
for DM using HAD structure with SNRTDS=-5dB . Here, the corresponding non-
robust method is used as a performance benchmark and the number of RF chains
is chosen as follows: $N_{RF}=\\{4,~{}8,~{}16\\}$. As can be seen from Fig. 9,
when DOA measurement error exists, the proposed robust DM method always
performs better than the corresponding non-robust one especially in the medium
and high SNR regions. In particular, with decrease in the number of RF chains,
the performance gain of the proposed robust method over nonrobust one become
more obvious. In such a situation, reducing the number of RF chains, under the
condition that other parameters are fixed like the receive SNR and a total
number of antennas at Alice, the DOA measurement error become large due to a
smaller number of spatial-time sampling points available. In other words, the
performance gain achieved by exploiting the density distribution of DOA
measurement error becomes more obvious in such a scenario. Otherwise, as the
number of RF chains increases, a larger number of spatial-time sampling points
is available to compute a high-precision variance matrix in ESPRIT. A higher-
precision DOA estimation is gotten. This implies a smaller measurement error,
which will accordingly decrease the performance gain achieved by the robust
one proposed by us over nonrobust one.

Fig. 10 shows the BER performance versus SNR of our proposed robust hybrid DM
with the corresponding non-robust one as a performance reference. Here,
$N=64$, and $M=\in\\{4,~{}8,~{}16\\}$. As shown in Fig. 10, we can find that
the proposed robust method achieves a better BER performance compared with the
corresponding non-robust one for $N_{RF}=\in\\{4,~{}8,~{}16\\}$. As the number
$N_{RF}$ of RF chains increases, the BER performance difference between robust
and non-robust ones disappears. The performance tendency in this figure is
very similar to that in Fig. 9.

##  VI Conclusion

In this paper, an I-HAD-ESPRIT DOA measurement method was proposed for HAD
receiver. Here, the phase ambiguity from HAD structure is addressed
successfully. The proposed I-HAD-ESPRIT can achieve the HAD CRLB. In the mean
time, we also proposed a SNR estimator. By using the proposed I-HAD-ESPRIT
method and histogram method in ML, we find an important fact that the measured
DOA or DOAME follows a Gaussian distribution. Subsequently, by applying the
proposed I-HAD-ESPRIT method repeatedly, we established two data sets: TDS and
RTS. After the maximum likelihood learning method was optimized over TDS and
RTS, we obtained the two high-precision learning outputs of DOA and variance.
To combine the two learning outputs and improve performance, three weight
methods was proposed to exploit the SNRs, numbers of TDS and RTS, or both. For
DOA, the proposed weight method III of exploiting both is the best one among
the three methods whereas the proposed weight method II is the best one among
the three methods. Finally, based on the above, a robust synthesis method for
HAD-based DM transmitter was proposed to exploit the density distribution of
DOA measurement error, with emphasis on the design of robust analog
beamforming vector. Compared to the corresponding non-robust method, the
proposed robust method can achieve a substantial SR performance gain.

.

## Appendix A: Derivation of Estimators of DOA mean and variance

From Section III, the measured value $\tilde{\theta}$ of DOA is shown to be a
Gaussian random variable with mean $\theta$ and variance $\sigma^{2}$, we have
its PDF as follows

$\displaystyle
f\left(\tilde{\theta}(n);\theta,\sigma^{2}\right)=\frac{1}{\sqrt{2\pi\sigma^{2}}}\exp\\{-\frac{1}{2\sigma^{2}}(\tilde{\theta}(n)-\theta)^{2}\\}$
(110)

Given $N_{DOA}$ measured values of angle of Bob:
$\left(\tilde{\theta}_{1},\cdots,\tilde{\theta}_{N_{DOA}}\right)$, we have
their likelihood

$\displaystyle
f\left(\tilde{\theta}_{1},\cdots,\tilde{\theta}_{N_{DOA}};\theta,\sigma^{2}\right)$
$\displaystyle=(\frac{1}{2\pi\sigma^{2}})^{\frac{N}{2}}\exp\\{-\frac{1}{2\sigma^{2}}\sum_{i=0}^{N-1}(\tilde{\theta}(n)-\theta)^{2}\\},$
(111)

which gives the corresponding log-likelihood function

$\displaystyle\ln
f\left(\tilde{\theta}_{1},\cdots,\tilde{\theta}_{N_{DOA}};\theta,\sigma^{2}\right)$
$\displaystyle=-\frac{N_{DOA}}{2}\ln(2\pi\sigma^{2})-\frac{1}{2\sigma^{2}}\sum_{i=0}^{N_{DOA}-1}(\tilde{\theta}_{i}-\theta)^{2}.$
(112)

Then, we compute the first derivative of $\ln P(\mathbf{\theta_{r}}|\theta)$
with respect to $\theta$ as follows

$\displaystyle\frac{\partial\ln
f\left(\tilde{\theta}_{1},\cdots,\tilde{\theta}_{N_{DOA}};\theta,\sigma^{2}\right)}{\partial\theta}=\frac{1}{\sigma^{2}}\sum_{i=0}^{N_{DOA}-1}(\tilde{\theta}_{i}-\theta)$
(113)

which is set to zero, then we can obtain the estimated value of $\theta$

$\displaystyle\hat{\theta}=\frac{1}{N_{DOA}}\sum_{i=0}^{N_{DOA}}\tilde{\theta}_{i}$
(114)

In the same way, we set the first derivative of $\ln
P(\mathbf{\theta_{r}}|\theta)$ with respect to $\sigma^{2}$ to zero

$\displaystyle\frac{\partial\ln
f\left(\tilde{\theta}_{1},\cdots,\tilde{\theta}_{N_{DOA}};\theta,\sigma^{2}\right)}{\partial\sigma^{2}}$
$\displaystyle=-\frac{N_{DOA}-1}{2\sigma^{2}}+\frac{1}{2\sigma^{4}}\sum_{i=0}^{N_{DOA}-1}(\tilde{\theta}_{i}-\theta)^{2}=0,$
(115)

which yields

$\displaystyle\hat{\sigma^{2}}=\frac{1}{N_{DOA}}\sum_{i=1}^{N_{DOA}}(\tilde{\theta}_{i}-\theta)^{2}.$
(116)

Since the true value of $\theta$ is unknown, replacing $\theta$ with its
estimated value $\hat{\theta}$ yields the bootstrap solution

$\displaystyle\hat{\sigma^{2}}=\frac{1}{N_{DOA}}\sum_{i=1}^{N_{DOA}}(\tilde{\theta}_{i}-\hat{\theta})^{2}$
(117)

in accordance with [43]. However, the above solution is a biased estimator.
Instead, we prefer the following unbiased estimator

$\displaystyle\hat{\sigma^{2}}=\frac{1}{N_{DOA}-1}\sum_{i=1}^{N_{DOA}}(\tilde{\theta}_{i}-\hat{\theta})^{2}$
(118)

with the same form as that of (117).

.

## Appendix B: Derivation of RAB $\mathbf{\tilde{V}_{RF}}$

Considering the DOA measurement error, the nonzero element of
$\mathbf{V}_{RF,RAB}$ can be represented as

$\displaystyle\hat{v}_{k,m}$
$\displaystyle=\int_{-\Delta\theta_{\mathrm{{max}}}}^{\Delta\theta_{\mathrm{{max}}}}e^{-j\frac{2\pi
d}{\lambda}[(k-1)M+m-\frac{N+1}{2}]\cos(\hat{\theta}-\Delta\theta)}\cdot
p(\Delta\theta)d(\Delta\theta)$
$\displaystyle=\int_{-\Delta\theta_{\mathrm{{max}}}}^{\Delta\theta_{\mathrm{{max}}}}e^{-j\alpha_{k,m}\cos(\hat{\theta}-\Delta\theta)}\cdot
p(\Delta\theta)d(\Delta\theta)$
$\displaystyle=\int_{-\Delta\theta_{\mathrm{{max}}}}^{\Delta\theta_{\mathrm{{max}}}}e^{-j\alpha_{k,m}[\cos(\hat{\theta})\cos(\Delta\theta)+\sin(\hat{\theta})\sin(\Delta\theta)]}\cdot
p(\Delta\theta)d(\Delta\theta)$

where

$\displaystyle\alpha_{k,m}=\frac{2\pi d}{\lambda}[(k-1)M+m-\frac{N+1}{2}].$
(120)

By utilizing the second-order Taylor expansion, we can expand
$\cos(\Delta\theta)$ and $\sin(\Delta\theta)$ at point $\Delta\theta=0$ as

$\displaystyle\cos(\Delta\theta)\approx 1-\frac{1}{2}(\Delta\theta)^{2}$
$\displaystyle\sin(\Delta\theta)\approx\Delta\theta.$ (121)

Substituting (Appendix B: Derivation of RAB $\mathbf{\tilde{V}_{RF}}$) in
(LABEL:69) yields

$\displaystyle\hat{v}_{k,m}$
$\displaystyle=\int_{-\Delta\theta_{\mathrm{{max}}}}^{\Delta\theta_{\mathrm{{max}}}}e^{-j\alpha_{k,m}[\cos(\hat{\theta})\cos(\Delta\theta)+\sin(\hat{\theta})\sin(\Delta\theta)]}\cdot
p(\Delta\theta)d(\Delta\theta)$
$\displaystyle=\int_{-\Delta\theta_{\mathrm{{max}}}}^{\Delta\theta_{\mathrm{{max}}}}e^{-j\alpha_{k,m}[\cos(\hat{\theta})-\cos(\hat{\theta})\cdot\frac{1}{2}(\Delta\theta)^{2}+\sin(\hat{\theta})(\Delta\theta)]}p(\Delta\theta)d(\Delta\theta)$
$\displaystyle=\xi_{k,m}\int_{-\Delta\theta_{\mathrm{{max}}}}^{\Delta\theta_{\mathrm{{max}}}}e^{j\alpha_{k,m}[\cos(\hat{\theta})\cdot\frac{1}{2}(\Delta\theta)^{2}-\sin(\hat{\theta})(\Delta\theta)]}p(\Delta\theta)d(\Delta\theta)$
$\displaystyle=\xi_{k,m}\int_{-\Delta\theta_{\mathrm{{max}}}}^{\Delta\theta_{\mathrm{{max}}}}\\{\cos(\alpha_{k,m}\psi)+j\sin(\alpha_{k,m}\psi)\\}p(\Delta\theta)d(\Delta\theta)$
$\displaystyle=\zeta_{k,m}+j\eta_{k,m}$ (122)

where

$\displaystyle\xi_{k,m}=e^{-j\alpha_{k,m}\cos(\hat{\theta})},$ (123)

and

$\displaystyle\psi=[\cos(\hat{\theta})\cdot\frac{1}{2}(\Delta\theta)^{2}-\sin(\hat{\theta})(\Delta\theta)].$
(124)

Similarly, by utilizing the second-order Taylor expansion,
$\cos(\alpha_{k,m}\psi)$ can be represented as

$\displaystyle\cos(\alpha_{k,m}\psi)=cos(\alpha_{k,m}[\cos(\hat{\theta})\cdot\frac{1}{2}(\Delta\theta)^{2}-\sin(\hat{\theta})(\Delta\theta)])$
$\displaystyle\approx
1-\frac{1}{2}\alpha_{k,m}^{2}[\cos(\hat{\theta})\cdot\frac{1}{2}(\Delta\theta)^{2}-\sin(\hat{\theta})(\Delta\theta)]^{2}$
$\displaystyle=1-\frac{1}{2}\alpha_{k,m}^{2}$
$\displaystyle\cdot[\frac{1}{4}\cos(\hat{\theta})^{2}(\Delta\theta)^{4}+\sin(\hat{\theta})^{2}(\Delta\theta)^{2}-\cos(\hat{\theta})\sin(\hat{\theta})\cdot(\Delta\theta)^{3}]$
$\displaystyle=1-\frac{1}{8}\alpha_{k,m}^{2}\cos(\hat{\theta})^{2}(\Delta\theta)^{4}-\frac{1}{2}\alpha_{k,m}^{2}\sin(\hat{\theta})^{2}(\Delta\theta)^{2}$
$\displaystyle+\frac{1}{2}\alpha_{k,m}^{2}\cos(\hat{\theta})\sin(\hat{\theta})(\Delta\theta)^{3}$
(125)

Since the last term of (Appendix B: Derivation of RAB
$\mathbf{\tilde{V}_{RF}}$) is an odd function of $\Delta\theta$, then

$\displaystyle\int_{-\Delta\theta_{\mathrm{{max}}}}^{\Delta\theta_{\mathrm{{max}}}}\frac{1}{2}\alpha_{k,m}^{2}\cos(\hat{\theta})\sin(\hat{\theta})(\Delta\theta)^{3}\cdot
p(\Delta\theta)d(\Delta\theta)=0.$ (126)

Now, let us define

$\displaystyle\chi_{1}=\int_{-\Delta\theta_{\mathrm{{max}}}}^{\Delta\theta_{\mathrm{{max}}}}(\Delta\theta)^{4}\cdot
p(\Delta\theta)d(\Delta\theta)$
$\displaystyle=\frac{2}{K_{d}\sqrt{2\pi\sigma^{2}}}\\{-\sigma^{2}\cdot\Delta\theta_{\mathrm{{max}}}^{3}\
e^{-\frac{\Delta\theta_{\mathrm{{max}}}^{2}}{2\sigma^{2}}}-3\sigma^{4}\Delta\theta_{\mathrm{{max}}}e^{-\frac{\Delta\theta_{\mathrm{{max}}}^{2}}{2\sigma^{2}}}$
$\displaystyle+\frac{3\sqrt{2\pi}}{2}\sigma^{5}erf(\frac{\Delta\theta_{\mathrm{{max}}}}{\sqrt{2}\sigma})\\},$
(127)

and

$\displaystyle\chi_{2}=\int_{-\Delta\theta_{\mathrm{{max}}}}^{\Delta\theta_{\mathrm{{max}}}}(\Delta\theta)^{2}\cdot
p(\Delta\theta)d(\Delta\theta)$
$\displaystyle=\frac{2}{K_{d}\sqrt{2\pi\sigma^{2}}}\\{-\sigma^{2}\Delta\theta_{\mathrm{{max}}}\cdot
e^{-\frac{\Delta\theta_{\mathrm{{max}}}^{2}}{2\sigma^{2}}}+\frac{\sqrt{2\pi}}{2}\sigma^{3}erf(\frac{\Delta\theta_{\mathrm{{max}}}}{\sqrt{2}\sigma})\\}.$
(128)

Then the real part $\zeta_{k,m}$ of (Appendix B: Derivation of RAB
$\mathbf{\tilde{V}_{RF}}$) can be expressed as

$\displaystyle\zeta_{k,m}=\xi_{k,m}(K_{d}-\frac{1}{8}\alpha_{k,m}^{2}\cos(\hat{\theta})^{2}\chi_{1}-\frac{1}{2}\alpha_{k,m}^{2}\sin(\hat{\theta})^{2}\chi_{2})$
(129)

In the same manner, we can have the $\eta_{k,m}$ as follows

$\displaystyle\eta_{k,m}=\frac{1}{2}\xi_{k,m}\cdot\alpha_{k,m}cos(\hat{\theta})\chi_{2}$
(130)

Until now, we complete the derivation of $\hat{v}_{k,m}$. That is,
$\hat{v}_{k,m}=\zeta_{k,m}+j\eta_{k,m}$. Due to the special structure of
analog precoder, we only need the phase of $\hat{v}_{k,m}$. Therefore we can
reformulated the analog $v_{k,m}$ as

$\displaystyle v_{RAB,k,m}=\frac{1}{\sqrt{M}}\exp(j*\angle(\hat{v}_{k,m})).$
(131)

which is what we need.

## References

  * [1] L. Godara, “Application of antenna arrays to mobile communications. II. Beam-forming and direction-of-arrival considerations,” _Proc. IEEE_ , vol. 85, no. 8, pp. 1195–1245, Aug. 1997.
  * [2] I. F. Gorodnitsky and B. D. Rao, “Sparse signal reconstruction from limited data using FOCUSS: a re-weighted minimum norm algorithm,” _IEEE Trans. Signal Process._ , vol. 45, no. 3, pp. 600–616, Mar. 1997.
  * [3] J. C. Chen, K. Yao, and R. E. Hudson, “Source localization and beamforming,” _IEEE Signal Process. Mag._ , vol. 19, no. 2, pp. 30–39, Mar. 2002.
  * [4] P. Stoica, P. Babu, and J. Li, “SPICE: A sparse covariance-based estimation method for array processing,” _IEEE Trans. Signal Process._ , vol. 59, no. 2, pp. 629–638, Feb. 2011.
  * [5] Z. Yang, L. Xie, and C. Zhang, “Off-grid direction of arrival estimation using sparse bayesian inference,” _IEEE Trans. Signal Process._ , vol. 61, no. 1, pp. 38–43, Jan. 2013.
  * [6] A. Hassanien and S. A. Vorobyov, “Transmit energy focusing for DOA estimation in MIMO radar with colocated antennas,” _IEEE Trans. Signal Process._ , vol. 59, no. 6, pp. 2669–2682, Jun. 2011.
  * [7] X. Zhang, L. Xu, L. Xu, and D. Xu, “Direction of departure (DOD) and direction of arrival (DOA) estimation in MIMO radar with reduced-dimension MUSIC,” _IEEE Commun. Lett._ , vol. 14, no. 12, pp. 1161–1163, Dec. 2010.
  * [8] R. Shafin, L. Liu, J. Zhang, and Y. Wu, “DOA estimation and capacity analysis for 3-D millimeter wave massive-MIMO/FD-MIMO OFDM systems,” _IEEE Trans. Wireless Commun._ , vol. 15, no. 10, pp. 6963–6978, Oct. 2016\.
  * [9] L. Wan, G. Han, J. Jiang, J. J. P. C. Rodrigues, N. Feng, and T. Zhu, “DOA estimation for coherently distributed sources considering circular and noncircular signals in massive MIMO systems,” _IEEE Syst. J._ , vol. 11, no. 1, pp. 41–49, Mar. 2017.
  * [10] H. Huang, J. Yang, H. Huang, Y. Song, and G. Gui, “Deep learning for super-resolution channel estimation and DOA estimation based massive MIMO system,” _IEEE Trans. Veh.Technol._ , vol. 67, no. 9, pp. 8549–8560, Aug. 2018.
  * [11] T. E. Tuncer and B. Friedlander, “Classical and modern direction-of-arrival estimation,” _New York, NY, USA:Elsevier_ , 2009.
  * [12] R. O. Schmidt, “Multiple emitter location and signal parameter estimation,” _IEEE Trans. Antennas Propag._ , vol. AP-34, no. 3, pp. 276–280, Mar. 1986\.
  * [13] R. Roy and T. Kailath, “ESPRIT-estimation of signal parameters via rotational invariance techniques,” _IEEE Trans. Acoust. Speech, Signal Processing_ , vol. 37, no. 7, pp. 984–995, July. 1989.
  * [14] D. Malioutov, M. Cetin, and A. S. Willsky, “A sparse signal reconstruct on perspective for source localization with sensor arrays,” _IEEE Trans. Signal Process._ , vol. 53, no. 8, pp. 3010–3022, Aug. 2005.
  * [15] M. M. Hyder and K. Mahata, “Direction-of-arrival estimation using a mixed$\ell_{2,0}$norm approximation,” _IEEE Trans. Signal Process._ , vol. 58, no. 9, pp. 4646–4655, Sep. 2010.
  * [16] Q. Li, X. Zhang, and H. Li, “Online direction of arrival estimation based on deep learning,” _IEEE International Conference on Acoustics, Speech and Signal Processing(ICASSP)_ , Apr 2018.
  * [17] S. Chakrabarty and E. A. P. Habets, “Broadband doa estimation using convolutional neural networks trained with noise signals,” _IEEE Workshop on Applications of Signal Processing to Audio and Acoustics(WASPAA)_ , pp. 1–3, Oct 2017.
  * [18] F. Shu, Y. Qin, T. Liu, L. Gui, Y. Zhang, J. Li, and Z. Han, “Low-complexity and high-resolution DOA estimation for hybrid analog and digital massive MIMO receive array,” _IEEE Trans. Commun._ , vol. 66, no. 6, pp. 2487–2501, June. 2018.
  * [19] M. R. M. Maina, D. K. Lang at, and D. P. K. Kihato, “An analysis of MUSIC, ESPRIT and root-MUSIC direction of arrival estimation techniques in smart antennas,” _Proceedings of the Sustainable Research and Innovation Conference_ , May. 2017.
  * [20] N. D. Sidiropoulos, R. Bro, and G. B. Giannakis, “Parallel factor analysis in sensor array processing,” _IEEE Trans. Signal Process._ , vol. 48, no. 8, pp. 2377–2388, 2000.
  * [21] A. Babakhani, D. B. Rutledge, and A. Hajimiri, “Transmitter architectures based on near-field direct antenna modulation,” _IEEE Journal of Solid-State Circuits_ , vol. 43, no. 12, pp. 2674 – 2692, Dec. 2008.
  * [22] M. P. Daly and J. T. Bernhard, “Directional modulation technique for phased arrays,” _IEEE Trans. Antennas Propag._ , vol. 57, no. 9, pp. 2633–2640, Sep. 2009.
  * [23] H. Shi and A. Tennant, “Enhancing the security of communication via directly modulated antenna arrays,” _IET Microwaves, Antennas and Propagation_ , vol. 7, no. 8, pp. 606–611, June. 2013.
  * [24] Y. Ding and V. F. Fusco, “A vector approach for the analysis and synthesis of directional modulation transmitters,” _IEEE Trans. Antennas Propag._ , vol. 62, no. 1, p. 361 C370, Jan. 2014.
  * [25] J. Hu, F. Shu, and J. Li, “Robust synthesis method for secure directional modulation with imperfect direction angle,” _IEEE Commun. Lett._ , vol. 20, no. 6, pp. 1084–1087, Jun. 2016.
  * [26] F. Shu, X. Wu, J. Li, R. Chen, and B. Vucetic, “Robust beamforming scheme for secure multi-beam directional modulation in broadcasting systems,” _IEEE Access_ , vol. 4, pp. 6614–6623, Oct. 2016.
  * [27] F. Shu, W. Zhu, X. Zhou, J. Li, and J. Lu, “Robust secure transmission of using main-lobe-integration-based leakage beamforming in directional modulation MU-MIMO systems,” _IEEE Syst. J._ , vol. 99, pp. 1–11, Nov. 2017.
  * [28] W. Zhu, F. Shu, T. Liu, and X. Zhou, “Secure precise transmission with multi-relay-aided directional modulation,” _2017 9th International Conference on Wireless Communications and Signal Processing (WCSP)_ , vol. 99, pp. 1–5, Dec. 2017.
  * [29] X. Zhou, J. Li, F. Shu, Q. Wu, and Y. Wu, “Secure swipt for directional modulation aided af relaying networks,” _arXiv preprint_ , vol. arXiv:1803.05278, 2018.
  * [30] J. Hu, S. Yan, F. Shu, J. Wang, J. Li, and Y. Zhang, “Artificial-noise-aided secure transmission with directional modulation based on random frequency diverse arrays,” _IEEE Access_ , vol. 5, pp. 1658–1667, Jan. 2017.
  * [31] F. Shu, X. Wu, J. Hu, J. Li, R. Chen, and J. Wang, “Secure and precise wireless transmission for random-subcarrier-selection-based directional modulation transmit antenna array,” _IEEE J. Sel. Areas Commun._ , vol. 36, no. 4, pp. 890–904, April. 2018.
  * [32] S., F. Shu, J. Lu, G. Gui, J. Wang, G. Xia, Y. Zhang, J. Li, and J. Wang, “Power allocation strategy of maximizing secrecy rate for secure directional modulation networks,” _IEEE Access._ , vol. 6, pp. 38 794–38 801, 2018\.
  * [33] H., S. Wan, W. Cai, L. Xu, X. Zhou, J. Wang, Y. Wu, F. Shu, J. Wang, and J. Wang, “GPI-based secrecy rate maximization beamforming scheme for wireless transmission with AN-aided directional modulation,” _IEEE Access._ , vol. 6, pp. 12 044–12 051, 2018.
  * [34] X. Zhang, A. F. Molisch, and S.-Y. Kung, “Variable-phase-shift-based rf-baseband codesign for MIMO antenna selection,” _IEEE Trans. Signal Process._ , vol. 53, no. 11, pp. 4091–4103, Nov. 2005.
  * [35] F. Sohrabi and W. Yu, “Hybrid analog and digital beamforming for mmwave OFDM large-scale antenna arrays,” _IEEE J. Sel. Areas Commun._ , vol. 35, no. 7, pp. 1432–1443, July 2017.
  * [36] ——, “Hybrid digital and analog beamforming design for large-scale antenna arrays,” _IEEE J. Sel. Topics Signal Process._ , vol. 10, no. 3, pp. 501–513, 2016.
  * [37] X. Yu, J. Shen, J. Zhang, and K. B. Letaief, “Alternating minimization algorithms for hybrid precoding in millimeter wave MIMO systems,” _IEEE J. Sel. Topics Signal Process._ , vol. 10, no. 3, pp. 485–500, Apr. 2016.
  * [38] X. Gao, L. Dai, S. Han, C. I, and R. W. Heath, “Energy-efficient hybrid analog and digital precoding for mmwave MIMO systems with large antenna arrays,” _IEEE J. Sel. Areas Commun._ , vol. 34, no. 4, pp. 998–1009, Apr. 2016.
  * [39] Y. R. Ramadan, H. Minn, and A. S. Ibrahim, “Hybrid analog cdigital precoding design for secrecy mmwave MISO-OFDM systems,” _IEEE Trans. Commun._ , vol. 65, no. 11, pp. 5009–5026, Nov. 2017.
  * [40] R. W. Heath, N. Gonz lez-Prelcic, S. Rangan, W. Roh, and A. M. Sayeed, “An overview of signal processing techniques for millimeter wave MIMO systems,” _IEEE J. Sel. Topics Signal Process._ , vol. 10, no. 3, pp. 436–453, 2016.
  * [41] R. A. Horn and C. R. Johnson, _Pattern Recognition And Machine Learning_. New York, USA: Springer Science+Business Mesia, LLC, 2013.
  * [42] S. Feng, W. Siming, Y. Shihao, W. Qian, W. Yongpeng, C. Riqing, L. Jun, and L. Jinhui, “Secure directional modulation to enhance physical layer security in iot networks,” _arXiv preprint_ , vol. arXiv:1712.02104, 2018.
  * [43] L. Wasserman, _All of Statistics_. Springer New York: Springer Texts in Statistics, 2003.