MUSIC Algorithm for Direction Finding

[ecosierra51]

weetabixharry:
I have read many articles about MUSIC, but I can´t to work alone.
I began to undestand with your example. Thank you, Very much!
I will begin to work with your second example.
I am working now with 3 channel interferometer. Antennae are placed in right angle, number 0 to North, number 1 as reference, and number 2 to the East. I measure amplitud and difference of fase of signal from 0, respected to 1; and amplitud and difference of fase of signal from 2, respected to 1.
When I finished to study your second example, I will submit the code changed to orthogonal placed antennae, in order receive your advice, for apply MUSIC algorithm.
Thanks, again!

ecosierra51

 

[weetabixharry]

I am working now with 3 channel interferometer. Antennae are placed in right angle, number 0 to North, number 1 as reference, and number 2 to the East. I measure amplitud and difference of fase of signal from 0, respected to 1; and amplitud and difference of fase of signal from 2, respected to 1.

Here are a couple of points to consider:

1) Please note that the MUSIC algorithm is very sensitive to calibration errors. You must know the locations of the antennas precisely (and you must model their gains and phases exactly). It is very easy to make my simple Matlab example work with synthetic signals, but much, much (much) more difficult using real-world equipment. Before trying MUSIC, you should make sure you can get reasonable results using more simple methods (such as conventional beamforming, or a Capon beamformer).

2) The number of signals that the basic MUSIC algorithm can handle is limited by the number of antennas. If you have N antennas, then you can handle at most N-1 signals (since the noise subspace must have at least 1 dimension). However, if I understand correctly, your interferometer setup receives a different polarisation at each antenna. I don't have much experience with polarised signals, but I believe this can help. I believe the MUSIC algorithm has been adapted specifically for polarised signals, so you may want to look into that.

 

[ecosierra51]

weetabixharry:
Thanks for your attention!

I have worked with your 2D example very well. Also, I modified array configuration, to 3 antennae orthogonally placed, as:
M = 2
N = 3
rx = [0, 0 ,1];
ry = [1, 0 ,0];
rz = [0, 0 ,0];
r = [rx; ry; rz];
Also, I tried with contour plot, in order to know Az, El coordinates. I don´t know to solve analitically.
Your 2D example works well with orthogonal array.
About your points to consider:
1) I am working with a self made interferometer in the HF band. It has precision of location of the antennas for 1 cm, (10 meter separation to reference antenna). Gain and phases are calibrated to 99%, 0.1º for every measurement.
I got reasonable results with clasic interferometer, correlative interferometer and conventional beamforming.
In HF band I have another problem: Signals can arrive from two DOAs: ground wave and ionosferic reflected wave.
I read about MUSIC don´t work with correlated signals.

Please, advise me about this problem.

2) This interferometer setup use 3 identical vertical polarized antennas. It measures phase diference, using one as reference, with 3 sinchronous calibrated receivers.
I will test your 2D example, ortogonal, with real signals next week.

Thanks in advance!
ecosierra51

 

[weetabixharry]

I modified array configuration, to 3 antennae orthogonally placed, as:
M = 2
N = 3
rx = [0, 0 ,1];
ry = [1, 0 ,0];
rz = [0, 0 ,0];
r = [rx; ry; rz];

With azimuth and elevation defined as in my example, rx, ry and rz must be column vectors (not row vectors). Therefore, r =[rx, ry, rz].

Signals can arrive from two DOAs: ground wave and ionosferic reflected wave.
I read about MUSIC don´t work with correlated signals.

Please, advise me about this problem.

You are correct - standard MUSIC generally does not work with fully correlated signals (i.e. where the magnitude of the complex correlation coefficient is equal to 1). The reason for this is that fully correlated signals together span only 1 dimension of the complex observation space. There are a few methods to solve this problem, but they may or may not be suitable for you:

1) Conventional spatial smoothing. Many people think spatial smoothing can only be applied to uniform linear arrays. This is not true. Spatial smoothing can be applied to any array which has two or more identical subarrays (i.e. they shift onto one another without rotation). These subarrays can be overlapping. For example, if you had a (4 x 4) uniform grid array like this:

x x x x
x x x x
x x x x
x x x x

Then you could view this as four identical overlapping (3 x 3) grid arrays:

x x x .      . x x x      . . . .      . . . .
x x x .      . x x x      x x x .      . x x x
x x x .      . x x x      x x x .      . x x x
. . . .      . . . .      x x x .      . x x x

Then, if you take the average of the covariance matrices of these subarrays, the resulting signal correlation can generally be seen to decrease. I can describe the maths behind this in more detail if you want. Of course, one disadvantage of this approach is that the resulting array is smaller (in my example, 9 antennas instead of 16). You always lose at least 1 antenna with conventional spatial smoothing.

