phased.GCCEstimator

Wideband direction of arrival estimation

Description

The phased.GCCEstimator System object™ creates a direction of arrival estimator for wideband signals. This System object estimates the direction of arrival or time of arrival among sensor array elements using the generalized cross-correlation with phase transform algorithm (GCC-PHAT). The algorithm assumes that all signals propagate from a single source lying in the array far field so the direction of arrival is the same for all sensors. The System object first estimates the correlations between all specified sensor pairs using GCC-PHAT and then finds the largest peak in each correlation. The peak identifies the delay between the signals arriving at each sensor pair. Finally, a least-squares estimate is used to derive the direction of arrival from all estimated delays.

To compute the direction of arrival for pairs of element in the array:

1. Create the phased.GCCEstimator object and set its properties.

2. Call the object with arguments, as if it were a function.

To learn more about how System objects work, see What Are System Objects?

Creation

Syntax

sGCC = phased.GCCEstimator

sGCC = phased.GCCEstimator(Name=Value)

Description

sGCC = phased.GCCEstimator creates a GCC direction of arrival estimator System object, sGCC. This object estimates the direction of arrival or time of arrival between sensor array elements using the GCC-PHAT algorithm.

sGCC = phased.GCCEstimator(Name=Value) returns a GCC direction of arrival estimator object, sGCC, with the specified property Name set to the specified Value. You can specify several name-value pair arguments in any order as Name1=Value1,...,NameN=ValueN.

Properties

Unless otherwise indicated, properties are nontunable, which means you cannot change their values after calling the object. Objects lock when you call them, and the release function unlocks them.

If a property is tunable, you can change its value at any time.

For more information on changing property values, see System Design in MATLAB Using System Objects.

sGCC — GCC-PHAT estimator

phased.GCCEstimator System object

GCC-PHAT estimator, specified as a phased.GCCEstimator System object.

Example: phased.GCCEstimator

SensorArray — Sensor array

phased.ULA System object (default) | Phased Array System Toolbox™ sensor array

Sensor array, specified as a Phased Array System Toolbox System object. The array can also consist of subarrays. If you do not specify this property, the default sensor array is a phased.ULA System object with default array property values.

Example: phased.URA

PropagationSpeed — Signal propagation speed (meters/second)

physconst("LightSpeed") (default) | real-valued positive scalar

Signal propagation speed, specified as a real-valued positive scalar. Units are in meters per second. The default propagation speed is the value returned by physconst("LightSpeed").

Example: 3e8

Data Types: single | double

SampleRate — Signal sample rate

1e6 (default) | positive real-valued scalar

Signal sample rate, specified as a positive real-valued scalar. Units are in hertz.

Example: 1e6

Data Types: single | double

SensorPairSource — Source of sensor pairs

"Auto" (default) | "Property"

Source of sensor pairs, specified as either "Auto" or "Property".

• "Auto" — choose this property value to compute correlations between the first sensor and all other sensors. The first sensor of the array is the reference channel.

• "Property" — choose this property value when you want to explicitly specify the sensor pairs to be used for computing correlations. Set the sensor pair indices using the SensorPair property. You can view the array indices using the viewArray method of any array System object.

Example: "Auto"

Data Types: char

SensorPair — Sensor pairs

[2;1] (default) | 2-by-N positive integer valued matrix

Sensor pairs used to compute correlations, specified as a 2-by-N positive integer-valued matrix. Each column of the matrix specifies a pair of sensors between which the correlation is computed. The second row specifies the reference sensors. Each entry in the matrix must be less than the number of array sensors or subarrays. To use the SensorPair property, you must also set the SensorPairSource value to "Property".

Example: [1,3,5;2,4,6]

Data Types: double

DelayOutputPort — Option to enable delay output

false (default) | true

Option to enable output of time delay values, specified as a Boolean. Set this property to true to output the delay values as an output argument of the step method. The delays correspond to the arrival angles of a signal between sensor pairs. Set this property to false to disable the output of delays.

Example: false

Data Types: logical

CorrelationOutputPort — Option to enable correlation output

false (default) | true

Option to enable output of correlation values, specified as a Boolean. Set this property to true to output the correlations and lags between sensor pairs as output arguments of the step method. Set this property to false to disable output of correlations.

Example: false

Data Types: logical

Usage

Syntax

