18. DOA & Beamforming

En este cap铆tulo cubrimos los conceptos de beamforming, direcci贸n de llegada (DOA) y arreglos en fase en general. Se analizan t茅cnicas como Capon y MUSIC, utilizando ejemplos de simulaci贸n de Python. Cubrimos beamforming frente a DOA y analizamos los dos tipos diferentes de matrices en fase (matriz pasiva escaneada electr贸nicamente o PESA y matriz activa escaneada electr贸nicamente o AESA).

Descripci贸n general y terminolog铆a

Un Phase Array, tambi茅n conocido como conjunto dirigido electr贸nicamente, es un conjunto de antenas que se pueden usar en el lado de transmisi贸n o recepci贸n para formar haces en una o m谩s direcciones deseadas. Se utilizan tanto en comunicaciones como en radar, y los encontrar谩 en tierra, en aire y en sat茅lites.

Los arreglos en fase se pueden dividir en tres tipos:

  1. Passive electronically scanned array (PESA), tambi茅n conocida como matriz en fase anal贸gica o tradicional, donde se utilizan desfasadores anal贸gicos para dirigir el haz. En el lado de recepci贸n, todos los elementos se suman despu茅s del cambio de fase (y opcionalmente, ganancia ajustable) y se convierten en un canal de se帽al que se convierte y se recibe. En el lado de transmisi贸n ocurre lo contrario; se emite una 煤nica se帽al digital desde el lado digital, y en el lado anal贸gico se utilizan desfasadores y etapas de ganancia para producir la salida que va a cada antena. Estos desfasadores digitales tendr谩n un n煤mero limitado de bits de resoluci贸n y controlar谩n la latencia.

  2. Active electronically scanned array (AESA), tambi茅n conocidos como matrices totalmente digitales, donde cada elemento tiene su propia interfaz de RF y la formaci贸n del haz se realiza completamente en el dominio digital. Este es el enfoque m谩s caro, ya que los componentes de RF son caros, pero proporciona mucha m谩s flexibilidad y velocidad que los PESA. Los arreglos digitales son populares entre los SDR, aunque la cantidad de canales de recepci贸n o transmisi贸n del SDR limita la cantidad de elementos en su arreglo.

  3. Hybrid array, que consta de subconjuntos que se asemejan individualmente a los PESA, donde cada subconjunto tiene su propia interfaz de RF al igual que los AESA. Este es el enfoque m谩s com煤n para los arreglos en fase modernos, ya que ofrece lo mejor de ambos mundos.

A continuaci贸n se muestra un ejemplo de cada uno de los tipos.

Example of phased arrays including Passive electronically scanned array (PESA), Active electronically scanned array (AESA), Hybrid array, showing Raytheon's MIM-104 Patriot Radar, ELM-2084 Israeli Multi-Mission Radar, Starlink User Terminal, aka Dishy

En este cap铆tulo nos centramos principalmente en DSP para arreglos totalmente digitales, ya que son m谩s adecuados para la simulaci贸n y DSP, pero en el siguiente cap铆tulo nos ponemos manos a la obra con el arreglo 鈥淧haser鈥 y SDR de Analog Devices que tiene 8 desfasadores anal贸gicos. alimentando a un Pluto.

Normalmente nos referimos a las antenas que forman un conjunto como elementos y, a veces, al conjunto se le llama 鈥渟ensor鈥. Estos elementos del conjunto suelen ser antenas omnidireccionales, igualmente espaciadas en una l铆nea o en dos dimensiones.

Un formador de haces es esencialmente un filtro espacial; filtra se帽ales de todas las direcciones excepto las direcciones deseadas. En lugar de taps, utilizamos pesos (tambi茅n conocidos como coeficientes) aplicados a cada elemento de un arreglo. Luego manipulamos los pesos para formar los aces (s) del arreglo, 隆de ah铆 el nombre formador de haces (beamforming)! Podemos dirigir estos aces (y nulos) extremadamente r谩pido; debe ser m谩s r谩pido que las antenas con card谩n mec谩nico, que pueden considerarse una alternativa a los conjuntos en fase. Una sola matriz puede rastrear electr贸nicamente m煤ltiples se帽ales a la vez y al mismo tiempo anular las interferencias, siempre que tenga suficientes elementos. Normalmente analizaremos la formaci贸n de haces en el contexto de un enlace de comunicaciones, donde el receptor pretende recibir una o m谩s se帽ales con la SNR m谩s alta posible.

Los enfoques de beamforming suelen dividirse en convencionales y adaptativos. Con la formaci贸n de haces convencional se supone que ya se conoce la direcci贸n de llegada de la se帽al de inter茅s, y el formador de haces implica elegir pesos para maximizar la ganancia en esa direcci贸n. Esto se puede utilizar en el lado de recepci贸n o transmisi贸n de un sistema de comunicaci贸n. La formaci贸n de haz adaptativa, por otro lado, implica ajustar constantemente los pesos en funci贸n de la salida del formador de haz, para optimizar algunos criterios, lo que a menudo implica anular una fuente de interferencia. Debido al circuito cerrado y la naturaleza adaptativa, la formaci贸n de haz adaptativa generalmente solo se usa en el lado de recepci贸n, por lo que la 鈥渟alida del formador de haz鈥 es simplemente la se帽al recibida, y la formaci贸n de haz adaptativa implica ajustar los pesos en funci贸n de las estad铆sticas de los datos recibidos.

La direcci贸n de llegada (DOA) dentro de DSP/SDR se refiere al proceso de utilizar un conjunto de antenas para estimar las direcciones de llegada de una o m谩s se帽ales recibidas por ese conjunto (a diferencia de la formaci贸n de haces, que se centra en el proceso de recibir una se帽al mientras rechaza la mayor cantidad de ruido e interferencia). Aunque DOA ciertamente se incluye dentro del tema de beamforming, los t茅rminos pueden resultar confusos. Algunas t茅cnicas como MVDR/Capon se aplicar谩n tanto a DOA como a la formaci贸n de haces, porque la misma t茅cnica utilizada para la formaci贸n de haces se utiliza para realizar DOA barriendo el 谩ngulo de inter茅s y realizando la operaci贸n de formaci贸n de haces en cada 谩ngulo, luego buscando picos en el resultado ( cada pico es una se帽al, pero no sabemos si es la se帽al de inter茅s, una interferencia o incluso un rebote multitrayecto de la se帽al de inter茅s). Puede pensar en estas t茅cnicas DOA como una envoltura alrededor de un formador de haz espec铆fico. Existen t茅cnicas DOA como MUSIC y ESPIRT que son estrictamente para el prop贸sito de DOA. Debido a que la mayor铆a de las t茅cnicas de formaci贸n de haces suponen que usted conoce el 谩ngulo de llegada de la se帽al de inter茅s, si el objetivo se est谩 moviendo o el conjunto se est谩 moviendo, tendr谩 que realizar DOA continuamente como paso intermedio, incluso si su objetivo principal es recibir y demodular la se帽al de inter茅s.

