Direction of Arrival Estimation Using the MUSIC Algorithm and Analog Devices FMCOMMS5

This example shows how to estimate the direction of arrival of a signal by using the Communications Toolbox™ Support Package for Xilinx® Zynq®-Based Radio and an Analog Devices™ FMCOMMS5 radio. The example transmits a 10 kHz signal from a suitable transmitter and calculates the direction of arrival of the signal by using a functionality from the Phased Array System Toolbox™.

Hardware requirements for running this example:

1. ZC706 and FMCOMMS5 radio hardware for the receive side.

2. Another radio hardware for the transmit signal source. By default, the example uses an Analog Devices ADALM-Pluto radio. You can also use any supported Xilinx Zynq-based radio hardware from this support package. For more information, see Hardware and Software Requirements.

3. Five antennas and suitable cables to connect to the hardware.

4. A 1-to-4 power splitter, or three 1-to-2 power splitters.

5. A means to set up the four receive antennas in a fixed linear array, separated by the correct distance.

Introduction

This example demonstrates how to find the direction of arrival of a signal using a phased.MUSICEstimator System object from Phased Array System Toolbox and a ZC706 and FMCOMMS5 radio hardware. The process consists of two steps.

1.Phase calibration between the four FMCOMMS5 radio channels: The signal transmit source is connected to all four receiving channels. This connection allows you to measure and to correct the phase offset between the channels.

2. Estimation of direction of arrival: The four receive channels are connected to the four channel antenna array and direction of arrival estimation is estimated. To find the direction of arrival of a signal, the example uses a phased.MUSICEstimator (Phased Array System Toolbox) System object. The phased.MUSICEstimator System object implements the incoming signal's spatial spectrum by using a narrowband MUSIC beamformer for a uniform linear array (ULA).

Setup

This example requires Phased Array System Toolbox. Before running the example, perform these steps:

  1. Configure your host computer to work with the Xilinx Zynq-based radio support package. For more information, see Install Support Package for Xilinx Zynq-Based Radio.

  2. If you are using a Pluto radio as the transmit signal source, you must also configure your host computer to work with the ADALM-Pluto radio support package. For more information, see Install Support Package for Analog Devices ADALM-PLUTO Radio (Communications Toolbox Support Package for Analog Devices ADALM-Pluto Radio).

  3. Make sure that you have both the transmitter script zynqRadioDOAEstimationTransmitterExample and the receiver script zynqRadioDOAEstimationFMCOMMS5ML open, with each configured to run on its own SDR hardware in its own instance of MATLAB.

  4. For calibration at the receiver side, connect the power splitter from the transmitter to the four receiver channels of FMCOMMS5.

  5. Set up a linear antenna array of four antennas with adjacent antennas at a distance of lambda/2. Since we are using a center frequency of 2.4 GHz, a spacing of 6.25 cm is required between two adjacent antennas.

Running the Example

First, start the transmitter by running zynqRadioDOAEstimationTransmitterExample. You can run the rest of this example by executing the zynqRadioDOAEstimationFMCOMMS5ML script.

Capture RF Signal

To capture the RF tone signal into MATLAB, create an SDR receiver System object and configure the object to receive samples at the baseband rate.

DOARx.SDRDeviceName = 'FMCOMMS5';
DOARx.IPAddress = '192.168.3.2';

% Radio parameters

RadioBasebandRate = 1e6;
CenterFrequency = 2.4e9;
ToneFrequency = 10e3;
RadioFrameLength = 4000;

% Angles by which each channel needs to be rotated
angle1 = 0;
angle2 = 0;
angle3 = 0;

sdrReceiver = sdrrx( DOARx.SDRDeviceName, ...
    'IPAddress',        DOARx.IPAddress, ...
    'CenterFrequency',  CenterFrequency, ...
    'BasebandSampleRate', RadioBasebandRate,...
    'GainSource', 'Manual', ...
    'Gain', 20, ...
    'SamplesPerFrame', RadioFrameLength, ...
    'ChannelMapping',[1 2 3 4],...
    'OutputDataType',   'double', ...
    'ShowAdvancedProperties', true, ...
    'BypassUserLogic', true);

To visualize the received signal in frequency and time domain, use the Spectrum Analyzer and the Time Scope System objects.

spectrumScope = dsp.SpectrumAnalyzer('SampleRate', RadioBasebandRate);
timeScope = timescope('TimeSpanSource','property','TimeSpan',4/ToneFrequency,'SampleRate',RadioBasebandRate);

To find the direction of arrival of a signal, use the MUSICEstimator System object.

lambda = physconst('LightSpeed')/CenterFrequency;
array = phased.ULA('NumElements',4,'ElementSpacing',lambda/2);
estimator = phased.MUSICEstimator('SensorArray',array,...
             'OperatingFrequency',CenterFrequency,'ForwardBackwardAveraging',true,...
             'DOAOutputPort',true,'NumSignalsSource','Property',...
             'NumSignals',1,'SpatialSmoothing',1);