ang = sGCC(X)

[ang,tau] = sGCC(X)

[ang,R,lag] = sGCC(X)

[ang,tau,R,lag] = sGCC(X)

Description

ang = sGCC(X) returns the direction of arrival, ang, of an input signal X. The argument X is a matrix specifying the received signals at the elements of the array specified in the SensorArray property. Signals propagate from a single source. Each column in X corresponds to the elements in the array (if an array is used) or the number of subarrays (if a subarray is used). Each row of X represents a single time snapshot.

[ang,tau] = sGCC(X) returns the time delays, tau, estimated from the correlations between pairs of sensors. To use this syntax, set the DelayOutputPort property to true. The argument tau is a P-element row vector, where P is the number of sensor pairs, and where P = N(N-1).

[ang,R,lag] = sGCC(X) returns the estimated correlations, R, between pairs of sensors, when you set the CorrelationOutputPort property to true. R is a matrix with P columns where P is the number of sensor pairs. Each column in R contains the correlation for the corresponding pair of sensors. lag is a column vector containing the time lags corresponding to the rows of the correlation matrix. The time lags are the same for all sensor pairs.

You can combine optional input arguments when their enabling properties are set. Optional inputs must be listed in the same order as their enabling properties. For example,[ang,tau,R,lag] = sGCC(X) is valid when you set both DelayOutputPort and CorrelationOutputPort to true.

Note

The object performs an initialization the first time the object is executed. This initialization locks nontunable properties and input specifications, such as dimensions, complexity, and data type of the input data. If you change a nontunable property or an input specification, the System object issues an error. To change nontunable properties or inputs, you must first call the release method to unlock the object.

Input Arguments

X — Received signal

M-by-N complex-valued matrix

Received signal, specified as an M-by-N complex-valued matrix. The quantity M is the number of sample values (snapshots) of the signal and N is the number of sensor elements in the array. For subarrays, N is the number of subarrays.

The size of the first dimension of the input matrix can vary to simulate a changing signal length. A size change can occur, for example, in the case of a pulse waveform with variable pulse repetition frequency.

Example: [[0;1;2;3;4;3;2;1;0],[1;2;3;4;3;2;1;0;0]]

Data Types: single | double
Complex Number Support: Yes

Output Arguments

ang — Direction of arrival

2-by-1 real-valued column vector | scalar

Direction of arrival of a signal, returned as a 2-by-1 real-valued column vector in the form [azimuth;elevation]. If the array is a uniform linear array, ang is a scalar representing the broadside angle. Angle units are in degrees, defined with respect to the local coordinate system of the array.

tau — Time delays of arrival

1-by-P real-valued row vector

Time delays of arrival, returned as 1-by-P real-valued row vector. P is the number of sensor pairs selected from the array.

• When SensorPairSource is set to "Auto", P = N - 1. N is the number of elements in the array.

• When SensorPairSource is set to "Property", P is the number of sensor pairs specified by the SensorPair property.

Time units are seconds. This output is enabled when you set the DelayOutputPort property to true.

R — Estimated cross-correlation

(2M+1)-by-P complex-valued matrix

Estimated cross-correlation between pairs of sensors, returned as a (2M+1)-by-P complex-valued matrix, where P is the number of sensor pairs selected from the array.

• When SensorPairSource is set to "Auto", P = N - 1. N is the number of elements in the array. The columns in R contain the correlations between the first sensor and all other sensors.

• When SensorPairSource is set to "Property", P is the number of sensor pairs specified by the SensorPair property. Each column in R contains the correlation for the corresponding pair of sensors.

M is the number of time samples in the input signal. This output is enabled when you set the CorrelationOutputPort property to true.

lag — Time lags

M-by-1 real-valued column vector

Time lags, returned as an (2M+1)-by-1 real-valued column vector. M is the number of time samples in the input signal. Each time lag applies to the corresponding row in the cross-correlation matrix.

Object Functions

To use an object function, specify the System object as the first input argument. For example, to release system resources of a System object named obj, use this syntax:

release(obj)

Common to All System Objects

step	Run System object algorithm
release	Release resources and allow changes to System object property values and input characteristics
reset	Reset internal states of System object

Examples

GCC Estimate of Direction of Arrival at Microphone Array

