pilotcalib

Array calibration using pilot sources

Syntax

estpos = pilotcalib(nompos,x,pilotang)

[estpos,esttaper] = pilotcalib(nompos,x,pilotang)

[estpos,esttaper] = pilotcalib(nompos,x,pilotang,nomtaper)

[estpos,esttaper] = pilotcalib(nompos,x,pilotang,nomtaper,uncerts)

Description

estpos = pilotcalib(nompos,x,pilotang) returns the estimated element positions, estpos, of a sensor array. The argument nompos represents the relative nominal positions of the sensor array before calibration. The nominal position is relative to the first element of the array. The argument x represents the signals received by the array coming from the pilot sources. The argument pilotang contains the known directions of each of the pilot sources. Three or more pilot sources are required in this case.

[estpos,esttaper] = pilotcalib(nompos,x,pilotang) also returns the estimated array taper, esttaper. Each element of esttaper contains the estimated taper value of the corresponding array element. In this case, the prior nominal taper is one for each element. Four or more pilot sources are required in this case.

[estpos,esttaper] = pilotcalib(nompos,x,pilotang,nomtaper) specifies nomtaper as the nominal taper of the array. Four or more pilot sources are required in this case.

[estpos,esttaper] = pilotcalib(nompos,x,pilotang,nomtaper,uncerts) specifies uncerts as the configuration settings to use for calibrating array uncertainties. Configuration settings determine which parameters to estimate.

Examples

Estimate ULA Element Positions Using Pilot Calibration

Construct a 7-element ULA array of isotropic antenna elements spaced one-half wavelength apart. Assume the array is geometrically perturbed in three dimensions. Perform pilot calibration on the array using 4 pilot sources at azimuth and elevation angles of (-60,0), (10,-40), (40,0), and (120,45) degrees. For the calibration process, pilot signals have an SNR of 30 dB. Each pilot signal contains 10,000 samples. Assume the signals have a frequency of 600 MHz.

Set up the ULA with nominal parameters

fc = 600e6;
c = physconst('LightSpeed');
lam = c/fc;
d = 0.5*lam;
sIso = phased.IsotropicAntennaElement('FrequencyRange',[100,900]*1e6);
Nelem = 7;
NominalTaper = ones(1,Nelem);
sULA = phased.ULA('Element',sIso,'NumElements',Nelem,'ElementSpacing',d,...
    'Taper',NominalTaper);

Create the pilot signals

Randomly perturb the element positions with a gaussian distribution having 0.1 wavelength standard deviation. Do not perturb the position of the first element or the tapers.

posstd = 0.1;
rng default
NominalElementPositions = getElementPosition(sULA)/lam;
ReferenceElement = NominalElementPositions(:,1);
PositionPert = [zeros(3,1),posstd*randn(3,Nelem-1)];
ActualElementPositions = NominalElementPositions + PositionPert;
ActualTaper = NominalTaper;

Generate the signals using the actual positions and tapers.