Set the simulation to capture 1000 seconds of data.

timeCounter = 0; % seconds
CalibrationTime = 1;  % Time for which Phase Calibration is done
RadioFrameTime  = (RadioFrameLength / RadioBasebandRate); % seconds

Phase Calibration Using Power Splitter

Estimate and correct the phase offset between all the receiver channels. To phase sync all signals, perform IQ rotation. The time scope shows the corresponding phase synced signals after IQ rotation.

try
    % Loop until the timeCounter reaches the target calibration time.
    while timeCounter < CalibrationTime

        [data, valid, ~] = sdrReceiver();

        if valid

                data = zynqRadioDOAEstimationAmplitudeNormalization(data);

                % Calculating the phase offset between 2 channels
                [angle1,angle2,angle3] = zynqRadioDOAEstimationFindPhaseDifference4Channels(data);

                % IQ Rotation to reduce the phase offset between all channels to 0
                data = zynqRadioDOAEstimationRotateIQ4Channels(data,angle1,angle2,angle3);

                % Visualize in time domain
                timeScope([real(data), imag(data)]);

                % Visualize frequency spectrum
                spectrumScope(data);

                % Set the limits in scopes
                dataMaxLimit = max(max(abs([real(data); imag(data)])));
                timeScope.YLimits = [ -dataMaxLimit*2, dataMaxLimit*2 ];
                timeCounter = timeCounter + RadioFrameTime;
        end
    end
catch ME
    rethrow(ME);
end

Change the Connections

  1. Carefully disconnect the power splitter from the transmitter and the receiver without powering off the board, and connect the FOUR channels of the FMCOMMS5 to the antenna array using equal length SMA cables.

  2. Connect the Pluto transmitter port to an antenna.

You can now move the transmitting antenna and see the corresponding change in arrival angle on the semicircular gauge. You may need to increase the gain on the transmitter object. You can do this by releasing the transmit object, modifying the script where necessary, and by running the script again.

disp("Please Remove the Power splitter and connect to the antenna array for DOA Estimation");

Display Direction of Arrival

Use the attached app to display the direction of arrival of the transmitted signal and to display the time and frequency plots of the received signal. Use the Time Scope to make sure that the signal strength is sufficient after connecting to the antenna array. Use the Spectrum Scope to monitor the power of the signal and the frequency of the signal received by the antenna array. Click the stop button to stop estimation of direction of arrival of the signal.

app = [];
if isempty(app)
   app = zynqRadioDOAEstimationApp();
end
StopTime = 1000; % Time to stop the simulation in seconds

Estimation of Angle of Arrival

Connect the [1 2 3 4] receiver channels of FMCOMMS5 to the antennas in the antenna array from left to right.

try
    % Loop until the timeCounter reaches the target stop time.

    while (timeCounter < StopTime) && (app.StopEstimatingDOAButton.Value ~= 1 )

        [data, valid, ~] = sdrReceiver();

        if valid

                data = zynqRadioDOAEstimationAmplitudeNormalization(data);

                % Rotating IQ data
                data = zynqRadioDOAEstimationRotateIQ4Channels(data,angle1,angle2,angle3);

                % Calculating the direction of arrival
                [y,doas] = estimator(data);
                pause(0.01);

                % Display angle in UI
                if isscalar(doas)
                    app.displayAngle(doas);
                end

                % Visualize in time domain
                if(app.DisplayTimeScopeCheckBox.Value)
                    timeScope([real(data), imag(data)]);
                    dataMaxLimit = max(max(abs([real(data); imag(data)])));
                    timeScope.YLimits = [ -dataMaxLimit*1.5, dataMaxLimit*1.5 ];
                end

                % Visualize frequency spectrum
                if(app.DisplaySpectrumScopeCheckBox.Value)
                    spectrumScope(data);
                end

                timeCounter = timeCounter + RadioFrameTime;
        end
    end
catch ME
    rethrow(ME);
end

Visualize Signal

When you select Display Time Scope, the app displays real and imaginary sinusoidal signals in the time domain. Use this information to check whether all signals are in synchronization in the initial condition (while using power splitter).

You can see the estimated direction of arrival angle in the semicircular gauge of the app.

Clear Up

Release the SDR receiver and visualization scopes, and delete the app handle.

release(sdrReceiver);
release(timeScope);
release(estimator);
release(spectrumScope);
delete(app);

Special Considerations

  1. To avoid multipath reflections, test this example in an open environment. Hold the transmitting and receiving antenna array at the same elevation. If you are working in a closed space, keep the transmitter in the near field of the antenna array.

  2. Make sure that elevation angle is near to zero. As the MUSIC algorithm estimates the broadside angle, the azimuth angle is only estimated when the elevation angle is near to zero.

  3. The script estimates the angle continuously from -70 to 70 degrees. When the transmitter is moved beyond these limits, discontinuities appear in the estimated angle.

  4. To increase the accuracy of the estimate, you must calibrate the antenna array setup.

See Also

Objects