Estimate the direction of arrival of a signal using the GCC-PHAT algorithm. The receiving array is a 5-by-5-element URA microphone array with elements spaced 0.25 meters apart. The arriving signal is a sequence of wideband chirps. The signal arrives from 17° azimuth and 0° elevation. Assume the speed of sound in air is 340 m/s.

Load the chirp signal.

load chirp;
c = 340.0;

Create the 5-by-5 microphone URA.

d = 0.25;
N = 5;
mic = phased.OmnidirectionalMicrophoneElement;
array = phased.URA([N,N],[d,d],Element=mic);

Simulate the incoming signal using the WidebandCollector System object™.

arrivalAng = [17;0];
collector = phased.WidebandCollector(Sensor=array,PropagationSpeed=c, ...
    SampleRate=Fs,ModulatedInput=false);
signal = collector(y,arrivalAng);

Estimate the direction of arrival.

estimator = phased.GCCEstimator(SensorArray=array, ...
    PropagationSpeed=c,SampleRate=Fs);
ang = estimator(signal)

ang = 2×1

   16.4538
   -0.7145

Plot Correlations from GCC Estimator

Estimate the direction of arrival of a signal using GCC-PHAT. The receiving array is a 5-by-5-element URA microphone array with elements spaced 25 cm apart. Choose 10 element pairs to compute the arrival angle. Assume the speed of sound in air is 340 m/sec. The arriving signal is a sequence of wideband sounds. Assume the signal arrives from ${54}^{\circ }$ azimuth and $5^{\circ }$ elevation. Estimate the arrival angle, and then plot the correlation function versus lag for one pair of elements.

Load the signal and extract a small portion (5000 samples) for computation.

load gong
y1 = y(1:5000);

Set up the receiving array.

N = 5;
d = 0.25;
mic = phased.OmnidirectionalMicrophoneElement;
array = phased.URA([N,N],[d,d],Element=mic);

Simulate the arriving plane wave using the WidebandCollector System object™.

c = 340.0;
arrivalAng = [54;5];
coll = phased.WidebandCollector(Sensor=array, ...
    PropagationSpeed=c, ...
    SampleRate=Fs, ...
    ModulatedInput=false);
signal = coll(y1,arrivalAng);

Estimate direction of arrival. Choose 10 sensors to correlate with the first element of the URA.

sensorpairs = [[2,4,6,8,10,12,14,16,18,20];ones(1,10)];
estimator = phased.GCCEstimator(SensorArray=array, ...
    PropagationSpeed=c,SampleRate=Fs, ...
    SensorPairSource="Property", ...
    SensorPair=sensorpairs, ...
    DelayOutputPort=true,CorrelationOutputPort=true);
[estimatedAng,taus,R,lags] = estimator(signal);

The estimated angle is:

disp(estimatedAng)

   53.7028
    4.9552

Plot the correlation between sensor 8 and sensor 1. This pair corresponds to the fourth column of the correlation matrix. The estimated value of tau (in msec ) for this pair is:

disp(1000*taus(4))

   -0.4757

plot(1000*lags,real(R(:,4)))
xlabel("Time lags (msec)")
ylabel("Correlation")

Algorithms

GCC-PHAT Cross-Correlation Algorithm

You can use generalized cross-correlation to estimate the time difference of arrival (TDOA) of a signal at two different sensors.

To estimate TDOA, find the locations of the peak of the cross-correlation function between the received signal at one sensor and the received signal at a second sensor. The peak of the cross-correlation function corresponds to the time delay. Both the phased.GCCEstimator and phased.TDOAEstimator System objects compute correlation peaks using the Generalized Cross-Correlation with Phase-Transformation (GCC-PHAT) algorithm. Peak locations are restricted to multiples of the sample interval which is the inverse of the signal sample rate. You can specify the signal sample rate by the SampleRate property which is common to both System objects.

To compute the time delay between two signals, x₁(t) and x₂(t), start with their time-domain representations in the presence of noise:

$$\begin{array}{l}{x}_1(t)=s(t)+n_1(t)\\ x_2(t)=\alpha s(t-D)+n_2(t)\end{array}$$

where n₁(t) and n₂(t) are the noises added to each signal. D is the time-delay.

In the frequency domain, the signals are represented as

$$\begin{array}{l}{X}_1(\omega )=S(\omega )+N_1(\omega )\\ X_2(\omega )=\alpha S(\omega )e^{-j\omega D}+N_2(\omega )\end{array}$$