2) Backward spatial smoothing. Again, many people think backward spatial smoothing can only be applied to linear arrays and, again, this is not true. Backward spatial smoothing is somewhat similar to conventional spatial smoothing, but instead of needing identical shifted subarrays, you need an array (or subarrays) that are a reflection of themself through the origin. In other words, it must be possible to transform the array manifold vector into its complex conjugate by simply reordering its elements. For example, the (4 x 4) grid array above has this property. Or, to keep the numbers simple, imagine a (2 x 2) grid array with:

[rx, ry, rz] =    -0.5   -0.5         0
                   0.5   -0.5         0
                  -0.5    0.5         0
                   0.5    0.5         0

If we write the locations in reverse order, we get:

reversed =         0.5    0.5         0
                  -0.5    0.5         0
                   0.5   -0.5         0
                  -0.5   -0.5         0

which is equal to -[rx,ry,rz]. Therefore, if we average the two resulting covariance matrices, the imaginary component of the correlation coefficient is eliminated. The array does not get smaller, but only the imaginary component is removed.

Of course, if your array has identical subarrays which reflect through their own centroids, it is possible to apply both conventional and backward spatial smoothing ("forward-backward" spatial smoothing).

3) Other signal diversity. If possible, you should try to increase the dimension of the observation space by exploiting other signal diversity. For example, if the incoming signals have different Doppler shifts/carrier offsets, then this can be exploited. Alternatively, if the signals arrive with different time delays (and the nature of the signal waveform is known in advance) then this can also be exploited. These methods require manipulation of the received data and corresponding modification to the received signal model. It is then sometimes possible to apply similar "smoothing" methods to the modified signal model, but this would take more time to describe.

4) Other methods (not MUSIC). If you can't get a big array with nice subarray symmetry and there is no other signal diversity, then a subspace-based algorithm like MUSIC may not be suitable. Try another approach.

 

[ecosierra51]

weetabixharry:
Thanks for your help and guides!
1) I solved problem about mistake in rows and columns as:
rx = [0; 0; 1];
ry = [1; 0; 0];
rz = [0; 0; 0];
r = [rx, ry, rz];

2) About methods to solve correlated signals:

Conventional spatial smoothing.
At this moment I can´t lose at least 1 antenna, because I only have 3 channels ready.

Other signal diversity.
It is very common that signals arrive with different time delays. For example, signal 1 propagates as ground wave, for 20 km; and signal 2 propagates as ionospheric wave, reflected in F2 layer at 300 km height.These rays have a difference near of 580 km. (2 miliseconds).
Please, explain how to manipulate the received data and received signal model, in order to apply “smoothing” methods.

3) How to analytically solve az/el for 2 peaks?

4) How to modify code to search Az for 360 degrees, in correspondence with orthogonal array?

Thanks in advance!
ecosierra51

 

[hrpreet]

hi weetabixharry!
i want to find the BER for the applied music algorithm. but the music algorithm has the coded sequence and hence output bits are more.. please help me with this

 

[weetabixharry]

It is very common that signals arrive with different time delays. For example, signal 1 propagates as ground wave, for 20 km; and signal 2 propagates as ionospheric wave, reflected in F2 layer at 300 km height.These rays have a difference near of 580 km. (2 miliseconds).
Please, explain how to manipulate the received data and received signal model, in order to apply “smoothing” methods.

In this case, it is not always best to approach the problem from an array processing perspective (since anything other than direction finding is essentially conventional signal processing and can also be solved using conventional methods). If the signal has good autocorrelation properties, then the delayed signal will have low cross-correlation with the direct signal. This means the two are fundamentally relatively easy to separate.

From an array processing perspective, you can treat this (after correctly rearranging the discrete signal data) as two signals that exist in orthogonal (or "almost" orthogonal) subspaces, making them easily separable using a subspace projection. From a more general perspective, you can separate the signals easily using matched filtering.

However, if your signals do not have good (delta-like) autocorrelation properties, then it is fundamentally more difficult to separate the direct signal from the delayed version. Essentially, this means you are pretty screwed, unless you can find some other way of suppressing the reflected component.

How to analytically solve az/el for 2 peaks?

I don't really understand this question. If you are looking for a closed-form solution, you may want to read about "root MUSIC". This was originally only formulated for linear arrays, but has been extended to arbitrary array geometries (at the expense of much greater computational complexity). See, for example: this paper and this paper.

4) How to modify code to search Az for 360 degrees, in correspondence with orthogonal array?

AzSearch = (0:1:360).';

 

[LoSfrek]

