phased.GLRTDetector

Generalized likelihood ratio detector

Since R2023b

Description

The generalized likelihood ratio test detector (GLRT) can detect signals with unknown parameters in the presence of noise. Unknown parameters include signal amplitude, phase, frequency, and arrival times. The detector replaces unknown parameters with their maximum likelihood estimates under the signal absent hypothesis H₀ or the alternative signal present hypothesis H₁ and then uses the LRT detector to output detection results.

1. Create the phased.GLRTDetector 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

detector = phased.GLRTDetector

detector = phased.GLRTDetector(Name = Value)

Description

detector = phased.GLRTDetector creates a generalized likelihood ratio test detector System object™ with default properties.

detector = phased.GLRTDetector(Name = Value) creates a likelihood ratio test detector System object with the specified property Name set to the specified Value. You can specify additional 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.

DataComplexity — Data complexity

'Complex' (default) | 'Real'

Data complexity, specified as 'Complex' or 'Real'. Complexity defines the format used to report output data. When DataComplexity is 'Complex', data is complex-valued. When DataComplexity is 'Real', data is real-valued.

Example: 'Real'

Data Types: char | string

ProbabilityFalseAlarm — Probability of false alarm

1e-1 (default) | nonnegative scalar

Probability of false alarm, specified as a positive scalar between 0 and 1, inclusive.

Example: 1e-6

Data Types: single | double

OutputFormat — Format of output data

'Detection result' (default) | 'Detection index'

Format of output data, specified as 'Detection result' or 'Detection index'. Output data is returned in the Y argument.

Example: 'Detection index'

Data Types: char | string

ThresholdOutputPort — Output detection statistics and threshold

false (default) | true

Output detection statistics and detection threshold, specified as false or true. Set this property to true to output the detection statistics in the stat argument and the detection threshold in the th argument. Set this property to false to suppress the output of detection statistics and threshold.

Data Types: logical

NoiseInputPort — Enable input of known noise power

false (default) | true

Enable the input of noise power, specified as false or true. You can input known or estimated noise power from a reference source. When true, use the input argument ncov to specify noise power. When false, the noise input port is disabled and the noise power is estimated by the GLRT detector.

Data Types: logical

SignalParameterOutputPort — Output maximum likelihood estimate of unknown signal parameters

false (default) | true

Output the maximum likelihood estimates of the unknown signal parameters under the alternative hypothesis (H₁), specified as false or true. Parameters are output in the estparam argument.

Data Types: logical

NoisePowerOutputPort — Output estimated noise power

false (default) | true

Output the estimated noise power, under the alternative hypothesis. Noise power is output in the estnoise argument.

Dependencies

To enable this property, set the NoiseInputPort property to false.

Data Types: logical

Usage

Syntax

Y = detector(X,hyp,obs)

[y,stat,th] = detector(___)

[Y,stst,th,estparam] = detector(___)

[Y,stst,th,estparam,estnoise] = detector(___)

Y = detector(X,hyp,obs,ncov)

Description

Y = detector(X,hyp,obs) performs generalized likelihood ratio test detection (GLRT) on the input data X. hyp is the augmented linear equality constraint matrix. obs is the observation matrix for a linear deterministic signal model. Detection results are returned in Y.

[y,stat,th] = detector(___) also returns the detection statistics, stat, and detection thresholds, th, of X. To enable this syntax, set the ThresholdOutputPort property to true. The format of stat depends on the OutputFormat property.

[Y,stst,th,estparam] = detector(___) also returns the maximum likelihood estimates of the unknown signal parameters estparam under the alternative hypothesis. To enable this syntax, set the SignalParameterOutputPort property to true.

[Y,stst,th,estparam,estnoise] = detector(___) also returns the estimated noise power, estnoise. To enable this syntax, set the NoisePowerOutputPort to true.

Y = detector(X,hyp,obs,ncov) also specifies the known noise power, ncov. To enable this syntax, set the NoiseInputPort property to true.

