using DSP
using FFTW
using Gadfly
using LinearAlgebra
using Statistics
using ColorSchemes
using Printf
function plot_signals(X::Matrix{ComplexF64}, n::Int64, fs::Float64)
t = 1/fs:1/fs:n * fs
n_elements = size(X)[1]
p_i = plot(Guide.title("In Phase"), Guide.ylabel("Amplitude"), Guide.xlabel("Time [s]"))
p_q = plot(Guide.title("Quadrature"), Guide.ylabel("Amplitude"), Guide.xlabel("Time [s]"))
#colors = [colorant"blue", colorant"red", colorant"green", colorant"cyan", colorant"yellow"]
colors = ColorSchemes.tab10.colors
for i in 1:n_elements
push!(p_i, layer(x=t[1:n], y=real(X[i, 1:n]), Geom.line, color=[colors[i]]))
push!(p_q, layer(x=t[1:n], y=imag(X[i, 1:n]), Geom.line, color=[colors[i]]))
end
stack = vstack(p_i, p_q)
end;
Direction Finding with MUSIC
This post simulates direction finding with a 4 element uniform linear array (ULA) using the MUSIC algorithm.
There are three steps:
- Generate a signal for each antenna with the appropriate phase delay for its angle of arrivale (AoA)
- Apply the MUSIC algorithm
- Evaluate the performance of the MUSIC algorithm
Step 1 - Generate Simulated Signals
For this simulation, there are two parts to generating the signals that are needed.
The first is generating four signals that are identical other than their phase. The phase for each signal gets delayed according to the intended angle of arrival (AoA) and the spacing of the array.
The second step is simply adding additive white Gaussian noise (AWGN).
power(X) = 10 .* log10.(mean(abs.(X) .^ 2, dims=2))
function add_noise!(noise_power::Float64, X::Matrix{ComplexF64})
noise = sqrt((10 .^ (noise_power./10))) * randn(ComplexF64, size(X))
X + noise
end;
function add_noise(power::Float64, X::Matrix{ComplexF64})
noise = sqrt((10 .^ (noise_power./10))) * randn(ComplexF64, size(X))
X + noise
end;
function gen_vectors(ϕ::Float64, fs::Float64, f::Float64, duration::Float64, n_elements::Int64)::Matrix{ComplexF64}
t = 1/fs:1/fs:duration
X = zeros(ComplexF64, n_elements, length(t))
Φ = ϕ .* collect(0:1:n_elements-1)
for i in 1:n_elements
temp = im * sin.(2π * f .* t .+ Φ[i]) + cos.(2π * f .* t .+ Φ[i]);
X[i, :] = temp
end
X;
end;
# Generate the four signals as if they are coming from an emitter that's at 34°
θ = 34.0 # true angle, degrees
ϕ = π * cos(deg2rad(θ)) # phase difference, radians
fs = 30e3 # sampling frquency, Hz
f = 4e3 # tone frequency, Hz
duration = 1.0 # signal duration, s
c = 299792458.0 # speed of light, m/s
λ = c / f; # wavelength, m
n_elements = 4; # number of elements
X = gen_vectors(ϕ, fs, f, duration, n_elements);
plot_signals(X, 50, fs)
<?xml version=”1.0” encoding=”UTF-8”?>
noise_power = -3.0 # dB, W
noise_power_linear = 10 .^ (noise_power ./ 10) # W
snr = 3 # SNR, log-scale, unitless
snr_linear = 10 .^ (snr ./ 10) # SNR, linear-scale, unitless
X_pwr = power(X) # Signal power, dB
X_pwr_linear = 10 .^ (X_pwr ./ 10) # Signal power, linear
noise = sqrt.(X_pwr_linear ./ snr_linear) .* randn(ComplexF64, size(X))
Xn = X + noise
n_pwr = power(noise)
Xn_pwr = power(Xn)
snr = X_pwr .- n_pwr
@printf "Avg power of X = [%2.2f %2.2f %2.2f %2.2f] dB\n" X_pwr[1] X_pwr[2] X_pwr[3] X_pwr[4]
@printf "Avg power of n = [%2.2f %2.2f %2.2f %2.2f] dB\n" n_pwr[1] n_pwr[2] n_pwr[3] n_pwr[4]
@printf "Avg power of Xn = [%2.2f %2.2f %2.2f %2.2f] dB\n" Xn_pwr[1] Xn_pwr[2] Xn_pwr[3] Xn_pwr[4]
@printf "SNR = [%2.2f %2.2f %2.2f %2.2f] dB\n" snr[1] snr[2] snr[3] snr[4]
plot_signals(Xn, 50, fs)
Avg power of X = [0.00 0.00 0.00 0.00] dB
Avg power of n = [-3.00 -3.03 -3.03 -2.99] dB
Avg power of Xn = [1.75 1.77 1.74 1.77] dB
SNR = [3.00 3.03 3.03 2.99] dB
<?xml version=”1.0” encoding=”UTF-8”?>
Step 2 - Apply the MUSIC Algorithm
There are three main components to the MUSIC algorith: Compute the covariance matrix, obtain the eigenvalues and eigenvectors, apply the eigenvector with the largest eigenvalue to the array manifold.
The first component is computing the covariance matrix. In this example there are four signals of some length, N. When the covariance matrix is computed we are left with a 4x4 symmetric matrix that shows the covariance between each the signal that each antenna receives.
Next, obtain the eigenvalues and eigenvector pairs. The key to MUSIC is separating the signal subspace and the noise subspace by discarding the eigenvector with the smallest eigenvalue (as this eigenvector corresponds to the noise subspace).
Finally, take the remaining eigenvectors and create the MUSIC pseudospectrum by applying the array manifold to the eigenvectors. The array manifold is what the array expects to receive for a given angle of arrival under ideal circumstances.
function array_manifold_vector(array::Vector{Float64}, θ::Float64)::Vector{ComplexF64}
exp.(-2π * im * cos(θ) .* array)
end;
function MUSIC(X::Matrix{ComplexF64}, n_targets::Int64, normalized_spacing::Float64)::Vector{Float64}
n_elements = size(X)[1]
element_locations = normalized_spacing .* 0.5 .* (n_elements .- 1 .- 2 .* collect(0:(n_elements-1)))
pseudo_spectrum_length = 1024;
vii = zeros(ComplexF64, n_elements, pseudo_spectrum_length);
thetas = range(0, stop=π, length=pseudo_spectrum_length)
for i in 1:pseudo_spectrum_length
vii[:, i] = array_manifold_vector(element_locations, thetas[i])
end
# Get the covariance matrix, R
R = Statistics.cor(X, X, dims=2);
# Get the eigenvalues and vectors of R
eigen_values = eigvals(R);
eigen_vectors = eigvecs(R);
# Sort the eigenvalues
i = sortperm(abs.(eigen_values));
eigen_vectors = eigen_vectors[:, i]
# Use the largest n_elements - 1 eigenvectors to make U_n
U_N = eigen_vectors[:, 1:n_elements-1];
U_N_sq = U_N * U_N';
# Generate the pseudospectrum
pseudo_spectrum = zeros(Float64, 1, pseudo_spectrum_length)
q = 0.0
for i in 1:pseudo_spectrum_length
q = vii'[[i], :] * U_N_sq * vii[:, [i]]
pseudo_spectrum[i] = abs(1.0/real(q[1]))
end
# Normalize the pseudospectrum to 0dB
pseudo_spectrum = 10 .* log10.(pseudo_spectrum/maximum(pseudo_spectrum));
# Return the pseudospectrum
vec(pseudo_spectrum)
end;
pseudo_spectrum = MUSIC(X, 1, 0.5)
(m, i) = findmax(pseudo_spectrum);
θ̂ = rad2deg(i/1024.0 * π);
error = abs(θ-θ̂)
@printf "θ = %2.2f°, θ̂ = %2.2f°, error = %2.4f°\n" θ θ̂ error
θ = 34.00°, θ̂ = 34.10°, error = 0.1016°
plot(x=range(0, stop=180, length=length(pseudo_spectrum)),
y=pseudo_spectrum,
Geom.line,
Coord.cartesian(xmin=0, xmax=180),
Guide.title("MUSIC Pseudospectrum"),
Guide.xlabel("Angle of Arrival (°)"),
Guide.ylabel("MUSIC Response (dB)"))
Step 3 - MUSIC Evaluation
The cell below sweeps through every angle in 1 degree increments and records the AoA where the array and MUSIC perform the worst (\Delta max). This will likely be around 1 degree or 179 degrees due to the grating lobes from the ULA.
emitter_angles = 1.0:1.0:179.0
Δ = zeros(size(emitter_angles))
n_elements = 4
for (index, θ) in enumerate(emitter_angles)
ϕ = π * cos(deg2rad(θ))
X = gen_vectors(ϕ, fs, f, duration, n_elements)
noise = noise = sqrt((10 .^ (noise_power./10))) * randn(ComplexF64, size(X))
pseudo_spectrum = MUSIC(X, 1, 0.5)
(m, i) = findmax(pseudo_spectrum);
θ̂ = rad2deg(i/1024.0 * π);
Δ[index] = abs(θ - θ̂)
end
(Δmax, Δmax_ind) = findmax(Δ)
Δmax_angle = emitter_angles[Δmax_ind]
println("Δmax = $Δmax at θ = $Δmax_angle")
Δmax = 0.25 at θ = 11.0
plot(x=emitter_angles, y=Δ,
Coord.cartesian(xmin=0, xmax=180),
Geom.line,
Guide.title("MUSIC Error as a function of AoA"),
Guide.xlabel("Angle of Arrival (°)"),
Guide.ylabel("Error (°)"))
