adaptive_beamforming.py

"""
Developed by Vladimir Fadeev
(https://github.com/kirlf)
Kazan, 2017 / 2020

Python 3.7

The result is uploaded in 
https://commons.wikimedia.org/wiki/File:AdaptiveBeamForming.png
"""

import numpy as np
import matplotlib.pyplot as plt

"""
Received signal model:
X = A*S + N

where 

A = [a(theta_1) a(theta_2) ... a(theta_d)] 
is the matrix of steering vectors 
(dimension is M x d, 
M is the number of sensors, 
d is the number of signal sources),

A steering vector represents the set of phase delays 
a plane wave experiences, evaluated at a set of array elements (antennas). 
The phases are specified with respect to an arbitrary origin.

theta is Direction of Arrival (DoA), 

S = 1/sqrt(2) * (X + iY)
is the transmit (modulation) symbols matrix 
(dimension is d x T, 
T is the number of snapshots)
(X + iY) is the complex values of the signal envelope,

N = sqrt(N0/2)*(G1 + jG2)
is additive noise matrix (AWGN)
(dimension is M x T),
N0 is the noise spectral density,
G1 and G2 are the random Gussian distributed values.

"""

M = 9 # number of antenna elements (sensors)

""" Correlation matrix of the information symbols:
Rss = S*S^H = I_d (try with QPSK, for example) """
Rss = np.eye((2)) 

""" Correlation matrix of additive noise:
Rnn = N*N^H = sigma_N^2 * I_M), 
where sigma_N^2 is the noise variance (power) """
Rnn = 0.1*np.eye((M)) # 

""" Let us consider 2 sources of the signals """
theta_1 = 0*(np.pi/180)  
theta_2 = 50*(np.pi/180)

""" Spatial frequency (some equivalent of a DoA):
mu = (2*pi / lambda ) * delta * sin(theta)

where 

delta is the antenna spacing 
(distance between antenna elements), and

lambda is the electromagnetic wave length.

Let us (delta = lambda / 2) then: 

"""
mu_1 = np.pi*np.sin(theta_1)
mu_2 = np.pi*np.sin(theta_2)

""" Steering vectors """
a_1 = np.exp(1j*mu_1*np.arange(M))
a_2 = np.exp(1j*mu_2*np.arange(M))
A = (np.array([a_1, a_2])).T

""" 
Correlation matrix of the received signals 
R_xx = X*X^H = A*R_ss*A^H + R_nn 
"""
R = A @ Rss @ np.conj(A).T + Rnn


""" Let us theta_1 is the signal, and theta_2 is the interferer """
g = np.array([1, 0]) # the first DoA is "switched on", the second DoA is "switched off". 

def calc_w_capon(A_i): 
    """ Capon's method (MVDR) 
    w_Capon = R^(-1) * A * (A^H * R^(-1) * A)^(-1) * g """
    w = (np.linalg.inv(R) @ A_i @ 
            np.linalg.inv( np.conj(A_i).T @ np.linalg.inv(R) @ A_i ) @ g).T 
    return w

def calc_power(w, a):
    """ P(theta) = |w_(opt)^H * a(theta)|^2 """
    P = (np.abs( (np.conj(w).T @ a) )**2).item()
    return P


""" Bartlett's method (сonventional beamforming)
w_Bart = a_1 / M """
w_bart = (a_1 / M).reshape((M,1))


""" Simulation loop.
Main idea: 
        1) We have the Rxx matrix from the receiver.
        2) We know the DoA of the information signal and
           DoA of the interference (e.g., based on frequency estimation methods)
        3) We should calculate optimal weight vector which will suppress interference.
        4) This should make SINR (Signal to Interference + Noise Ratio) better.
        5) Interference DoA can changes, but estimated Rxx should be the same!
"""
sinr_thetas = np.arange(1, 91)*(np.pi/180) # degrees (from 1 to 90) -> radians

SINR_Capon = np.empty(len(sinr_thetas), dtype = complex)
SINR_Bart = np.empty(len(sinr_thetas), dtype = complex)

for idx, theta_i in enumerate(sinr_thetas):
    
    """ Let's try to simulate changing of the interference picture! 
    For this redefine DoA of intereference. """
    mu_2 = np.pi*np.sin(theta_i)
    a_2 = np.exp(1j*mu_2*np.arange(M))
    
    A_sinr = (np.array([a_1, a_2])).T
    
    """ Capon's (MVDR) method: """
    w_capon = calc_w_capon(A_sinr)

    signal_capon = calc_power(w_capon, a_1)
    interf_capon = calc_power(w_capon, a_2)

    """ P_noise = w^H * Rnn * w """
    noise_capon = (np.conj(w_capon).T @ Rnn @ w_capon).item()
    
    """ SINR - Signal to Interference + Noise Ratio """
    SINR_Capon[idx] =  signal_capon / (interf_capon + noise_capon)

    """ Bartlett's method 
    (uses the same weight vector for every cases - not adaptive): """
    signal_bart = calc_power(w_bart, a_1)
    interf_bart = calc_power(w_bart, a_2)
    noise_bart = (np.conj(w_bart).T @ Rnn @ w_bart).item()
    SINR_Bart[idx] =  signal_bart / (interf_bart + noise_bart)