You can combine optional input and output arguments when their enabling properties are set. Optional inputs and outputs must be listed in the same order as the order of the enabling properties. For example, [Y,stat,th,estparam] = detector(X,hyp,obs,ncov).

Input Arguments

X — Input data

N-by-1 real-valued vector | N-by-1 complex-valued vector | N-by-M real-valued matrix | N-by-M complex-valued matrix

Input data, specified as a real-valued or complex-valued N-by-1 vector, or an N-by-M real-valued or complex-valued matrix. N is the signal length and M is the number of data channels. Detection is performed along the columns of X. The size of each row M cannot change during simulation. When M = 1, X represents a single channel of data. When M > 1, X can represent N samples from M data channels. Data streams can later be combined, for example, by beamforming.

The input data has a general interpretation. For example, the data can be interpreted as:

• Time series – N samples of a time series

• Sensor - N represents a snapshot of data samples from a set of sensors

Data Types: single | double
Complex Number Support: Yes

hyp — Augmented linear equality constraint matrix

real-valued R-by-(P + 1) matrix | complex-valued R-by-(P + 1) matrix

Augmented linear equality constraint matrix, specified as a real-valued or complex-valued R-by-(P + 1) matrix. The matrix takes the form [A,b] representing the equation:

$$A\Theta =b$$

where the unknown parameters are represented by ϴ. A has the rank R ≥ 1. The augmented linear equality constraint matrix expresses the null hypothesis H₀.

For this signal model, the GLRT detector determines whether to reject the null hypothesis which is expressed in the form A*θ = b where A is an R-by-P matrix with R ≤ P and rank R ≥ 1. b is an R-by-1 vector. A and b are carried in the augmented linear equality constraint matrix hyp = [A,b]. Because there are D signal models, the GLRT detector outputs D detection results for each column of X

Data Types: single | double
Complex Number Support: Yes

obs — Observation matrix

N-by-P-by-D array

Observation matrix for the linear deterministic signal model, specified as an N-by-P-by-D array N > p and rank P, D is the number of signal models, and the white Gaussian noise is a N-by-1 vector determined by the covariance ncov argument. The observation matrix is defined by X = obs*param + noise,

Example: 20.0

Data Types: single | double
Complex Number Support: Yes

ncov — Known noise power

positive scalar

Known noise power, specified as a scalar or a row vector of length D. When ncov is a scalar, it represents equal known noise power for D models. When ncov is a row vector of length D, it represents the known noise power for the D models, respectively.

Example: 20.0

Dependencies

To enable this argument, set the NoiseInputPort property to true.

Data Types: single | double

Output Arguments

Y — Detection results

D-by-M logical-valued vector | 1-by-L integer-valued vector | 2-by-L integer-valued matrix

Detection results of D models for M independent data samples returned as a D-by-M logical-valued vector, 1-by-L integer-valued vector, or 2-by-L integer-valued matrix. The format of Y depends on how the OutputFormat property is specified. By default, the OutputFormat property is set to 'Detection result'.

When OutputFormat is 'Detection result', Y is a D-by-M matrix containing logical detection results, where D is the number of signal models and M is the number of columns of X. For each row, Y is true in a column if there is a detection in the corresponding column of arg. Otherwise, Y is false.

When OutputFormat is 'Detection index', Y is a 1-by-L vector or a 2-by-L matrix containing detection indices, where L is the number of detection found in the M data samples and D models. When X is a column vector, Y is a 1-by-L vector and contains the index of the detections found in the D models. When X is a matrix, Y is a 2-by-L matrix, and each column of Y has the form [detrow;detcol], where detrow is the index of the model and detcol is the column index of X.

stat — detection statistic

N-by-M (default) | 1-by-L

Detection statistics, returned as a N-by-M matrix or 1-by-L vector. The format of stat depends on the OutputFormat property.

• When OutputFormat is 'Detection result', stat has the same size as Y.