Los arreglos en fase y el beamforming/DOA se utilizan en todo tipo de aplicaciones, aunque con mayor frecuencia los ver谩 utilizados en m煤ltiples formas de radar, comunicaci贸n mmWave dentro de 5G, comunicaciones por sat茅lite e interferencias. Cualquier aplicaci贸n que requiera una antena de alta ganancia, o que requiera una antena de alta ganancia que se mueva r谩pidamente, es buena candidata para los arreglos en fase.

Requerimientos SDR

Como se mencion贸, los arreglos en fase anal贸gicos implican un desfasador anal贸gico (y generalmente una ganancia ajustable) por canal, lo que significa que un arreglo en fase anal贸gico es una pieza de hardware dedicada que debe ir junto con un SDR. Por otro lado, cualquier SDR que contenga m谩s de un canal se puede utilizar como una matriz digital sin hardware adicional, siempre que los canales sean coherentes en fase y se muestreen utilizando el mismo reloj, lo que suele ser el caso de los SDR que tienen m煤ltiples canales. Recibir canales a bordo. Hay muchos SDR que contienen dos canales de recepci贸n, como el Ettus USRP B210 y Analog Devices Pluto (el segundo canal est谩 expuesto mediante un conector uFL en la propia placa). Desafortunadamente, ir m谩s all谩 de dos canales implica entrar en el segmento de SDR de m谩s de 10.000 d贸lares, al menos a partir de 2023, como el USRP N310. El principal problema es que los SDR de bajo costo normalmente no se pueden 鈥渆ncadenar鈥 para escalar el n煤mero de canales. La excepci贸n son KerberosSDR (4 canales) y KrakenSDR (5 canales), que utilizan m煤ltiples RTL-SDR que comparten un LO para formar una matriz digital de bajo costo; la desventaja es la frecuencia de muestreo (hasta 2,56 MHz) y el rango de sintonizaci贸n muy limitados (hasta 1766 MHz). A continuaci贸n se muestran la placa KrakenSDR y un ejemplo de configuraci贸n de antena.

The KrakenSDR

En este cap铆tulo no utilizamos ning煤n SDR espec铆fico; en lugar de eso, simulamos la recepci贸n de se帽ales usando Python y luego pasamos por el DSP y observamos el desempe帽o del Beamforming/DOA para arreglos digitales.

Introducci贸n a Matrix Math en Python/NumPy

Python tiene muchas ventajas sobre MATLAB, como ser gratuito y de c贸digo abierto, diversidad de aplicaciones, comunidad vibrante, 铆ndices que comienzan desde 0 como cualquier otro lenguaje, uso dentro de AI/ML y parece haber una biblioteca para cualquier cosa que se pueda imaginar. Pero lo que se queda corto es c贸mo se codifica/representa la manipulaci贸n de matrices (en t茅rminos de computaci贸n/velocidad, es bastante r谩pido, con funciones implementadas bajo un lenguaje m谩s eficiente como C/C++). No ayuda que haya m煤ltiples formas de representar matrices en Python, con el m茅todo np.matrix obsoleto en favor de np.ndarray. En esta secci贸n proporcionamos una breve introducci贸n a c贸mo hacer c谩lculos matriciales en Python usando NumPy, para que cuando lleguemos a los ejemplos de DOA se sienta m谩s c贸modo.

Comencemos saltando a la parte m谩s molesta de las matem谩ticas matriciales en NumPy; Los vectores se tratan como matrices 1D, por lo que no hay forma de distinguir entre un vector de fila y un vector de columna (se tratar谩 como un vector de fila de forma predeterminada), mientras que en MATLAB un vector es un objeto 2D. En Python puedes crear un nuevo vector usando a = np.array([2,3,4,5]) o convertir una lista en un vector usando mylist = [2, 3, 4 , 5] entonces a = np.asarray(mylist), pero tan pronto como quieras hacer c谩lculos matriciales, la orientaci贸n importa y estos se interpretar谩n como vectores de fila. Intentando hacer una transposici贸n en este vector, p.e. usando a.T, no lo cambiar谩 a un vector de columna. La forma de hacer un vector de columna a partir de un vector normal a es usar a = a.reshape(-1,1). El -1 le dice a NumPy que calcule el tama帽o de esta dimensi贸n autom谩ticamente, mientras mantiene la longitud de la segunda dimensi贸n 1. Lo que esto crea es t茅cnicamente una matriz 2D, pero la segunda dimensi贸n tiene una longitud 1, por lo que sigue siendo esencialmente 1D de una perspectiva matem谩tica. Es s贸lo una l铆nea adicional, pero realmente puede alterar el flujo del c贸digo matem谩tico matricial.

Ahora veamos un ejemplo r谩pido de matem谩ticas matriciales en Python; multiplicaremos una matriz 3x10 por una matriz 10x1. Recuerde que 10x1 significa 10 filas y 1 columna, lo que se conoce como vector de columna porque es solo una columna. Desde nuestros primeros a帽os escolares sabemos que esta es una multiplicaci贸n de matrices v谩lida porque las dimensiones internas coinciden y el tama帽o de la matriz resultante son las dimensiones externas, o 3x1. Usaremos np.random.randn() para crear 3x10 y np.arange() para crear 10x1, por conveniencia:

A = np.random.randn(3,10) # 3x10
B = np.arange(10) # 1D array of length 10
B = B.reshape(-1,1) # 10x1
C = A @ B # matrix multiply
print(C.shape) # 3x1
C = C.squeeze() # see next subsection
print(C.shape) # 1D array of length 3, easier for plotting and other non-matrix Python code

Despu茅s de realizar c谩lculos matriciales, es posible que el resultado se parezca a: [[ 0. 0.125 0.251 -0.376 -0.251 ...]] que claramente tiene solo una dimensi贸n de datos, pero si va a trazarlo, obtendr谩 un error o un gr谩fico que no muestra nada. Esto se debe a que el resultado es t茅cnicamente una matriz 2D y necesitas convertirla a una matriz 1D usando a.squeeze(). La funci贸n squeeze() elimina cualquier dimensi贸n de longitud 1 y resulta 煤til al realizar c谩lculos matriciales en Python. En el ejemplo anterior, el resultado ser铆a [ 0. 0.125 0.251 -0.376 -0.251 ...] (observe los segundos corchetes que faltan), que se puede trazar o usar en otro c贸digo Python que espere algo 1D .

Al codificar matem谩ticas matriciales, la mejor comprobaci贸n de cordura que puede hacer es imprimir las dimensiones (usando A.shape) para verificar que sean lo que espera. Considere pegar la forma en los comentarios despu茅s de cada l铆nea para referencia futura, de modo que sea f谩cil asegurarse de que las dimensiones coincidan al realizar multiplicaciones matriciales o por elementos.

A continuaci贸n se muestran algunas operaciones comunes tanto en MATLAB como en Python, como una especie de hoja de referencia a la que recurrir:

Operation

MATLAB

Python/NumPy

Create (Row) Vector, size 1 x 4

a = [2 3 4 5];

a = np.array([2,3,4,5])

Create Column Vector, size 4 x 1

a = [2; 3; 4; 5]; or a = [2 3 4 5].'

a = np.array([[2],[3],[4],[5]]) or
a = np.array([2,3,4,5]) then
a = a.reshape(-1,1)