"""
Capon's method is more stable, 
Bartlett's method cannot well mitigate changed interference.
"""

plt.subplots(figsize=(10, 5), dpi=150)
plt.plot(sinr_thetas*(180/np.pi), 10*np.log10(np.real(SINR_Capon)), color='green', label='Capon')
plt.plot(sinr_thetas*(180/np.pi), 10*np.log10(np.real(SINR_Bart)), color='red', label='Bartlett')
plt.grid(color='r', linestyle='-', linewidth=0.2)
plt.xlabel('Azimuth angles θ (degrees)')
plt.ylabel('SINR (dB)')
plt.legend()
plt.show()

""" References
1. Haykin, Simon, and KJ Ray Liu. 
Handbook on array processing and sensor networks. 
Vol. 63. John Wiley & Sons, 2010. pp. 102-107


2. Haykin, Simon S. 
Adaptive filter theory. 
Pearson Education India, 2008. pp. 422-427

3. Richmond, Christ D. 
"Capon algorithm mean-squared error threshold 
SNR prediction and probability of resolution." IEEE 
"""

mvdr_music.py

"""
Developed by Vladimir Fadeev
(https://github.com/kirlf)

Kazan, 2017 / 2020
Python 3.7

The result is uploaded in 
https://commons.wikimedia.org/wiki/File:MUSIC_MVDR.png
"""

import numpy as np
import matplotlib.pyplot as plt

"""
Received signal model:

	X = A*S + W

where 
A = [a(theta_1) a(theta_2) ... a(theta_d)] 
is the matrix of steering vectors 
(dimension is M x d, 

M is the number of sensors, 

d is the number of signal sources),

A steering vector represents the set of phase delays 
a plane wave experiences, evaluated at a set of array elements (antennas). 

The phases are specified with respect to an arbitrary origin.
theta is Direction of Arrival (DoA), 

S = 1/sqrt(2) * (X + iY)
is the transmit (modulation) symbols matrix 
(dimension is d x T, 

T is the number of snapshots)
(X + iY) is the complex values of the signal envelope,

W = sqrt(N0/2)*(G1 + jG2)
is additive noise matrix (AWGN)
(dimension is M x T),

N0 is the noise spectral density,

G1 and G2 are the random Gaussian distributed values.
"""

M = 10 # number of sensors 
SNR = 10 # Signal-to-Noise ratio (dB) 
d = 3 # number sources of EM waves
N = 50 # number of snapshots

""" Signal matrix """

S = ( np.sign(np.random.randn(d,N)) + 1j * np.sign(np.random.randn(d,N)) ) / np.sqrt(2) # QPSK

""" Noise matrix 

Common formula:

	AWGN = sqrt(N0/2)*(G1 + jG2), 

where G1 and G2 - independent Gaussian processes.

Since Es(symbol energy) for QPSK is 1 W, noise spectral density: 
	
	N0 = (Es/N)^(-1) = SNR^(-1) [W] (let SNR = Es/N0); 

or in logarithmic scale::

	SNR_dB = 10log10(SNR) => N0_dB = -10log10(SNR) = -SNR_dB [dB]; 

We have SNR in logarithmic (in dBs), convert to linear:
	SNR = 10^(SNR_dB/10) => sqrt(N0) = (10^(-SNR_dB/10))^(1/2) = 10^(-SNR_dB/20) 
"""

W = ( np.random.randn(M,N) + 1j * np.random.randn(M,N) ) / np.sqrt(2) * 10**(-SNR/20) # AWGN

mu_R = 2*np.pi / M  # standard beam width