• When OutputFormat is 'Detection index', stat is a 1-by-L vector containing detection statistics for each corresponding detection in Y.

th — Detection threshold

scalar

Detection threshold, returned as a scalar.

Dependencies

To enable this argument, set the ThresholdOutputPort property to true.

estparam — Maximum likelihood estimates of signal parameters

P-by-M-by-D array

Maximum likelihood estimates (MLE) of unknown signal parameters, returned as a P-by-M-by-D array.

Dependencies

To enable this argument, set the SignalParameterOutputPort property to true.

estnoise — Estimated noise power

positive scalar

Estimated noise power, returned as a positive scalar. When the OutputFormat property is 'Detection result', estnoise has the same size as Y. When OutputFormat property is 'Detection index', estnoise returns a noise power estimate of size 1-by-L for each corresponding detection in Y.

Dependencies

To enable this argument, set the NoisePowerOutputPort property to true.

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

Generalized Likelihood Ratio Test Detection on Real Data

Perform GLRT detection on a real Gaussian noise matrix with a desired probability of false alarm of 0.1. Assume that the signal dimension is 4, the noise power is unknown, and there are 2 unknown signal parameters with true values 0. Use a 3-model observation matrix. Perform the detection on all samples of the input and evaluate the probability of false alarm.

rng(10031);
glrt = phased.GLRTDetector( ...
    DataComplexity = 'Real', ...
    ProbabilityFalseAlarm = 0.1);
N = 4;
M = 1000;
p = 2;
D = 3;
x = randn(N,M);
hyp = [eye(p) zeros(p,1)];
obs = randn(N,p,D);
dresult = glrt(x,hyp,obs);
Pfa = sum(dresult,2)/M

Pfa = 3×1

    0.0970
    0.1030
    0.1110

Generalized Likelihood Ratio Threshold Detection on Complex Data

Perform Generalized Likelihood Ratio Threshold detection on complex Gaussian noise matrix with a desired probability of false alarm of 0.1. Assume that the signal dimension is 5. Also assume that the noise power is unknown. Assume there is only one unknown signal parameter with a true value zero. Use a single-model observation matrix. Perform the detection on all samples of the input and evaluate the probability of false alarm.

rng default
glrt = phased.GLRTDetector(DataComplexity = 'Complex', ...
    ProbabilityFalseAlarm = 0.1);
N = 5;
M = 1000;
x = 1/sqrt(2)*(randn(N,M)+1i*randn(N,M));
hyp = [1 0];
obs = ones(N,1);
dresult = glrt(x,hyp,obs);
Pfa = sum(dresult)/M

Pfa = 
0.1060

Perform Generalized Likelihood Ratio Test Detection

Perform GLRT detection on a complex Gaussian noise matrix with a desired probability of false alarm of 0.1. Assume that the signal dimension is 2 and that the noise power is unknown. There is only 1 unknown signal parameter with true value 0. Use a single-model observation matrix. Perform the detection on all samples of the input and evaluate the probability of false alarm.

Create the 1000 samples of Gaussian random data.

rng default
N = 2;
M = 1000;
x = 1/sqrt(2)*(randn(N,M) + 1i*randn(N,M));

Create the GLRT detector System object™.

glrt = phased.GLRTDetector(DataComplexity = 'complex', ...
    ProbabilityFalseAlarm = 0.1);

Specify the observation and hypotheses matrices.

hyp = [1 0];
obs = ones(N,1);

Solve and display the first 15 of the detection results.

dresult = glrt(x,hyp,obs);
disp(dresult(1:15))

   0   0   0   0   1   0   0   0   0   1   0   0   0   0   0

Estimate the probability of false alarm.

Pfa = sum(dresult)/M

Pfa = 
0.0930

The probability of false alarm is close to the desired probability of 0.1.

Perform Generalized Likelihood Ratio Test Detection For Two Parameters