Create 2D Matrix

A = [1 2; 3 4; 5 6];

A = np.array([[1,2],[3,4],[5,6]])

Get Size

size(A)

A.shape

Transpose a.k.a. A^T

A.'

A.T

Complex Conjugate Transpose
a.k.a. Conjugate Transpose
a.k.a. Hermitian Transpose
a.k.a. A^H

A'

A.conj().T

(unfortunately there is no A.H for ndarrays)

Elementwise Multiply

A .* B

A * B or np.multiply(a,b)

Matrix Multiply

A * B

A @ B or np.matmul(A,B)

Dot Product of two vectors (1D)

dot(a,b)

np.dot(a,b) (never use np.dot for 2D)

Concatenate

[A A]

np.concatenate((A,A))

Factor de Arreglo Matem谩tico

Para llegar a la parte divertida, tenemos que hacer un poco de matem谩ticas, pero la siguiente secci贸n se ha escrito para que las matem谩ticas sean extremadamente simples y tengan diagramas que las acompa帽en, solo se utilizan las propiedades trigonom茅tricas y exponenciales m谩s b谩sicas. . Es importante comprender las matem谩ticas b谩sicas detr谩s de lo que haremos en Python para realizar DOA.

Considere una matriz 1D de tres elementos espaciados uniformemente:

Diagram showing direction of arrival (DOA) of a signal impinging on a uniformly spaced antenna array, showing boresight angle and distance between elements or apertures

En este ejemplo, una se帽al llega desde el lado derecho, por lo que llega primero al elemento m谩s a la derecha. Calculemos el retraso entre el momento en que la se帽al llega a ese primer elemento y el momento en que llega al siguiente elemento. Podemos hacer esto formando el siguiente problema trigonom茅trico. Intenta visualizar c贸mo se form贸 este tri谩ngulo a partir del diagrama de arriba. El segmento resaltado en rojo representa la distancia que la se帽al tiene que recorrer despu茅s de haber alcanzado el primer elemento, antes de llegar al siguiente.

Trig associated with direction of arrival (DOA) of uniformly spaced array

Si recuerdas SOH CAH TOA, en este caso estamos interesados en el lado 鈥渁dyacente鈥 y tenemos la longitud de la hipotenusa (d), por lo que necesitamos usar un coseno:

\cos(90 - \theta) = \frac{\mathrm{adjacent}}{\mathrm{hypotenuse}}

Debemos resolver para la adyacente, ya que eso es lo que nos dir谩 qu茅 tan lejos debe viajar la se帽al entre el primer y el segundo elemento, para que se vuelva adyacente = d \cos(90 - \theta). Ahora hay una identidad trigonom茅trica que nos permite convertir esto en = d \sin(\theta) adyacente. Sin embargo, esto es solo una distancia, necesitamos convertir esto a un tiempo, usando la velocidad de la luz: tiempo transcurrido = d \sin(\theta) / c [segundos]. Esta ecuaci贸n se aplica entre cualquier elemento adyacente de nuestra matriz, aunque podemos multiplicar todo por un n煤mero entero para calcular entre elementos no adyacentes, ya que est谩n espaciados uniformemente (lo haremos m谩s adelante).

Ahora, para conectar estas matem谩ticas de activaci贸n y velocidad de la luz con el mundo del procesamiento de se帽ales. Denotemos nuestra se帽al de transmisi贸n en banda base s(t) y se est谩 transmitiendo en alguna portadora, f_c , por lo que la se帽al de transmisi贸n es s(t) e^{2j \pi f_ct}. Digamos que esta se帽al llega al primer elemento en el momento t = 0, lo que significa que llega al siguiente elemento despu茅s de d \sin(\theta) / c [segundos] como calculamos anteriormente. Esto significa que el segundo elemento recibe:

s(t - \Delta t) e^{2j \pi f_c (t - \Delta t)}

\mathrm{where} \quad \Delta t = d \sin(\theta) / c

recuerde que cuando tiene un cambio de tiempo, se resta del argumento de tiempo.

Cuando el receptor o SDR realiza el proceso de conversi贸n descendente para recibir la se帽al, esencialmente la multiplica por la portadora pero en la direcci贸n inversa, por lo que despu茅s de realizar la conversi贸n descendente el receptor ve:

s(t - \Delta t) e^{2j \pi f_c (t - \Delta t)} e^{-2j \pi f_c t}

= s(t - \Delta t) e^{-2j \pi f_c \Delta t}

Ahora podemos hacer un peque帽o truco para simplificar esto a煤n m谩s; considere c贸mo cuando tomamos muestras de una se帽al, se puede modelar sustituyendo t por nT donde T es el per铆odo de muestra y n es solo 0, 1, 2 , 3鈥 Sustituyendo esto obtenemos s(nT - \Delta t) e^{-2j \pi f_c \Delta t}. Bueno, nT es mucho mayor que \Delta t que podemos deshacernos del primer t茅rmino \Delta t y nos queda s(nT ) e^{-2j \pi f_c \Delta t}. Si la frecuencia de muestreo alguna vez llega a ser lo suficientemente r谩pida como para acercarse a la velocidad de la luz en una distancia peque帽a, podemos revisar esto, pero recuerde que nuestra frecuencia de muestreo solo necesita ser un poco mayor que el ancho de banda de la se帽al de inter茅s.

Sigamos con estos c谩lculos, pero comenzaremos a representar las cosas en t茅rminos discretos para que se parezca mejor a nuestro c贸digo Python. La 煤ltima ecuaci贸n se puede representar de la siguiente manera, volvamos a conectar \Delta t:

s[n] e^{-2j \pi f_c \Delta t}

= s[n] e^{-2j \pi f_c d \sin(\theta) / c}

Ya casi hemos terminado, pero afortunadamente hay una simplificaci贸n m谩s que podemos hacer. Recuerde la relaci贸n entre la frecuencia central y la longitud de onda: \lambda = \frac{c}{f_c} o la forma que usaremos: f_c = \frac{c}{\lambda}. Al conectar esto obtenemos:

s[n] e^{-2j \pi \frac{c}{\lambda} d \sin(\theta) / c}

= s[n] e^{-2j \pi d \sin(\theta) / \lambda}

En DOA lo que nos gusta hacer es representar d, la distancia entre elementos adyacentes, como una fracci贸n de longitud de onda (en lugar de metros), el valor m谩s com煤n elegido para d durante el proceso de dise帽o del arreglo. es utilizar la mitad de la longitud de onda. Independientemente de lo que sea d, a partir de este momento vamos a representar d como una fracci贸n de longitud de onda en lugar de metros, simplificando la ecuaci贸n y todo nuestro c贸digo:

s[n] e^{-2j \pi d \sin(\theta)}

Esto es para elementos adyacentes, para el k鈥櫭﹕imo elemento solo necesitamos multiplicar d por k:

s[n] e^{-2j \pi d k \sin(\theta)}

隆Y hemos terminado! 隆Esta ecuaci贸n anterior es lo que ver谩 en los documentos e implementaciones de DOA en todas partes! Normalmente llamamos a ese t茅rmino exponencial el 鈥渇actor de matriz鈥 (a menudo indicado como a) y lo representamos como una matriz, una matriz 1D para una matriz de antenas 1D, etc. En Python a es:

a = [np.exp(-2j*np.pi*d*0*np.sin(theta)), np.exp(-2j*np.pi*d*1*np.sin(theta)), np.exp(-2j*np.pi*d*2*np.sin(theta)), ...] # note the increasing k
# or
a = np.exp(-2j * np.pi * d * np.arange(Nr) * np.sin(theta)) # where Nr is the number of receive antenna elements

Observe c贸mo el elemento 0 da como resultado 1+0j (porque e^{0}=1); Esto tiene sentido porque todo lo anterior era relativo a ese primer elemento, por lo que recibe la se帽al tal como est谩 sin ning煤n cambio de fase relativa. As铆 es puramente como funcionan las matem谩ticas; en realidad, cualquier elemento podr铆a considerarse como la referencia, pero como ver谩 en nuestro c贸digo/matem谩tico m谩s adelante, lo que importa es la diferencia en fase/amplitud recibida entre los elementos. Todo es relativo.

Recibiendo se帽ales

Usemos el concepto de factor de matriz para simular una se帽al que llega a una matriz. Para una se帽al de transmisi贸n solo usaremos un tono por ahora:

import numpy as np
import matplotlib.pyplot as plt

sample_rate = 1e6
N = 10000 # number of samples to simulate

# Create a tone to act as the transmitter signal
t = np.arange(N)/sample_rate # time vector
f_tone = 0.02e6
tx = np.exp(2j * np.pi * f_tone * t)

Ahora simulemos un conjunto que consta de tres antenas omnidireccionales en una l铆nea, con 1/2 longitud de onda entre las adyacentes (tambi茅n conocido como 鈥渆spaciado de media longitud de onda鈥). Simularemos la se帽al del transmisor que llega a esta matriz en un cierto 谩ngulo, theta. Comprender el factor de matriz a a continuaci贸n es la raz贸n por la que revisamos todos los c谩lculos anteriores.

d = 0.5 # half wavelength spacing
Nr = 3
theta_degrees = 20 # direction of arrival (feel free to change this, it's arbitrary)
theta = theta_degrees / 180 * np.pi # convert to radians
a = np.exp(-2j * np.pi * d * np.arange(Nr) * np.sin(theta)) # array factor
print(a) # note that it's 3 elements long, it's complex, and the first element is 1+0j

Para aplicar el factor de matriz tenemos que hacer una multiplicaci贸n matricial de a y tx, as铆 que primero conviertamos ambos a 2D, usando el enfoque que discutimos anteriormente cuando revisamos c贸mo hacer matem谩ticas matriciales en Python. Comenzaremos convirtiendo ambos en vectores de fila usando x.reshape(-1,1). Luego realizamos la multiplicaci贸n matricial, indicada por el s铆mbolo @. Tambi茅n tenemos que convertir tx de un vector de fila a un vector de columna usando una operaci贸n de transposici贸n (imag铆nelo girando 90 grados) para que la matriz multiplique las dimensiones internas coincidan.

a = a.reshape(-1,1)
print(a.shape) # 3x1
tx = tx.reshape(-1,1)
print(tx.shape) # 10000x1

# matrix multiply
r = a @ tx.T  # dont get too caught up by the transpose, the important thing is we're multiplying the array factor by the tx signal
print(r.shape) # 3x10000.  r is now going to be a 2D array, 1D is time and 1D is the spatial dimension

En este punto r es una matriz 2D, tama帽o 3 x 10000 porque tenemos tres elementos de la matriz y 10000 muestras simuladas. Podemos extraer cada se帽al individual y trazar las primeras 200 muestras; a continuaci贸n trazaremos solo la parte real, pero tambi茅n hay una parte imaginaria, como cualquier se帽al de banda base. Una parte molesta de las matem谩ticas matriciales en Python es la necesidad de agregar .squeeze(), que elimina todas las dimensiones con longitud 1, para volver a una matriz NumPy 1D normal que se espera para el trazado y otras operaciones.

plt.plot(np.asarray(r[0,:]).squeeze().real[0:200]) # the asarray and squeeze are just annoyances we have to do because we came from a matrix
plt.plot(np.asarray(r[1,:]).squeeze().real[0:200])
plt.plot(np.asarray(r[2,:]).squeeze().real[0:200])
plt.show()
../_images/doa_time_domain.svg

Tenga en cuenta los cambios de fase entre elementos como esperamos que sucedan (a menos que la se帽al llegue de forma paralela, en cuyo caso llegar谩 a todos los elementos al mismo tiempo y no habr谩 un cambio, establezca theta en 0 para ver). El elemento 0 parece llegar primero, y los dem谩s se retrasan ligeramente. Intente ajustar el 谩ngulo y vea qu茅 sucede.

Como paso final, agreguemos ruido a esta se帽al recibida, ya que cada se帽al con la que trataremos tiene cierta cantidad de ruido. Queremos aplicar el ruido despu茅s de aplicar el factor de matriz, porque cada elemento experimenta una se帽al de ruido independiente (podemos hacer esto porque AWGN con un cambio de fase aplicado sigue siendo AWGN):

n = np.random.randn(Nr, N) + 1j*np.random.randn(Nr, N)
r = r + 0.5*n # r and n are both 3x10000
../_images/doa_time_domain_with_noise.svg

DOA convencional

Ahora procesaremos estos ejemplos r, fingiendo que no conocemos el 谩ngulo de llegada, y realizaremos DOA, que implica estimar los 谩ngulos de llegada con un DSP y algo de c贸digo Python. Como se analiz贸 anteriormente en este cap铆tulo, el beamforming y realizar DOA son muy similares y a menudo se basan en las mismas t茅cnicas. A lo largo del resto de este cap铆tulo investigaremos diferentes 鈥渂eamfirming鈥, y para cada uno comenzaremos con el c贸digo/matem谩tica del beamforming que calcula los pesos, w. Estos pesos se pueden 鈥渁plicar鈥 a la se帽al entrante r mediante la ecuaci贸n simple w^H r, o en Python w.conj().T @ r. En el ejemplo anterior, r es una matriz 3x10000, pero despu茅s de aplicar los pesos nos queda 1x10000, como si nuestro receptor solo tuviera una antena, y podemos utilizar un DSP RF normal para procesar la se帽al. Despu茅s de desarrollar el beamforming, lo aplicaremos al problema de DOA.

Comenzaremos con el enfoque del beamforming 鈥渃onvencional鈥, tambi茅n conocido como beamforming de retardo y suma. Nuestro vector de pesos w necesita ser una matriz 1D para una matriz lineal uniforme; en nuestro ejemplo de tres elementos, w es una matriz 3x1 de pesos complejos. Con la formaci贸n de haces convencional dejamos la magnitud de los pesos en 1 y ajustamos las fases para que la se帽al se sume constructivamente en la direcci贸n de nuestra se帽al deseada, a la que nos referiremos como \theta. 隆Resulta que estos son exactamente los mismos c谩lculos que hicimos arriba!

w_{conventional} = e^{-2j \pi d k \sin(\theta)}

o en Python:

w = np.exp(-2j * np.pi * d * np.arange(Nr) * np.sin(theta)) # Conventional, aka delay-and-sum, beamformer
r = w.conj().T @ r # example of applying the weights to the received signal (i.e., perform the beamforming)

donde Nr es el n煤mero de elementos en nuestra matriz lineal uniforme con espaciado de d fracciones de longitud de onda (m谩s a menudo ~0,5). Como puede ver, los pesos no dependen de nada m谩s que de la geometr铆a de la matriz y el 谩ngulo de inter茅s. Si nuestra matriz implicara calibrar la fase, tambi茅n incluir铆amos esos valores de calibraci贸n.

Pero 驴c贸mo sabemos el 谩ngulo de inter茅s theta? Debemos comenzar realizando DOA, que implica escanear (muestrear) todas las direcciones de llegada desde -蟺 a +蟺 (-180 a +180 grados), por ejemplo, en incrementos de 1 grado. En cada direcci贸n calculamos los pesos usando un formador de haces; Comenzaremos usando el formador de haz convencional. Aplicar los pesos a nuestra se帽al r nos dar谩 una matriz 1D de muestras, como si la recibi茅ramos con 1 antena direccional. Luego podemos calcular la potencia en la se帽al tomando la varianza con np.var(), y repetir para cada 谩ngulo en nuestro escaneo. Trazaremos los resultados y los veremos con nuestros ojos/cerebro humanos, pero lo que hace la mayor铆a de los DSP RF es encontrar el 谩ngulo de potencia m谩xima (con un algoritmo de b煤squeda de picos) y llamarlo estimaci贸n de DOA.

theta_scan = np.linspace(-1*np.pi, np.pi, 1000) # 1000 different thetas between -180 and +180 degrees
results = []
for theta_i in theta_scan:
   w = np.exp(-2j * np.pi * d * np.arange(Nr) * np.sin(theta_i)) # Conventional, aka delay-and-sum, beamformer
   r_weighted = w.conj().T @ r # apply our weights. remember r is 3x10000
   results.append(10*np.log10(np.var(r_weighted))) # power in signal, in dB so its easier to see small and large lobes at the same time
results -= np.max(results) # normalize

# print angle that gave us the max value
print(theta_scan[np.argmax(results)] * 180 / np.pi) # 19.99999999999998

plt.plot(theta_scan*180/np.pi, results) # lets plot angle in degrees
plt.xlabel("Theta [Degrees]")
plt.ylabel("DOA Metric")
plt.grid()
plt.show()
../_images/doa_conventional_beamformer.svg

隆Encontramos nuestra se帽al! Probablemente est茅 empezando a darse cuenta de d贸nde entra el t茅rmino matriz dirigida el茅ctricamente. Intente aumentar la cantidad de ruido para llevarlo al l铆mite; es posible que necesite simular la recepci贸n de m谩s muestras para SNR bajas. Intente tambi茅n cambiar la direcci贸n de llegada.

If you prefer viewing angle on a polar plot, use the following code:

fig, ax = plt.subplots(subplot_kw={'projection': 'polar'})
ax.plot(theta_scan, results) # MAKE SURE TO USE RADIAN FOR POLAR
ax.set_theta_zero_location('N') # make 0 degrees point up
ax.set_theta_direction(-1) # increase clockwise
ax.set_rlabel_position(55)  # Move grid labels away from other labels
plt.show()
Example polar plot of performing direction of arrival (DOA) showing the beam pattern and 180 degree ambiguity

Seguiremos viendo este patr贸n de bucles sobre 谩ngulos y teniendo alg煤n m茅todo para calcular los pesos del beamforming y luego aplicarlos a la se帽al recibida. En el pr贸ximo m茅todo de beamforming (MVDR), usaremos nuestra se帽al recibida r como parte de los c谩lculos de peso, lo que la convierte en una t茅cnica adaptativa. Pero primero investigaremos algunas cosas interesantes que suceden con los arreglos en fase, incluido por qu茅 tenemos ese segundo pico a 160 grados.

Ambiguedad de 180 grados

Hablemos de por qu茅 hay un segundo pico a 160 grados; el DOA que simulamos fue de 20 grados, pero no es casualidad que 180 - 20 = 160. Imag铆nese tres antenas omnidireccionales en l铆nea colocadas sobre una mesa. La mira del conjunto es de 90 grados con respecto al eje del conjunto, como se indica en el primer diagrama de este cap铆tulo. Ahora imagine el transmisor frente a las antenas, tambi茅n sobre la mesa (muy grande), de modo que su se帽al llegue en un 谩ngulo de +20 grados desde el punto de mira. Bueno, la matriz ve el mismo efecto ya sea que la se帽al llegue con respecto a su frente o atr谩s, el retraso de fase es el mismo, como se muestra a continuaci贸n con los elementos de la matriz en rojo y los dos posibles DOA del transmisor en verde. Por lo tanto, cuando realizamos el algoritmo DOA, siempre habr谩 una ambig眉edad de 180 grados como esta, la 煤nica forma de evitarla es tener una matriz 2D o una segunda matriz 1D colocada en cualquier otro 谩ngulo con respecto a la primera matriz. Quiz谩s se pregunte si esto significa que tambi茅n podr铆amos calcular solo -90 a +90 grados para ahorrar ciclos de c谩lculo, 隆y estar铆a en lo cierto!

../_images/doa_from_behind.svg

Orientaci贸n del arreglo

Para demostrar el siguiente concepto, intentemos barrer el 谩ngulo de llegada (AoA) de -90 a +90 grados en lugar de mantenerlo constante en 20:

Animation of direction of arrival (DOA) showing the broadside of the array

A medida que nos acercamos a la orientaci贸n de la matriz (tambi茅n conocido como endfire), que es cuando la se帽al llega al eje de la matriz o cerca de 茅l, el rendimiento disminuye. Vemos dos degradaciones principales: 1) el l贸bulo principal se ensancha y 2) obtenemos ambig眉edad y no sabemos si la se帽al proviene de la izquierda o de la derecha. Esta ambig眉edad se suma a la ambig眉edad de 180 grados discutida anteriormente, donde obtenemos un l贸bulo adicional en 180 theta, lo que hace que ciertos AoA conduzcan a tres l贸bulos de aproximadamente el mismo tama帽o. Sin embargo, esta ambig眉edad general tiene sentido, los cambios de fase que ocurren entre los elementos son id茅nticos ya sea que la se帽al llegue del lado izquierdo o derecho. el eje de la matriz. Al igual que con la ambig眉edad de 180 grados, la soluci贸n es utilizar una matriz 2D o dos matrices 1D en diferentes 谩ngulos. En general, la formaci贸n de haces funciona mejor cuando el 谩ngulo est谩 m谩s cerca del eje de punter铆a.

Cuando d no es 位/2