Hello everyone,
first of all, thank you Weetabixharry for your helpful code, very easy to read and understand.
I am reading a lot about MUSIC and I am trying to understand if it could be successfully used to estimate the DOA with a conformal, hexagonal array. The elements are directive and outward-facing, so the signal is not received from every element: the array must be divided into subarrays, and start from there to get it to work.
I have found some papers about it (for example here)
The hard part, though, is that I should then try to use it with an element spacing larger than lambda/2, and try to remove the subsequent ambiguity in DOAs. Do you think it would be too difficult to implement? Maybe I should think about some other algorithms?

 

[LeeYoung]

So, we don't need to change the code much. The most important thing is that we cannot use a linear array to estimate elevation; we must have a 2D or 3D array (this is obvious if you think about the rotational symmetry of a linear array). Therefore, I changed the array to a uniform circular array.

Then, we just have to calculate a separate 1D azimuth spectrum for each elevation. Here is the modified code:

---

Code Matlab M - [expand]
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
close all; clear all; clc;
 
 
% ======= (1) TRANSMITTED SIGNALS ======= %
 
% Signal source directions
az = [35;39;127]; % Azimuths
el = [63;14;57]; % Elevations
M = length(az); % Number of sources
 
% Transmitted signals
L = 200; % Number of data snapshots recorded by receiver
m = randn(M,L); % Example: normally distributed random signals
 
% ========= (2) RECEIVED SIGNAL ========= %
 
% Wavenumber vectors (in units of wavelength/2)
k = pi*[cosd(az).*cosd(el), sind(az).*cosd(el), sind(el)].';
 
% Number of antennas
N = 10; 
 
% Array geometry [rx,ry,rz] (example: uniform circular array)
radius = 0.5/sind(180/N);
rx = radius*cosd(360*(0:N-1).'/N);
ry = radius*sind(360*(0:N-1).'/N);
r = [rx, ry, zeros(N,1)];
 
% Matrix of array response vectors
A = exp(-1j*r*k);
 
% Additive noise
sigma2 = 0.01; % Noise variance
n = sqrt(sigma2)*(randn(N,L) + 1j*randn(N,L))/sqrt(2);
 
% Received signal
x = A*m + n;
 
% ========= (3) MUSIC ALGORITHM ========= %
 
% Sample covariance matrix
Rxx = x*x'/L;
 
% Eigendecompose
[E,D] = eig(Rxx);
[lambda,idx] = sort(diag(D)); % Vector of sorted eigenvalues
E = E(:,idx); % Sort eigenvalues accordingly
En = E(:,1:end-M); % Noise eigenvectors (ASSUMPTION: M IS KNOWN)
 
% MUSIC search directions
AzSearch = (0:1:180).'; % Azimuth values to search
ElSearch = (0:1:90); % Elevation values to search
 
% 2D MUSIC spectrum
Z = zeros(length(AzSearch),length(ElSearch));
for i = 1:length(ElSearch)
    % Elevation search value
    el = ElSearch(i);
    
    % Points on azimuth array manifold curve to search (for this el)
    kSearch = pi*[cosd(AzSearch)*cosd(el), ...
              sind(AzSearch)*cosd(el), ...
              ones(size(AzSearch))*sind(el)].';
    ASearch = exp(-1j*r*kSearch);
    
    % Compute azimuth spectrum for this elevation
    Z(:,i) = sum(abs(ASearch'*En).^2,2);
end
 
% Plot
figure();
surf(AzSearch, ElSearch, -10*log10(Z.'/N));
shading interp;
title('2D MUSIC Example');
xlabel('Azimuth (degrees)');
ylabel('Elevation (degrees)');
zlabel('MUSIC spectrum (dB)');
grid on; axis tight;

---

This produces a 2D spectrum like this:

View attachment 114265

Note that, for convenient viewing, I plot the negated spectrum. This means that we get peaks at the target directions (instead of nulls), which I find easier to view for the 2D case.

Of course, we don't have to stop at 2D. We can extend, for example, to 3D and 4D to estimate Doppler and time of arrival. This is interesting to experiment with in theory, but in practice a 3D or 4D exhaustive search is far too slow to compute. Instead, we normally find clever ways of reducing the problem (e.g. a 1D search for time of arrival, followed by a 2D azimuth-elevation search is much faster than an exhaustive 3D search).

hi, weetabixharry
I'm using the Linux 802.11 CSI tool to do the indoor positioning based WiFi.
I can obtain the csi trace, i.e, the received signal x from the userspace of Ubuntu by the modified firmware.
How can I do the estimation of Az and ToF by utilizing the csi matrix x? Of course I know MUSIC can do this, but I'm confused about how to construct the steering vector.
Could you please modify the matlab code for me?
Thank you in advance!