cases = [[-1., 0, 1.], [-0.5, 0, 0.5], [-0.3, 0, 0.3],] # resolutions 
for idxm, c in enumerate(cases):
    
    """ DoA (spatial frequencies) """
    mu_1 = c[0]*mu_R
    mu_2 = c[1]*mu_R
    mu_3 = c[2]*mu_R

    """ Steering vectors """
    a_1 = np.exp(1j*mu_1*np.arange(M))
    a_2 = np.exp(1j*mu_2*np.arange(M))
    a_3 = np.exp(1j*mu_3*np.arange(M))

    A = (np.array([a_1, a_2, a_3])).T # steering matrix 
    
    """ Received signal """
    X = np.dot(A,S) + W 

    """ Rxx """
    R = np.dot(X,np.matrix(X).H)

    U, Sigma, Vh = np.linalg.svd(X, full_matrices=True)
    U_0 = U[:,d:] # noise sub-space

    thetas = np.arange(-90,91)*(np.pi/180) # azimuths
    mus = np.pi*np.sin(thetas) # spatial frequencies
    
    a = np.empty((M, len(thetas)), dtype = complex)
    for idx, mu in enumerate(mus):
        a[:,idx] = np.exp(1j*mu*np.arange(M))

    # MVDR:
    S_MVDR = np.empty(len(thetas), dtype = complex)
    for idx in range(np.shape(a)[1]):
        a_idx =  (a[:, idx]).reshape((M, 1))
        S_MVDR[idx] = 1 / (np.dot(np.matrix(a_idx).H, np.dot(np.linalg.pinv(R),a_idx)))

    # MUSIC:
    S_MUSIC = np.empty(len(thetas), dtype = complex)
    for idx in range(np.shape(a)[1]):
        a_idx =  (a[:, idx]).reshape((M, 1))
        S_MUSIC[idx] = np.dot(np.matrix(a_idx).H,a_idx)\
        / (np.dot(np.matrix(a_idx).H, np.dot(U_0,np.dot(np.matrix(U_0).H,a_idx))))

    plt.subplots(figsize=(10, 5), dpi=150)
    plt.semilogy(thetas*(180/np.pi), np.real( (S_MVDR / max(S_MVDR))), color='green', label='MVDR')
    plt.semilogy(thetas*(180/np.pi), np.real((S_MUSIC/ max(S_MUSIC))), color='red', label='MUSIC')
    plt.grid(color='r', linestyle='-', linewidth=0.2)
    plt.xlabel('Azimuth angles (degrees)')
    plt.ylabel('Power (pseudo)spectrum (normalized)')
    plt.legend()
    plt.title('Case #'+str(idxm+1))
    plt.show()

""" References

1. Haykin, Simon, and KJ Ray Liu. Handbook on array processing and sensor networks. Vol. 63. John Wiley & Sons, 2010. pp. 102-107
2. Hayes M. H. Statistical digital signal processing and modeling. – John Wiley & Sons, 2009.
3. Haykin, Simon S. Adaptive filter theory. Pearson Education India, 2008. pp. 422-427
4. Richmond, Christ D. "Capon algorithm mean-squared error threshold SNR prediction and probability of resolution." IEEE Transactions on Signal Processing 53.8 (2005): 2748-2764.
5. S. K. P. Gupta, MUSIC and improved MUSIC algorithm to esimate dorection of arrival, IEEE, 2015.
"""

standard_esprit.py

import numpy as np
import scipy.linalg as LA
import matplotlib.pyplot as plt

def standard_esprit(Rxx, d, M):
    """
    Standard ESPRIT calulations (max-overlapping case)
    Derived from: 
    https://github.com/dengjunquan/DoA-Estimation-MUSIC-ESPRIT/blob/master/DoAEstimation.py
    
    Parameters
    __________
    
    Rxx: 2-D array
        Covariance matrix of the received signal
    d: int
        Number of signal sources
    M: int
        Number of sensors
        
    Returns
    _______
    
    estimated: 1-D array of floats
    
        Estimated frequencies
    """
    _, U = LA.eig(Rxx) # SVD
    Us = U[:, :d] # signal sub-space
    # Selection matrices: J1 = [Im 0] and J2 = [0 Im]:
    phis = LA.pinv(Us[:M-1]) @ Us[1:M] # invariance equation (least square solution)
    eigs, _ = LA.eig(phis) # eigen values
    estimated = np.angle(eigs) # frequency estimation
    return estimated
  
M = 10 # number of sensors
SNR = 10 # Signal-to-Noise ratio (dB) 
d = 3 # number of souces
N = 50 # number of snapshots

S = ( np.sign(np.random.randn(d, N)) + 1j * np.sign(np.random.randn(d, N)) ) / np.sqrt(2) # QPSK
W = ( np.random.randn(M, N) + 1j * np.random.randn(M, N) ) / np.sqrt(2) * 10**(-SNR / 20) # AWGN

mu_R = 2*np.pi / M 
mu_1, mu_2, mu_3 = [item*mu_R for item in [-1., 0, 1.]] # directions of arrival (DoAs)

# steering vectors
a_1 = np.exp(1j*mu_1*np.arange(M))
a_2 = np.exp(1j*mu_2*np.arange(M))
a_3 = np.exp(1j*mu_3*np.arange(M))

A = (np.array([a_1, a_2, a_3])).T # steering matrix
X = A @ S + W # received signal matrix

Rxx = X @ X.conjugate().T / N # covariance matrix

print('Actual DoAs:', np.sort(standard_esprit(Rxx, d, M)),'\n')
print('ESPRIT DoAs:', np.sort(res),'\n')