Hasta ahora hemos estado usando una distancia entre elementos, d, igual a la mitad de la longitud de onda. Entonces, por ejemplo, una matriz dise帽ada para WiFi de 2,4 GHz con espaciado 位/2 tendr铆a un espaciado de 3e8/2,4e9/2 = 12,5 cm o aproximadamente 5 pulgadas, lo que significa que una matriz de 4x4 elementos tendr铆a aproximadamente 15鈥 x 15鈥 x el Altura de las antenas. Hay ocasiones en las que es posible que una matriz no pueda lograr exactamente un espaciado 位/2, como cuando el espacio es restringido o cuando la misma matriz tiene que funcionar en una variedad de frecuencias portadoras.

Examinemos cu谩ndo el espaciado es mayor que 位/2, es decir, demasiado espaciado, variando d entre 位/2 y 4位. Eliminaremos la mitad inferior del gr谩fico polar ya que de todos modos es un espejo de la parte superior.

Animation of direction of arrival (DOA) showing what happens when distance d is much more than half-wavelength

Como puede ver, adem谩s de la ambig眉edad de 180 grados que discutimos anteriormente, ahora tenemos ambig眉edad adicional, y empeora a medida que d aumenta (forma de l贸bulos extra/incorrectos). Estos l贸bulos adicionales se conocen como l贸bulos de rejilla y son el resultado de un 鈥渁liasing espacial鈥. Como aprendimos en el cap铆tulo Muestreo IQ, cuando no tomamos muestras lo suficientemente r谩pido obtenemos alias. Lo mismo ocurre en el 谩mbito espacial; si nuestros elementos no est谩n lo suficientemente espaciados con respecto a la frecuencia portadora de la se帽al, obtenemos resultados basura en nuestro an谩lisis. 隆Puedes pensar en espaciar las antenas como espacio de muestreo! En este ejemplo podemos ver que los l贸bulos de la rejilla no se vuelven demasiado problem谩ticos hasta que d > 位, pero ocurrir谩n tan pronto como supere el espaciado 位/2.

Ahora bien, 驴qu茅 sucede cuando d es menor que 位/2, como cuando necesitamos ajustar la matriz en un espacio peque帽o? Repitamos la misma simulaci贸n:

Animation of direction of arrival (DOA) showing what happens when distance d is much less than half-wavelength

Si bien el l贸bulo principal se ensancha a medida que d disminuye, todav铆a tiene un m谩ximo a 20 grados y no hay l贸bulos reticulares, por lo que, en teor铆a, esto a煤n funcionar铆a (al menos con una SNR alta). Para comprender mejor qu茅 se rompe cuando d se vuelve demasiado peque帽o, repitamos el experimento pero con una se帽al adicional que llega desde -40 grados:

Animation of direction of arrival (DOA) showing what happens when distance d is much less than half-wavelength and there are two signals present

Una vez que llegamos a un nivel inferior a 位/4, no se puede distinguir entre los dos caminos diferentes y la matriz funciona mal. Como veremos m谩s adelante en este cap铆tulo, existen t茅cnicas de formaci贸n de haces que proporcionan haces m谩s precisos que la formaci贸n de haces convencional, pero mantener d lo m谩s cerca posible de 位/2 seguir谩 siendo un tema.

MVDR/Capon Beamformer

Ahora veremos un beamformer que es un poco m谩s complicado que la t茅cnica convencional/de retardo y suma, pero que tiende a funcionar mucho mejor, llamado Respuesta sin distorsi贸n de varianza m铆nima (MVDR) o Capon Beamformer. Recuerde que la varianza de una se帽al corresponde a la cantidad de potencia que hay en la se帽al. La idea detr谩s de MVDR es mantener la se帽al en el 谩ngulo de inter茅s con una ganancia fija de 1 (0 dB), minimizando al mismo tiempo la variaci贸n/potencia total de la se帽al formada en haz resultante. Si nuestra se帽al de inter茅s se mantiene fija, minimizar la potencia total significa minimizar las interferencias y el ruido tanto como sea posible. A menudo se lo denomina beamformer 鈥渆stad铆sticamente 贸ptimo鈥.

El MVDR/Capon beamformer puede resumirse en la siguiente ecuaci贸n:

w_{mvdr} = \frac{R^{-1} a}{a^H R^{-1} a}

donde R es la estimaci贸n de la matriz de covarianza basada en las muestras recibidas, calculada multiplicando r por la transpuesta compleja conjugada de s铆 mismo, es decir, R = r r^H, y el resultado ser谩 una matriz de tama帽o Nr x Nr (3x3 en los ejemplos que hemos visto hasta ahora). Esta matriz de covarianza nos dice qu茅 tan similares son las muestras recibidas de los tres elementos. El vector a es el vector de direcci贸n correspondiente a la direcci贸n deseada y se analiz贸 al principio de este cap铆tulo.

Si ya conocemos la direcci贸n de la se帽al de inter茅s, y esa direcci贸n no cambia, solo nos queda calcular los pesos una vez y simplemente usarlos para recibir nuestra se帽al de inter茅s. Aunque incluso si la direcci贸n no cambia, nos beneficia recalcular estos pesos peri贸dicamente, para tener en cuenta los cambios en la interferencia/ruido, raz贸n por la cual nos referimos a estos beamformer digitales no convencionales como beamforming 鈥渁daptativa鈥; utilizan informaci贸n de la se帽al que recibimos para calcular los mejores pesos. Solo como recordatorio, podemos realizar formaci贸n de haces usando MVDR calculando estos pesos y aplic谩ndolos a la se帽al con w.conj().T @ r, tal como lo hicimos en el m茅todo convencional, el 煤nico La diferencia es c贸mo se calculan los pesos.

Para realizar DOA utilizando el beamformer MVDR, simplemente repetimos el c谩lculo de MVDR mientras escaneamos todos los 谩ngulos de inter茅s. Es decir, actuamos como si nuestra se帽al viniera del 谩ngulo \theta, incluso si no es as铆. En cada 谩ngulo calculamos los pesos MVDR, luego los aplicamos a la se帽al recibida y luego calculamos la potencia en la se帽al. El 谩ngulo que nos da la mayor potencia es nuestra estimaci贸n de DOA, o mejor a煤n, podemos trazar la potencia en funci贸n del 谩ngulo para ver el patr贸n del haz, como hicimos anteriormente con el formador de haz convencional, de esa manera no necesitamos asumir cuantas se帽ales est谩n presentes.

En Python podemos implementar el beamformer MVDR/Capon de la siguiente manera, lo cual se har谩 como una funci贸n para que sea f谩cil de usar m谩s adelante:

# theta is the direction of interest, in radians, and r is our received signal
def w_mvdr(theta, r):
   a = np.exp(-2j * np.pi * d * np.arange(Nr) * np.sin(theta)) # steering vector in the desired direction theta
   a = a.reshape(-1,1) # make into a column vector (size 3x1)
   R = r @ r.conj().T # Calc covariance matrix. gives a Nr x Nr covariance matrix of the samples
   Rinv = np.linalg.pinv(R) # 3x3. pseudo-inverse tends to work better/faster than a true inverse
   w = (Rinv @ a)/(a.conj().T @ Rinv @ a) # MVDR/Capon equation! numerator is 3x3 * 3x1, denominator is 1x3 * 3x3 * 3x1, resulting in a 3x1 weights vector
   return w

Al utilizar este beamformer MVDR en el contexto de DOA, obtenemos el siguiente ejemplo de Python:

theta_scan = np.linspace(-1*np.pi, np.pi, 1000) # 1000 different thetas between -180 and +180 degrees
results = []
for theta_i in theta_scan:
   w = w_mvdr(theta_i, r) # 3x1
   r_weighted = w.conj().T @ r # apply weights
   power_dB = 10*np.log10(np.var(r_weighted)) # power in signal, in dB so its easier to see small and large lobes at the same time
   results.append(power_dB)
results -= np.max(results) # normalize

Cuando se aplica a la simulaci贸n del ejemplo DOA anterior, obtenemos lo siguiente:

../_images/doa_capons.svg

Parece funcionar bien, pero para compararlo realmente con otras t茅cnicas tendremos que crear un problema m谩s interesante. Configuremos una simulaci贸n con una matriz de 8 elementos que recibe tres se帽ales desde diferentes 谩ngulos: 20, 25 y 40 grados, con la de 40 grados recibida a una potencia mucho menor que las otras dos, como una forma de darle vida a las cosas. Nuestro objetivo ser谩 detectar las tres se帽ales, lo que significa que queremos poder ver picos notables (lo suficientemente altos como para que un algoritmo de b煤squeda de picos los extraiga). El c贸digo para generar este nuevo escenario es el siguiente:

Nr = 8 # 8 elements
theta1 = 20 / 180 * np.pi # convert to radians
theta2 = 25 / 180 * np.pi
theta3 = -40 / 180 * np.pi
a1 = np.exp(-2j * np.pi * d * np.arange(Nr) * np.sin(theta1)).reshape(-1,1) # 8x1
a2 = np.exp(-2j * np.pi * d * np.arange(Nr) * np.sin(theta2)).reshape(-1,1)
a3 = np.exp(-2j * np.pi * d * np.arange(Nr) * np.sin(theta3)).reshape(-1,1)
# we'll use 3 different frequencies.  1xN
tone1 = np.exp(2j*np.pi*0.01e6*t).reshape(1,-1)
tone2 = np.exp(2j*np.pi*0.02e6*t).reshape(1,-1)
tone3 = np.exp(2j*np.pi*0.03e6*t).reshape(1,-1)
r = a1 @ tone1 + a2 @ tone2 + 0.1 * a3 @ tone3
n = np.random.randn(Nr, N) + 1j*np.random.randn(Nr, N)
r = r + 0.05*n # 8xN

Puedes poner este c贸digo en la parte superior de tu script, ya que estamos generando una se帽al diferente a la del ejemplo original. Si ejecutamos nuestro beamformer MVDR en este nuevo escenario obtenemos los siguientes resultados:

../_images/doa_capons2.svg

Funciona bastante bien, podemos ver las dos se帽ales recibidas con solo 5 grados de diferencia, y tambi茅n podemos ver la tercera se帽al (a -40 o 320 grados) que se recibi贸 a una d茅cima parte de la potencia de las dem谩s. Ahora ejecutemos el formador de haces convencional en este nuevo escenario:

../_images/doa_complex_scenario.svg

Si bien puede tener una forma bonita, no encuentra las tres se帽ales en absoluto鈥 Al comparar estos dos resultados, podemos ver el beneficio de usar un beamformer m谩s complejo y 鈥渁daptativo鈥.

Como comentario breve para el lector interesado, en realidad existe una optimizaci贸n que se puede realizar al realizar DOA con MVDR, mediante un truco. Recuerde que calculamos la potencia de una se帽al tomando la varianza, que es la media de la magnitud al cuadrado (asumiendo que el valor promedio de nuestra se帽al es cero, lo que casi siempre es el caso de la RF de banda base). Podemos representar tomando la potencia en nuestra se帽al despu茅s de aplicar nuestros pesos como:

P_{mvdr} = \frac{1}{N} \sum_{n=0}^{N-1} \left| w^H_{mvdr} r_n \right|^2

Si reemplazamos la ecuaci贸n para los pesos MVDR obtenemos:

P_{mvdr} = \frac{1}{N} \sum_{n=0}^{N-1} \left| \left( \frac{R^{-1} a}{a^H R^{-1} a} \right)^H r_n \right|^2

  = \frac{1}{N} \sum_{n=0}^{N-1} \left| \frac{a^H R^{-1}}{a^H R^{-1} a} r_n \right|^2

 ... \mathrm{math}

  = \frac{1}{a^H R^{-1} a}

Lo que significa que no tenemos que aplicar los pesos en absoluto, esta ecuaci贸n final definida para la potencia se puede usar directamente en nuestro escaneo DOA, ahorr谩ndonos algunos c谩lculos:

def power_mvdr(theta, r):
    a = np.exp(-2j * np.pi * d * np.arange(r.shape[0]) * np.sin(theta)) # steering vector in the desired direction theta_i
    a = a.reshape(-1,1) # make into a column vector (size 3x1)
    R = r @ r.conj().T # Calc covariance matrix. gives a Nr x Nr covariance matrix of the samples
    Rinv = np.linalg.pinv(R) # 3x3. pseudo-inverse tends to work better than a true inverse
    return 1/(a.conj().T @ Rinv @ a).squeeze()

Para usar esto en la simulaci贸n anterior, dentro del bucle for, lo 煤nico que queda por hacer es tomar el 10*np.log10() y listo, no hay pesos que aplicar; 隆Nos saltamos el c谩lculo de los pesos!

Hay muchos m谩s formadores de haces por ah铆, pero a continuaci贸n nos tomaremos un momento para analizar c贸mo la cantidad de elementos afecta nuestra capacidad para realizar formaci贸n de haces y DOA.

Numero de elementos

Proximamente!

MUSIC

Ahora cambiaremos de tema y hablaremos de un tipo diferente de beamforming. Todos los anteriores han ca铆do en la categor铆a de 鈥渞etraso y suma鈥, pero ahora nos sumergiremos en los m茅todos 鈥渟ubespaciales鈥. Estos implican dividir el subespacio de se帽al y el subespacio de ruido, lo que significa que debemos estimar cu谩ntas se帽ales recibe la matriz para obtener un buen resultado. La clasificaci贸n de se帽ales m煤ltiples (MUSIC) es un m茅todo subespacial muy popular que implica calcular los vectores propios de la matriz de covarianza (que, por cierto, es una operaci贸n computacionalmente intensiva). Dividimos los vectores propios en dos grupos: subespacio de se帽al y subespacio de ruido, luego proyectamos vectores de direcci贸n en el subespacio de ruido y nos dirigimos hacia los nulos. Esto puede parecer confuso al principio, 隆y es parte del por qu茅 M脷SICA parece magia negra!

La ecuaci贸n central de MUSIC es la siguiente:

\hat{\theta} = \mathrm{argmax}\left(\frac{1}{a^H V_n V^H_n a}\right)