Perform GLRT detection on a given complex Gaussian noise matrix with a desired probability of false alarm of 0.1. Assume that the signal dimension is 4 and the noise power is unknown. There are 2 unknown signal parameters with true values 0. Use a 3-model observation matrix. Perform the detection on all samples of the input and evaluate the probability of false alarm.

Create the random data.

rng default
N = 4;
M = 1000;
P = 2;
D = 3;
x = 1/sqrt(2)*(randn(N,M)+1i*randn(N,M));

Create the phased.GLRTDetector System object™.

glrt = phased.GLRTDetector(DataComplexity = 'complex', ...
    ProbabilityFalseAlarm = 0.1);
hyp = [eye(P) zeros(P,1)]

hyp = 2×3

     1     0     0
     0     1     0

obs = randn(N,P,D) + 1i*randn(N,P,D);

Run the detector and compute the probability of false alarm.

dresult = glrt(x,hyp,obs);
pfa = sum(dresult,2)/M

pfa = 3×1

    0.0890
    0.0980
    0.0970

Perform Angle-Domain GLRT Detection of Two Targets

Perform angle-domain GLRT detection for two targets located at 11 and -47 degrees with a desired probability of false alarm of 0.01. Use a ULA of dimension 128 to receive a single-snapshot signal with complex white Gaussian noise. The noise power and target complex amplitudes are unknown. Use a 256-model observation matrix, where each model checks the existence of a target following that model. Perform the detection on the input and plot the detection result.

rng default
glrt = phased.GLRTDetector(DataComplexity = 'complex', ...
    ProbabilityFalseAlarm = 0.01);

Create two targets located at 11 and -47 degrees.

tarAng = [11,-47];
tarNum = length(tarAng);

Start with a ULA array with 128 elements having 1/2-wavelength element spacing.

N = 128;
fc = 3e8;
elementPos = (0:N-1)*physconst('LightSpeed')/fc/2;

Create single-snapshot received signal from two targets at the ULA.

tarAmp = 1/sqrt(2)*(randn(tarNum,1) + 1i*randn(tarNum,1));
signal = steervec(elementPos,tarAng)*tarAmp;

Single-snapshot of complex white Gaussian noise at the ULA.

noise = 1/sqrt(2)*(randn(N,1) + 1i*randn(N,1));

Add the noise to the signal.

x = signal + noise;

Partition angle into 256 angle bins

D = 256;
angGrid = asind(2*(-D/2:D/2-1)/D);

Perform GLRT detection on 256 signal models. Each model assumes one target. % N-by-1-by-D

obs = permute(steervec(elementPos,angGrid),[1 3 2]);
hyp = [1,0];
dresult = glrt(x,hyp,obs);

Plot the GLRT detection results at the 256 angle bins.

plot(angGrid,dresult)
xlabel('Angle (degrees)')
ylabel('Detection Result')
grid on

Perform Angle-Doppler-Domain GLRT Detection

Perform angle-Doppler-domain GLRT detection for a target located at 28 degree at a Doppler of 10 Hz with a desired probability of false alarm of 1e-6. Use a ULA of dimension 128 to receive a 32-pulse signal under complex white Gaussian noise. The noise power and target complex amplitude are unknown. Use an observation matrix with 256-by-64 angle-Doppler models, where each model checks the existence of a target following that model. Perform the detection on the input and plot the spatial-temporal detection map.

rng default
glrt = phased.GLRTDetector('DataComplexity','complex',...
    'ProbabilityFalseAlarm',1e-6,...
    'ThresholdOutputPort',true);

Target is located at 28 degree at a Doppler of 10 Hz.

tarAng = -28;
tarDop = 10;

Create a ULA with 128 elements and with elements spaced at 1/20-wavelength.

elementNum = 128;
fc = 3e8;
elementPos = (0:elementNum-1)*physconst('LightSpeed')/fc/2;

ULA observes 32 pulses at a PRF of 50 Hz

fPRF = 50;
pulseNum = 32;

Create the spatial-temporal signal from the target at the ULA.