Define the cross-correlation function between two signals in the time domain and in the frequency domain.

$$\begin{array}{c}{R}_{x_1,x_2}(\tau )=\text{E}[x_1(t)x_2(t-\tau )]\\ R_{x_1,x_2}(\omega )=X_1(\omega )X_2^{*}(\omega )\\ =\alpha {\left|S(\omega )\right|}^2{e}^{j\omega \tau }\end{array}$$

To identify the time delay, locate the peak of the cross-correlation function τ_max. When the signal-to-noise ratio (SNR) is large, the correlation peak τ_max corresponds to the actual time delay D. When the correlation function is more sharply peaked, performance improves. You can sharpen the cross correlation peak by using a weighting function that whitens the input signals. This technique is called generalized cross-correlation (GCC). One particular weighting function normalizes the signal spectral density with the spectrum magnitude, leading to the generalized cross-correlation phase transform method (GCC-PHAT). The Generalized Cross-Correlation function with Phase Transform (GCC-PHAT) at lag τ is defined as

$$\begin{array}{c}{\tilde{R}}_{x_1{x}_2}(\tau )=\frac{1}{2\pi }{\displaystyle {\int }_{-\pi }^{\pi }\frac{R_{x_1{x}_2}(\omega )}{\left|R_{x_1{x}_2}(\omega )\right|}}{e}^{i\omega \tau }d\omega \\ =\frac{1}{2\pi }{\displaystyle {\int }_{-\pi }^{\pi }{e}^{i\omega (D+\tau )}}d\omega \\ =\delta (D-\tau )\end{array}$$

where δ is the Dirac-delta function which is non-zero at D=τ, which is the time-difference of arrival (TDOA). Typically, the lag is calculated in samples, so to convert that to seconds, divide by the sampling rate. You can estimate the time delay by finding the time lag that maximizes the cross-correlation between the two signals.

From the TDOA, you can estimate the broadside arrival angle of the plane wave with respect to the line connecting the two sensors. For two sensors separated by distance L, the broadside arrival angle, Broadside Angles, is related to the time lag by

$$\sin \beta =\frac{c\tau }{L}$$

where c is the propagation speed in the medium.

If you use just one sensor pair, you can only estimate the broadside angle of arrival. However, if you use multiple pairs of non-collinear sensors, for example, in a URA, you can estimate the arrival azimuth and elevation angles of the plane wave using least-square estimation. For N sensors, you can write the delay time τ_kj of a signal arriving at the k^th sensor with respect to the j^th sensor by

$$\begin{array}{l}c{\tau }_{kj}=-\left({\overrightarrow{x}}_k-{\overrightarrow{x}}_j\right)\cdot \overrightarrow{u}\\ \overrightarrow{u}=\cos \alpha \sin \theta \widehat{i}+\sin \alpha \sin \theta \widehat{j}+\cos \theta \widehat{k}\end{array}$$

where u is the unit propagation vector of the plane wave. The angles α and θ are the azimuth and elevation angles of the propagation vector. All angles and vectors are defined with respect to the local axes. You can solve the first equation using least-squares to yield the three components of the unit propagation vector. Then, you can solve the second equation for the azimuth and elevation angles.

Data Precision

This System object supports single and double precision for input data, properties, and arguments. If the input data X is single precision, the output data is single precision. If the input data X is double precision, the output data is double precision. The precision of the output is independent of the precision of the properties and other arguments.

References

[1] Knapp, C. H. and G.C. Carter, “The Generalized Correlation Method for Estimation of Time Delay.” IEEE Transactions on Acoustics, Speech and Signal Processing. Vol. ASSP-24, No. 4, Aug 1976.

[2] G. C. Carter, “Coherence and Time Delay Estimation.” Proceedings of the IEEE. Vol. 75, No. 2, Feb 1987.

Extended Capabilities

C/C++ Code Generation

Generate C and C++ code using MATLAB® Coder™.

Usage notes and limitations:

See System Objects in MATLAB Code Generation (MATLAB Coder).

This System object supports single and double precision for input data, properties, and arguments. If the input data X is single precision, the output data is single precision. If the input data X is double precision, the output data is double precision. The precision of the output is independent of the precision of the properties and other arguments.

Version History

Introduced in R2015b