donde V_n es la lista de vectores propios subespaciales de ruido que mencionamos (una matriz 2D). Se encuentra calculando primero los vectores propios de R, lo cual se hace simplemente con w, v = np.linalg.eig(R) en Python, y luego dividiendo los vectores (:code :w) en funci贸n de cu谩ntas se帽ales creemos que est谩 recibiendo la matriz. Hay un truco para estimar el n煤mero de se帽ales del que hablaremos m谩s adelante, pero debe estar entre 1 y Nr - 1. Es decir, si est谩s dise帽ando un arreglo, al momento de elegir el n煤mero de elementos debes tener uno m谩s que el n煤mero de se帽ales anticipadas. Una cosa a tener en cuenta sobre la ecuaci贸n anterior es que V_n no depende del factor de matriz a, por lo que podemos precalcularlo antes de comenzar a recorrer theta. El c贸digo completo de MUSIC es el siguiente:

num_expected_signals = 3 # Try changing this!

# part that doesn't change with theta_i
R = r @ r.conj().T # Calc covariance matrix, it's Nr x Nr
w, v = np.linalg.eig(R) # eigenvalue decomposition, v[:,i] is the eigenvector corresponding to the eigenvalue w[i]
eig_val_order = np.argsort(np.abs(w)) # find order of magnitude of eigenvalues
v = v[:, eig_val_order] # sort eigenvectors using this order
# We make a new eigenvector matrix representing the "noise subspace", it's just the rest of the eigenvalues
V = np.zeros((Nr, Nr - num_expected_signals), dtype=np.complex64)
for i in range(Nr - num_expected_signals):
   V[:, i] = v[:, i]

theta_scan = np.linspace(-1*np.pi, np.pi, 1000) # -180 to +180 degrees
results = []
for theta_i in theta_scan:
    a = np.exp(-2j * np.pi * d * np.arange(Nr) * np.sin(theta_i)) # array factor
    a = a.reshape(-1,1)
    metric = 1 / (a.conj().T @ V @ V.conj().T @ a) # The main MUSIC equation
    metric = np.abs(metric.squeeze()) # take magnitude
    metric = 10*np.log10(metric) # convert to dB
    results.append(metric)

results /= np.max(results) # normalize

Al ejecutar este algoritmo en el complejo escenario que hemos estado usando, obtenemos los siguientes resultados muy precisos, que muestran el poder de MUSIC:

Example of direction of arrival (DOA) using MUSIC algorithm beamforming

Ahora, 驴qu茅 pasar铆a si no tuvi茅ramos idea de cu谩ntas se帽ales estaban presentes? Bueno, hay un truco; ordena las magnitudes de los valores propios de mayor a menor y graficalas (puede ser 煤til graficarlas en dB):

plot(10*np.log10(np.abs(w)),'.-')
../_images/doa_eigenvalues.svg

Los valores propios asociados con el subespacio de ruido ser谩n los m谩s peque帽os y todos tender谩n alrededor del mismo valor, por lo que podemos tratar estos valores bajos como un 鈥減iso de ruido鈥, y cualquier valor propio por encima del piso de ruido representa una se帽al. Aqu铆 podemos ver claramente que se est谩n recibiendo tres se帽ales y ajustar nuestro algoritmo de MUSIC en consecuencia. Si no tiene muchas muestras de IQ para procesar o las se帽ales tienen una SNR baja, es posible que la cantidad de se帽ales no sea tan obvia. Si茅ntase libre de jugar ajustando num_expected_signals entre 1 y 7; descubrir谩 que subestimar el n煤mero provocar谩 que falten se帽ales, mientras que sobreestimar solo afectar谩 ligeramente el rendimiento.

Otro experimento que vale la pena probar con MUSIC es ver qu茅 tan cerca pueden llegar dos se帽ales (en 谩ngulo) sin dejar de distinguirlas; Las t茅cnicas subespaciales son especialmente buenas en eso. La siguiente animaci贸n muestra un ejemplo, con una se帽al a 18 grados y otro 谩ngulo de llegada que se desplaza lentamente.

../_images/doa_music_animation.gif

ESPRIT

Proximamente!

Patr贸n de antena Inactivo

Recuerde que nuestro vector de direcci贸n lo seguimos viendo,

np.exp(-2j * np.pi * d * np.arange(Nr) * np.sin(theta))

encapsula la geometr铆a de la matriz y su 煤nico otro par谩metro es la direcci贸n hacia la que desea dirigirse. Podemos calcular y trazar el patr贸n de antena 鈥淚nactivo鈥 (respuesta del conjunto) cuando se dirige hacia una determinada direcci贸n, lo que nos indicar谩 la respuesta natural del conjunto si no realizamos ninguna formaci贸n de haz adicional. Esto se puede hacer tomando la FFT de los pesos conjugados complejos, no se necesita bucle for. La parte complicada es asignar los contenedores de la salida FFT a 谩ngulos en radianes o grados, lo que implica un arcoseno como se puede ver en el ejemplo completo a continuaci贸n:

N_fft = 512
theta = theta_degrees / 180 * np.pi # doesnt need to match SOI, we arent processing samples, this is just the direction we want to point at
w = np.exp(-2j * np.pi * d * np.arange(Nr) * np.sin(theta)) # steering vector
w = np.conj(w) # or else our answer will be negative/inverted
w_padded = np.concatenate((w, np.zeros(N_fft - Nr))) # zero pad to N_fft elements to get more resolution in the FFT
w_fft_dB = 10*np.log10(np.abs(np.fft.fftshift(np.fft.fft(w_padded)))**2) # magnitude of fft in dB
w_fft_dB -= np.max(w_fft_dB) # normalize to 0 dB at peak

# Map the FFT bins to angles in radians
theta_bins = np.arcsin(np.linspace(-1, 1, N_fft)) # in radians

# find max so we can add it to plot
theta_max = theta_bins[np.argmax(w_fft_dB)]

fig, ax = plt.subplots(subplot_kw={'projection': 'polar'})
ax.plot(theta_bins, w_fft_dB) # MAKE SURE TO USE RADIAN FOR POLAR
ax.plot([theta_max], [np.max(w_fft_dB)],'ro')
ax.text(theta_max - 0.1, np.max(w_fft_dB) - 4, np.round(theta_max * 180 / np.pi))
ax.set_theta_zero_location('N') # make 0 degrees point up
ax.set_theta_direction(-1) # increase clockwise
ax.set_rlabel_position(55)  # Move grid labels away from other labels
ax.set_thetamin(-90) # only show top half
ax.set_thetamax(90)
ax.set_ylim([-30, 1]) # because there's no noise, only go down 30 dB
plt.show()
../_images/doa_quiescent.svg

Resulta que este patr贸n coincidir谩 casi exactamente con el patr贸n que se obtiene al realizar DOA con el beamforming convencional (retraso y suma), cuando hay un solo tono presente en 鈥渢heta_grados鈥 y poco o ning煤n ruido. La gr谩fica puede verse diferente debido a qu茅 tan bajo llega el eje y en dB, o debido al tama帽o de la FFT utilizada para crear este patr贸n de respuesta inactivo. Intente modificar theta_ Degrees o el n煤mero de elementos Nr para ver c贸mo cambia la respuesta.

2D DOA

Proximamente!

Nulos de direcci贸n

Proximament!

Conclusi贸n y Referencias

Todo el c贸digo Python, incluido el c贸digo utilizado para generar las figuras/animaciones, se puede encontrar en la p谩gina de GitHub del libro.