Nsamp  = 10000;
ncov = 0.001;
PilotAng = [-60,10,40,120; 0,-40,0,45];
Npilot = size(PilotAng,2);
for n = 1:Npilot
    X(:,:,n) = sensorsig(ActualElementPositions,...
        Nsamp,PilotAng(:,n),ncov,'Taper',ActualTaper.');
end

Perform the pilot calibration

estpos = pilotcalib(NominalElementPositions - ReferenceElement*ones(1,Nelem),...
    X,PilotAng);

Add back the position of the reference sensor

estpos = estpos + NominalElementPositions(:,1)*ones(1,Nelem);

Examine the root mean squared (RMS) error of the calibrated parameters

Compute the RMS value of the initial position error.

numpos = 3*Nelem;
initposRMSE = sqrt(sum(PositionPert(:).^2)/numpos);

Compute the RMS value of the calibrated position error.

solvposErr = ActualElementPositions - estpos;
solvposRMSE = sqrt(sum(solvposErr(:).^2)/(numpos));

Compare the calibrated RMS position error to the initial position RMS error. The calibration reduces the RMS position error.

disp(solvposRMSE/initposRMSE)

   2.3493e-04

Estimate ULA Element Position and Taper Errors Using Pilot Calibration

Construct a 7-element ULA array of isotropic antenna elements spaced one-half wavelength apart. Assume the array is geometrically perturbed in three dimensions. Perform pilot calibration on the array using 4 pilot sources at azimuth and elevation angles of (-60,0), (10,80), (40,-40), and (-80,0) degrees. For the calibration process, pilot signals have an SNR of 30 dB. Each pilot signal contains 10,000 samples. Assume the signals have a frequency of 600 MHz.

Set up the ULA with nominal parameters

fc = 600e6;
c = physconst('LightSpeed');
lam = c/fc;
d = 0.5*lam;
sIso = phased.IsotropicAntennaElement('FrequencyRange',[100,900]*1e6);
Nelem = 7;
NominalTaper = ones(1,Nelem);
sULA = phased.ULA('Element',sIso,'NumElements',Nelem,'ElementSpacing',d,...
    'Taper',NominalTaper);

Create the pilot signals

Randomly perturb the element positions using a Gaussian distribution that has a standard deviation of 0.1 wavelength. Do not perturb the position of the first element.

posstd = 0.1;
rng default
NominalElementPositions = getElementPosition(sULA)/lam;
ReferenceElement = NominalElementPositions(:,1);
PositionPert = [zeros(3,1),posstd*randn(3,Nelem-1)];
ActualElementPositions = NominalElementPositions + PositionPert;

Perturb the taper in magnitude and phase. Do not perturb the first taper.

tapermagstd = 0.15;
taperphasestd = 0.15;
tapermagpert = tapermagstd*[0; randn(Nelem-1,1)];
ActualTaper = NominalTaper' + tapermagpert;
taperphasepert = taperphasestd*[0;randn(Nelem-1,1)];
ActualTaper = ActualTaper.*exp(1i*taperphasepert);

Generate the signals using the perturbed positions, tapers and four pilot sources.

Nsamp  = 10000;
ncov = 0.001;
PilotAng = [-60,10,40,-80; 10,80,-40,0];
Npilot = size(PilotAng,2);
for n = 1:Npilot
    X(:,:,n) = sensorsig(ActualElementPositions,Nsamp,...,
        PilotAng(:,n),ncov,'Taper',ActualTaper);
end

Perform the pilot calibration

[estpos,esttaper] = pilotcalib(...
    NominalElementPositions - ReferenceElement*ones(1,Nelem),...
    X,PilotAng);

Add back the position of the reference sensor

estpos = estpos + NominalElementPositions(:,1)*ones(1,Nelem);

Examine the root mean square (RMS) error of the calibrated parameters

Compute the RMS values of the initial taper perturbations.

tapermagpertRMSE = sqrt(tapermagpert'*tapermagpert/Nelem);
taperphasepertRMSE = sqrt(taperphasepert'*taperphasepert/Nelem);

Compute the RMS value of the calibrated taper magnitude error.

diff = abs(ActualTaper) - abs(esttaper);
diff2 = diff'*diff;
tapermagsolvRMSE = sqrt(diff2/Nelem);

Compare the calibrated RMS magnitude error to the initial RMS magnitude error. The calibration reduces the RMS magnitude error.

disp(tapermagsolvRMSE/tapermagpertRMSE)

   6.7715e-04

Compute the RMS value of the calibrated taper phase error.

diff = unwrap(angle(ActualTaper) - angle(esttaper));
diff2 = diff'*diff;
tapersolvphaseRMSE = sqrt(diff2/Nelem);

Compare the calibrated RMS phase error to the initial RMS phase error. The calibration reduces the RMS phase error.

disp(tapersolvphaseRMSE/taperphasepertRMSE)

    0.0021

% Compute the RMS value of the initial position error.
numpos = 3*Nelem;
initposRMSE = sqrt(sum(PositionPert(:).^2)/numpos);

Compute the RMS value of the calibrated position error.

solvposErr = ActualElementPositions - estpos;
solvposRMSE = sqrt(sum(solvposErr(:).^2)/(numpos));

Compare the calibrated RMS position error to the initial position RMS error. The calibration reduces the RMS position error.

disp(solvposRMSE/initposRMSE)

   3.6308e-04

Estimate URA Element Position Errors Using Pilot Calibration

Construct a 9-element URA of isotropic antenna elements spaced one-half wavelength apart. Assume the array has been geometrically perturbed in all directions except for the first element. Perform pilot calibration on the array using 5 pilot sources at azimuth and elevation angles of (-60,0), (10,-40), (40,0), (120,45), and (170,50) degrees. For the calibration process, pilot signals have an SNR of 30 dB. Each pilot signal contains 10,000 samples. Assume the signals have a frequency of 600 MHz.

Create the array

For convenience, use a phased.URA System object� to set the nominal position and taper values.

fc = 300e6;
c = physconst('LightSpeed');
lam = c/fc;
d = 0.5*lam;
sIso = phased.IsotropicAntennaElement('FrequencyRange',[100,900]*1e6);
sURA = phased.URA('Element',sIso,'Size',[3,3],...
    'ElementSpacing',d,'Taper',ones(3,3));
Nelem = getNumElements(sURA);
taper = getTaper(sURA);

Create the pilot signals

Randomly perturb the element positions using a Gaussian distribution that has a standard deviation of 0.1 wavelength. Do not perturb the position of the first element.

posstd = 0.1;
rng default
NominalElementPositions = getElementPosition(sURA)/lam;
ReferenceElement = NominalElementPositions(:,1);
PositionPert = [zeros(3,1),posstd*randn(3,Nelem-1)];
ActualElementPositions = NominalElementPositions + PositionPert;

Perturb the taper in magnitude and phase. Do not perturb the first taper.

NominalTaper = getTaper(sURA);
tapermagstd = 0.1;
taperphasestd = 0.1;
tapermagpert = tapermagstd*[0; randn(Nelem-1,1)];
ActualTaper = NominalTaper + tapermagpert;
taperphasepert = taperphasestd*[0;randn(Nelem-1,1)];
ActualTaper = ActualTaper.*exp(1i*taperphasepert);

Generate the pilot signals using the perturbed positions and tapers.

Nsamp  = 10000;
ncov = 0.001;
PilotAng = [-60,10,40,120,170; 0,-40,0,45,50];
Npilot = size(PilotAng,2);
for n = 1:Npilot
    X(:,:,n) = sensorsig(ActualElementPositions,Nsamp,...
        PilotAng(:,n),ncov,'Taper',ActualTaper);
end

Perform the pilot calibration

[estpos,esttaper] = pilotcalib(NominalElementPositions - ReferenceElement*ones(1,Nelem),...
    X,PilotAng,NominalTaper);

Add back the position of the reference sensor.

estpos = estpos + NominalElementPositions(:,1)*ones(1,Nelem);

Examine the root mean square (RMS) error of the calibrated parameters

Compute the RMS values of the initial taper perturbations to compare with the RMS values of the calibrated parameters.

tapermagpertRMSE = sqrt(tapermagpert'*tapermagpert/Nelem);
taperphasepertRMSE = sqrt(taperphasepert'*taperphasepert/Nelem);

Compute the RMS value of the calibrated taper magnitude error.

diff = abs(ActualTaper) - abs(esttaper);
diff2 = diff'*diff;
tapermagsolvRMSE = sqrt(diff2/Nelem);

Compare the calibrated RMS magnitude error to the initial RMS error. The calibration reduces the RMS magnitude error.

disp(tapermagsolvRMSE/tapermagpertRMSE)

    0.0014

Compute the RMS value of the calibrated taper phase error.

diff = unwrap(angle(ActualTaper) - angle(esttaper));
diff2 = diff'*diff;
tapersolvphaseRMSE = sqrt(diff2/Nelem);

Compare the calibrated RMS phase error to the initial RMS error. The calibration reduces the RMS phase error.

disp(tapersolvphaseRMSE/taperphasepertRMSE)

    0.0015

Compute the RMS value of the initial position error.

numpos = 3*Nelem;
initposRMSE = sqrt(sum(PositionPert(:).^2)/numpos);

Compute the RMS value of the calibrated position error.

solvposErr = ActualElementPositions - estpos;
solvposRMSE = sqrt(sum(solvposErr(:).^2)/(numpos));

Compare the calibrated RMS position error to initial position RMS error. The calibration reduces the RMS position error.

disp(solvposRMSE/initposRMSE)

   7.1582e-04

Estimate Selected ULA Parameters Using Pilot Calibration

Construct a 6-element ULA of isotropic antenna elements that are spaced one-half wavelength apart. Assume the array has been geometrically perturbed in the x-y plane and contains an unknown taper error. Perform pilot calibration on the array using four pilot sources at azimuth and elevation angles of (-60,0), (10,-40), (40,0), and (120,45) degrees. For the calibration process, pilot signals have an SNR of 30 dB. Each pilot signal contains 10,000 samples. Assume the signals have a frequency of 600 MHz.

Set up the ULA with nominal parameters

fc = 600e6;
c = physconst('LightSpeed');
lam = c/fc;
d = 0.5*lam;
sIso = phased.IsotropicAntennaElement('FrequencyRange',[100,900]*1e6);
Nelem = 6;
NominalTaper = ones(1,Nelem);
sULA = phased.ULA('Element',sIso,'NumElements',Nelem,'ElementSpacing',d,...
    'Taper',NominalTaper);

Create the pilot signals

Randomly perturb the element positions using a Gaussian distribution that has a standard deviation of 0.13 wavelength. Do not perturb the position of the first element.

posstd = 0.13;
rng default
NominalElementPositions = getElementPosition(sULA)/lam;
ReferenceElement = NominalElementPositions(:,1);
PositionPert = [zeros(3,1),posstd*randn(3,Nelem-1)];
ActualElementPositions = NominalElementPositions + PositionPert;

Perturb the taper in magnitude and phase. Do not perturb the first taper.

tapermagstd = 0.15;
taperphasestd = 0.15;
tapermagpert = tapermagstd*[0; randn(Nelem-1,1)];
ActualTaper = NominalTaper' + tapermagpert;
taperphasepert = taperphasestd*[0;randn(Nelem-1,1)];
ActualTaper = ActualTaper.*exp(1i*taperphasepert);

Generate the signals using the perturbed positions and tapers.

Nsamp  = 10000;
ncov = 0.001;
PilotAng = [-60,10,40,120; 0,-40,0,45];
Npilot = size(PilotAng,2);
for n = 1:Npilot
    X(:,:,n) = sensorsig(ActualElementPositions,Nsamp,...
        PilotAng(:,n),ncov,'Taper',ActualTaper);
end

Perform the pilot calibration

Turn off estimation of taper weights.

[estpos,esttaper] = pilotcalib(NominalElementPositions - ReferenceElement*ones(1,Nelem),...
    X,PilotAng,NominalTaper.',[1,1,1,0]');

Add back the position of the reference sensor

estpos = estpos + NominalElementPositions(:,1)*ones(1,Nelem);

Examine the root mean square (RMS) error of the calibrated parameters

Compute the RMS values of the initial taper perturbations to compare with the RMS values of the calibrated parameters.

tapermagpertRMSE = sqrt(tapermagpert'*tapermagpert/Nelem);
taperphasepertRMSE = sqrt(taperphasepert'*taperphasepert/Nelem);

Compute the RMS value of the calibrated taper magnitude error.

diff = abs(ActualTaper) - abs(esttaper);
diff2 = diff'*diff;
tapermagsolvRMSE = sqrt(diff2/Nelem);

Compare the calibrated RMS magnitude error to the initial RMS error. The calibration reduces the RMS magnitude error.

disp(tapermagsolvRMSE/tapermagpertRMSE)

    1.0000

Compute the RMS value of the calibrated taper phase error

diff = unwrap(angle(ActualTaper) - angle(esttaper));
diff2 = diff'*diff;
tapersolvphaseRMSE = sqrt(diff2/Nelem);

Compare the calibrated RMS phase error to the initial RMS error. The calibration reduces the RMS phase error.

disp(tapersolvphaseRMSE/taperphasepertRMSE)

     1

Compute the RMS value of the initial position error.

numpos = 3*Nelem;
initposRMSE = sqrt(sum(PositionPert(:).^2)/numpos);

Compute the RMS value of the calibrated position error.

solvposErr = ActualElementPositions - estpos;
solvposRMSE = sqrt(sum(solvposErr(:).^2)/(numpos));

Compare the calibrated RMS position error to initial position RMS error. The calibration reduces the RMS position error.

disp(solvposRMSE/initposRMSE)

    0.1502

Input Arguments

nompos — Nominal relative element positions

real-valued 3-by-N matrix

Nominal relative element positions, specified as a real-valued 3-by-N matrix. The dimension N is the number of elements in the sensor array. Elements positions are relative to the first element of the array and are specified in units of signal wavelength. Each column of nompos represents the [x;y;z] coordinates of the corresponding element. The nominal position of all sensors must be within one-half of a wavelength of their actual positions for successful calibration.

Data Types: double

x — Pilot signals

complex-valued L-by-N-by-M matrix

Pilot signals, specified as a complex-valued L-by-N-by-M matrix. The argument x represents the signals received by the array when pilot sources are transmitting. The dimension L is the number of snapshots of each pilot source signal. The dimension N is the number of array elements. The dimension M is the number of pilot sources.

Data Types: double
Complex Number Support: Yes

pilotang — Pilot angles

real-valued 2-by-M matrix

Pilot angles, specified as a real-valued 2-by-M matrix. The dimension M is the number of pilot sources. Each column contains the direction of the pilot source in the form [azimuth; elevation]. Angle units are in degrees. The azimuth angle must lie between -180° and 180° and the elevation angle must lie between -90° and 90°. The azimuth angle is measured from the x-axis to the projection of the source direction into the xy plane, positive toward the y-axis. The elevation angle is defined as the angle from the xy plane to the source direction, positive toward the z-axis. Calibration source directions must span sufficiently diverse azimuth and elevation angles.

Data Types: double

nomtaper — Nominal taper

1 (default) | complex-valued N-by-1 column vector

Nominal taper of array elements, specified as a complex-valued N-by-1 column vector. The dimension N is the number of array elements. Each component represents the nominal taper of the corresponding element.

Data Types: double
Complex Number Support: Yes

uncerts — Uncertainty estimation configuration

[1,1,1,1] (default) | 1-by-4 vector of ones and zeros

Uncertainty estimation configuration, specified as a 1-by-4 vector consisting of 0’s and 1’s. The vector uncerts determines which uncertainties to estimated. The vector takes the form of [xflag; yflag; zflag; taperflag]. Set xflag, yflag, or zflag to 1 to estimate uncertainties in the x, y, or z axes. Set taperflag to 1 to estimate uncertainties in the taper. The number of pilot sources must greater than or equal to the number of 1’s in the vector.

For example, set uncerts to [0;1;1;1] to estimate uncertainties in the y and z element position components and the taper simultaneously.

Data Types: double

Output Arguments

estpos — Estimated positions

real-valued 3-by-N matrix

Estimated element positions, returned as a real-valued 3-by-N matrix. Units are in signal wavelength. The dimension N is the number of array elements. Each column of estpos represents the [x;y;z] coordinates of the corresponding element.

esttaper — Estimated taper

complex-valued N-by-1 column vector

Estimated taper values, returned as a complex-valued N-by-1 column vector. The dimension N is the number of array elements. Each element of esttaper represents the taper of the corresponding sensor element.

Algorithms

This algorithm requires that the pilot sources be independent narrowband plane-wave signals incoming from the far field region of the array. In addition, signals must not exhibit multipath propagation effects or coherence. All elements in the sensor array are assumed to be isotropic.

The algorithm calibrates relative positions of the array sensors with respect to the first sensor. To use the algorithm, first subtract the position of the first element from each element, then pass the relative array into the function as the nominal position argument to produced the calibrated relative positions. Finally, add back the first element position to all the relative positions to create the fully calibrated array.

References

[1] N. Fistas and A. Manikas, "A New General Global Array Calibration Method", IEEE Proceedings of ICASSP, Vol. IV, pp. 73-76, April 1994.

Extended Capabilities

C/C++ Code Generation

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

Usage notes and limitations:

Does not support variable-size inputs.

Introduced in R2015a