tarAmp = 1/sqrt(2)*(randn + 1i*randn);
dopSv = dopsteeringvec(tarDop,pulseNum,fPRF);
angSv = steervec(elementPos,tarAng);
signal = kron(dopSv,angSv)*tarAmp;

Add the spatial-temporal complex white Gaussian noise at the ULA.

noisePow = 0.01;
N = elementNum*pulseNum;
noise = sqrt(noisePow/2)*(randn(N,1) + 1i*randn(N,1));
x = signal + noise;

Partition angle into 256 angle bins and Doppler into 64 Doppler bins

angNum = 256;
dopNum = 64;
angGrid = asind(2*(-angNum/2:angNum/2-1)/angNum);
dopGrid = fPRF*(-dopNum/2:dopNum/2-1)/dopNum; D = angNum*dopNum;

Perform GLRT detection on 256-by-64 angle-Doppler models, each assuming 1 target.

obs = zeros(N,1,D);
for dopBin = 1:dopNum
    dopGridSv = dopsteeringvec(dopGrid(dopBin),pulseNum,fPRF);
    for angBin = 1:angNum
        angGridSv = steervec(elementPos,angGrid(angBin));
        angDopGridSv = kron(dopGridSv,angGridSv);
        obs(:,:,(dopBin-1)*angNum+angBin) = angDopGridSv;
    end
end
hyp = [1,0];
[dresult,stat,th] = glrt(x,hyp,obs);
stat = reshape(stat,angNum,dopNum);

Plot a spatial-temporal GLRT detection map.

[dopMesh,angMesh] = meshgrid(dopGrid,angGrid);
surf(angMesh,dopMesh,pow2db(stat))
hold on
mesh(angMesh,dopMesh,pow2db(th)*ones(angNum,dopNum))
xlabel('Angle (degrees)')
ylabel('Doppler (Hz)')
zlabel('Power (dB)')

More About

Data Precision

This System object supports single and double precision for input data, properties and arguments. If input data is single precision, all the non-logical outputs are also single precision. If input data is double precision, the non-logical outputs are double precision. The precision of the non-logical output is independent of the precision of the properties and other arguments.

Signal Model

The GLRT detector assumes D linear deterministic signals in each column of X. The GLRT detector assumes that the data follows the linear deterministic signal model X = H*param + noise where X is an N-by-1 data vector, H the observation matrix given as an N-by-P-by-D array.

• param is a P-by-1 vector that contains P unknown signal parameters.

• H is an N-by-P-by-D array that summarizes D observation matrices of size N-by-P with N > p\P and rank P.

• noise is an N-by-1 noise vector assumed to be white Gaussian.

For this signal model, the GLRT detector determines whether to reject the null hypothesis. The null hypothesis is expressed in the form A*param = b where A is an R-by-P matrix with R ≤ P and rank R ≥ 1, and b is an R-by-1 vector. A and b are form the augmented linear equality constraint matrix hyp = [A,b]. Because there are D signal models, the GLRT detector outputs D detection results for each column of X.

H is the observation matrix for a linear deterministic signal model X = H*param + noise. obs is an N-by-P-by-D array with N > p and with rank p, D is the number of signal models, and noise is a white Gaussian noise N-by-1 vector.

The GLRT algorithm replaces the unknown parameters by their maximum likelihood estimates under the H₀ or H₁ hypotheses. Subsequently, the algorithm applies the likelihood ratio test. For a data vector x with unknown model parameters ϴ, the GLRT decides H₁ if the likelihood ratio is greater than the threshold η

$$\frac{p(x,{\widehat{\Theta }}_1,H_1)}{p(x,{\widehat{\Theta }}_0,H_0)}>\eta$$

for the estimated ϴ.

Detection Statistics and Threshold

For a real data model with known Gaussian white noise, the test statistic is given by

$$T(x)=\frac{{\left(A{\widehat{\theta }}_1-b\right)}^T{\left(A{(H^{T}H)}^{-1}{A}^T\right)}^{-1}\left(A{\widehat{\theta }}_1-b\right)}{{\sigma }^2}$$

where $${\widehat{\theta }}_1$$ represents the maximum likelihood estimate of θ given by

$${\widehat{\theta }}_1={\left(H^{T}H\right)}^{-1}{H}^{T}x$$

For a complex data model with known Gaussian noise, the test statistic is given by

$$T(x)=2\frac{{\left(A{\widehat{\theta }}_1-b\right)}^{†}{\left(A{(H^{†}H)}^{-1}{A}^{†}\right)}^{-1}\left(A{\widehat{\theta }}_1-b\right)}{{\sigma }^2}$$

where H^† denotes the adjoint of H where

$${\widehat{\theta }}_1={\left(H^{†}H\right)}^{-1}{H}^{†}x$$

For the real data model with unknown white Gaussian noise. the test statistic is

$$T(x)=\frac{1}{r}\frac{{\left(A{\widehat{\theta }}_1-b\right)}^T{\left(A{(H^{T}H)}^{-1}{A}^T\right)}^{-1}\left(A{\widehat{\theta }}_1-b\right)}{{\widehat{\sigma }}^2}$$

where again the maximum likelihood estimate of θ is represented by $${\widehat{\theta }}_1$$

$${\widehat{\theta }}_1={\left(H^{T}H\right)}^{-1}{H}^{T}x$$

. but here σ² is estimated by

$$\begin{array}{l}{\widehat{\sigma }}^2=\frac{1}{N-p}{\Vert P_H^{\perp }x\Vert }^2\\ P_H^{\perp }=I-H{(}^{H^T}{H}^T\end{array}$$

For complex data with unknown Gaussian noise, the test statistic is

$$T(x)=\frac{1}{r}\frac{{\left(A{\widehat{\theta }}_1-b\right)}^{†}{\left(A{(H^{†}H)}^{-1}{A}^{†}\right)}^{-1}\left(A{\widehat{\theta }}_1-b\right)}{{\widehat{\sigma }}^2}$$

where

$${\widehat{\theta }}_1={\left(H^{†}H\right)}^{-1}{H}^{†}x$$

and the estimate of σ² is

$$\begin{array}{l}{\widehat{\sigma }}^2=\frac{1}{N-p}{\Vert P_H^{\perp }x\Vert }^2\\ P_H^{\perp }=I-H{(}^{H^H}{H}^H\end{array}$$

The detection threshold for real data with known white Gaussian noise and a given theoretical probability of false alarm P_fa is calculated from the chi-square inverse cumulative distribution function (CDF) with r degrees of freedom.

$$\gamma =\text{chi2inv}(1-P_{\text{fa}},r)$$

For complex data with known white Gaussian noise and a given theoretical probability of false alarm P_fa, the theoretical detection threshold is calculated from the chi-square inverse cumulative distribution function (CDF) with 2r degrees of freedom

$$\gamma =\text{chi2inv}(1-P_{\text{fa}},2r)$$

For real data with unknown white Gaussian noise and a given theoretical probability of false alarm P_fa, the theoretical detection threshold calculated from the f-inverse cumulative distribution function (CDF) with r degrees of freedom

$$\gamma =\text{finv}(1-P_{\text{fa}},r,N-p)$$

For complex data with unknown white Gaussian noise and a given theoretical probability of false alarm P_fa, the theoretical detection threshold is calculated from the f-inverse cumulative distribution function (CDF) with 2r degrees of freedom

$$\gamma =\text{finv}(1-P_{\text{fa}},2r,2(N-p))$$

References

[1] Steven M. Kay, Fundamentals of Statistical Signal Processing, Detection Theory, Prentice-Hall PTR, 1993.

[2] Mark A. Richards, Fundamentals of Radar Signal Processing, Third edition, McGraw-Hill Education, 2022.

Extended Capabilities

C/C++ Code Generation
Generate C and C++ code using MATLAB® Coder™.

Version History

Introduced in R2023b