FMCW Radar Altimeter Simulation
This example shows how to model a radar altimeter and measure the performance by simulating two scenarios using a land surface and moving platform. You will see how to build a waveform level altimeter model using a radarTransceiver and custom signal processor class. You will also see how to build a radarScenario by generating a LandSurface from an imported elevation data file and a Platform with an imported flight trajectory. Finally, you will see how to simulate the radar scenario to measure the performance of the radar altimeter.
Radar Parameter Definition
Aircraft have radar altimeters to determine the altitude above the local terrain as depicted in the figure below. Radar altimeters typically operate at 4.3 GHz and employ either FMCW or pulsed waveforms. These radars typically have very wide beamwidths (40 - 60 degrees) to provide accurate altitude estimates as the aircraft pitch and roll angles fluctuate throughout the flight.
The altimeter model uses parameters from the first altimeter described in Table 2 of the Radar Altimeter Recommendation ITU Document [1].
fCenter = 4.3e9; % Center frequency (Hz) bandwidth = 150e6; % Bandwidth (Hz) prf = 143; % Pulse repetition frequency (Hz) antBeamwidth = 40; % Beamwidth (degrees) pTransmitter = 0.4; % Transmitter power (W) rMax = 1676; % Highest altitude required (m) sweepDir = 'Triangle'; % Sweep direction of the FMCW waveform noiseFigure = 8; % Receiver noise figure (dB)
Creating the Radar Altimeter Model
The radar altimeter model consists of two components. The first component is a radarTransceiver. This component models the transmission of the waveform and reception of the ground echo at the altimeter. The second component is a custom component created for this example. This component takes the IQ samples generated by the radarTransceiver and performs the signal processing necessary to calculate the altitude.
Antenna Design
A phased.GaussianAntennaElement is used to model the altimeter antenna element. The idealized antenna pattern of this antenna element has no sidelobes but will suffice for the purposes of this simulation.
% Generate the Gaussian antenna element frequencyRange = [fCenter-bandwidth/2 fCenter+bandwidth/2]; altimeterAntennaElement = phased.GaussianAntennaElement("Beamwidth",antBeamwidth,"FrequencyRange",frequencyRange); % Plot the beamwidth to verify that it matches the expected antenna beamwidth beamwidth(altimeterAntennaElement,fCenter);
The transmit and receive antennas are modeled as a phased.Radiator and phased.Collector. Each antenna contains the phased.GaussianAntennaElement.
altimeterTxAntenna = phased.Radiator('OperatingFrequency',fCenter,'Sensor',altimeterAntennaElement, 'SensorGainMeasure','dBi'); altimeterRxAntenna = phased.Collector('OperatingFrequency',fCenter,'Sensor',altimeterAntennaElement,'SensorGainMeasure','dBi');
Waveform Design
Altimeters commonly use FMCW technology for many reasons, including the ability to measure short distances, high range accuracy, and low power requirements as discussed in [1]. Therefore, this example uses a triangle FMCW waveform.
A key processing step for FMCW radar is dechirping. Dechirping is the process of mixing the transmitted waveform with the received signal. For dechirping to work properly in simulation, the sampling frequency of the waveform must be at least the signal bandwidth. The simulation differs from the real world in this regard, because in a real receiver, dechirping occurs in the analog domain. The A/D sampling rate in an FMCW receiver only needs to be high enough to capture the bandwidth of the dechirped signal, which is much lower than the signal bandwidth.
Therefore, two sampling rates are used for the purpose of this simulation. The first sampling rate called the waveformSamplingRate is based on the waveform bandwidth and is used to generate the IQ and dechirp the signal. The second sampling rate called the receiverSamplingRate is based on the maximum expected beat frequency in the dechirped signal. After dechirping, the signal is down-sampled to the receiverSamplingRate to perform the remainder of the signal processing tasks.
To generate the waveform, first calculate the sweep slope and set the waveformSamplingRate, which is based on the bandwidth. The sampling rate of the waveform is not exactly the bandwidth, because the waveform sampling rate times the sweep time must be an integer for the simulation to work properly.
tSweep = 1/(2*prf); % Sweep time (s) is half of the sweep period because this is a two sided FMCW waveform nSweep = 2; % There are two sweeps in one processing interval, 1 upsweep, 1 downsweep sweepSlope = bandwidth/tSweep; % Sweep slope (Hz / s) waveformSamplingRate = round(bandwidth*tSweep)/tSweep; % Sampling rate of the waveform (Hz)
Create the FMCW waveform once the sweep time and waveform sampling rate have been calculated.
altimeterWaveform = phased.FMCWWaveform('SweepBandwidth',bandwidth,'SampleRate',waveformSamplingRate,'SweepDirection',sweepDir,'SweepTime',tSweep,'OutputFormat','Sweeps');
Create the Radar Transceiver
The radarTransceiver can be fully modeled now that the antenna and waveform have been created.
Create the transmitter and receiver using previously defined system parameters and the selected sampling rate calculated in the previous section.
altimeterTransmitter = phased.Transmitter('PeakPower',pTransmitter); altimeterReceiver = phased.ReceiverPreamp('SampleRate', waveformSamplingRate,'NoiseFigure', noiseFigure);
Create the radarTransceiver. The mounting angle is set so that the altimeter is facing towards the surface in the North, East, Down (NED) reference frame.
altimeterTransceiver = radarTransceiver('Waveform',altimeterWaveform,... 'Transmitter',altimeterTransmitter, ... 'TransmitAntenna',altimeterTxAntenna, ... 'ReceiveAntenna',altimeterRxAntenna, ... 'Receiver',altimeterReceiver, ... 'MountingAngles', [0 -90 0], ... 'NumRepetitions', nSweep)
altimeterTransceiver =
radarTransceiver with properties:
Waveform: [1x1 phased.FMCWWaveform]
Transmitter: [1x1 phased.Transmitter]
TransmitAntenna: [1x1 phased.Radiator]
ReceiveAntenna: [1x1 phased.Collector]
Receiver: [1x1 phased.ReceiverPreamp]
MechanicalScanMode: 'None'
ElectronicScanMode: 'None'
MountingLocation: [0 0 0]
MountingAngles: [0 -90 0]
NumRepetitionsSource: 'Property'
NumRepetitions: 2
RangeLimitsSource: 'Property'
RangeLimits: [0 Inf]
RangeOutputPort: false
TimeOutputPort: false
Create the Signal Processor
Collect the upsweep and downsweep waveform samples from the altimeterTransceiver. The signal processor uses these samples to dechirp the received signals.
upsweepWaveformSamples = altimeterTransceiver.Waveform(); % Upsloping signal IQ downsweepWaveformSamples = altimeterTransceiver.Waveform(); % Downsloping signal IQ
Plot the spectrogram of the waveform to ensure it matches expectations for a triangular FMCW waveform.
figure; stft([upsweepWaveformSamples;downsweepWaveformSamples],waveformSamplingRate,'FrequencyRange','twosided') title('Spectrogram of Altimeter Waveform (BW = 150 MHz, PRF = 143 Hz)')
The IQ data is down-sampled to the receiverSamplingRate by the signal processor after dechirping. The receiverSamplingRate must be at least twice the maximum beat frequency. The maximum beat frequency is the sum of the maximum range delay frequency and the maximum Doppler shift frequency. Calculate the maximum beat frequency and receiver sampling rate.
frmax = range2beat(rMax,sweepSlope); % Maximum frequency (Hz) due to range vmax = 8; % Max rate of descent (m/s) h = freq2wavelen(fCenter); % Wavelength (m) fdmax = 2*speed2dop(vmax,h); % Maximum frequency shift (Hz) due to radial velocity fbmax = frmax+fdmax; % Maximum beat frequency (Hz) receiverSamplingRate = 2*fbmax; % Sampling rate of the receiver (Hz)
A custom class represents the altimeter signal processor for this example. Now that all required inputs have been generated, create the altimeterSignalProcessor. The concept of this signal processor is discussed more in a later section.
altimeterSignalProcessor = helperAltimeterSignalProcessor(upsweepWaveformSamples,downsweepWaveformSamples,sweepSlope,waveformSamplingRate,receiverSamplingRate)
altimeterSignalProcessor =
helperAltimeterSignalProcessor with properties:
UpsweepReference: [524476x1 double]
DownsweepReference: [524476x1 double]
Detector: [1x1 phased.CFARDetector]
RangeResponse: [1x1 phased.RangeResponse]
Create Radar Scenario with a Land Surface
The simulation of the altimeter model is composed of an altimeter mounted on a platform with a flight trajectory over a land surface. The first part of the scenario setup is creating a LandSurface from imported elevation data. This elevation data was generated from a Digital Terrain Elevation Data (DTED) file that contains real world terrain data from Chicago O'Hare Airport. For more information on generating surface clutter signals using real terrain data, see Introduction to Radar Scenario Clutter Simulation.
Initialize the Radar Scenario
Initialize the radarScenario before adding a LandSurface and ClutterGenerator. tStart and tStop refer to the start and stop time for the trajectory. tUpdate controls the size of the time step used by the simulation. This value is set to 5 seconds instead of the default update rate to speed up the simulation. The default update rate if unset would be the PRF of the radar.
tStart = 300; % Simulated trajectory start (s) tStop = 395; % Simulated trajectory stop (s) tUpdate = 5; % Simulated trajectory update interval (s) scene = radarScenario('StopTime',tStop-tStart,'UpdateRate',1/tUpdate,'IsEarthCentered',true);
Generate the Surface and Clutter
Elevation data is loaded to generate the surface in the boundary of interest. The reflectivity of the land is created using an APL model that supports a 90 degree grazing angle in an urban environment. The coordinates of the flight trajectory are used to choose the boundaries of the surface. A LandSurface and ClutterGenerator are added to the radarScenario to generate IQ data of the ground reflection as the simulation progresses.
load('surfaceHeightData.mat','surfaceHeight') % Load surface elevation data refl = surfaceReflectivityLand(Model = "APL", LandType = 'Urban'); % Land reflectivity model boundaryLatLim = [41.9799 41.9879]; % Set the boundary latitude limits (deg) boundaryLonLim = [-87.9000 -87.7750]; % Set the boundary longitude limits (deg) % Create the land surface and clutter generator surface = landSurface(scene,'Terrain',surfaceHeight.','Boundary',[boundaryLatLim; boundaryLonLim],'RadarReflectivity',refl); clutter = clutterGenerator(scene,altimeterTransceiver);
Landing Scenario Simulation
Altimeter performance requirements vary as a function of altitude and altimeters must continue to provide accurate altitude information as the aircraft roll angle changes. To verify that the altimeter meets these requirements, a landing scenario and a banking scenario are simulated. The landing scenario is discussed first followed by the banking scenario.
Generate Flight Trajectory
The landing scenario uses a real flight trajectory of an aircraft approaching Chicago O'Hare Airport. The last step in setting up the simulation is adding a platform to the radarScenario with the altimeter radarTransceiver as a sensor and a trajectory using the landing trajectory.
flightTrajectoryFile = 'ORD_FlightTraj.mat'; landingFlightTrajectory = helperGetFlightTrajectory(surface,flightTrajectoryFile,tStart,tStop,tUpdate); altimeterPlatform = platform(scene,'Sensors',altimeterTransceiver,'Trajectory',landingFlightTrajectory,'Dimensions',struct('Length',.00075,'Width',20,'Height',.007,'OriginOffset',[0 0 0]));
The animations below visualize the landing scenario. Execute the following line of code to generate the animations.
helperAltimeterGlobeAnimation(altimeterTransceiver,tStop,tStart,tUpdate,boundaryLatLim,boundaryLonLim,surfaceHeight,landingFlightTrajectory,'Landing')
Run Simulation
First set the random number generator seed for reproducible results and setup containers for storing data as the simulation proceeds.
rng("default"); % Set random number generator seed nSimulationPoints = floor((tStop-tStart)/tUpdate)+1; % Calculate the number of simulation points landingRadarSignal = cell(1,nSimulationPoints); % Store the raw radar signal landingMeasuredAltitude = zeros(1,nSimulationPoints); % Store measured altitude landingTruthAltitude = zeros(1,nSimulationPoints); % Store truth altitde landingRollAngle = zeros(1,nSimulationPoints); % Store the aircraft roll angle in degrees cnt = 1; % Initialize count variable
With the variables set up, run the simulation calculating the measured altitude, true altitude, and roll angle at each step. The clutter resolution is set at each simulation step so that the density of clutter scatterers is sufficient. A higher scatterer density is required as the aircraft descends.
while advance(scene) % Calculate truth altitude [landingTruthAltitude(cnt), landingRollAngle(cnt)] = helperCalculateGroundTruth(surface,altimeterPlatform); % Calculate required clutter resolution and set accordingly clutter.Resolution = helperDetermineClutterResolution(landingTruthAltitude(cnt)); % Receive the IQ and track aircraft position and simulation time iqsig = receive(scene); radarSignal = iqsig{:}; % save the first step of iq data for later analysis if cnt == 1 firstLandingSignal = radarSignal; end % Calculate measured and truth altitude landingMeasuredAltitude(cnt) = altimeterSignalProcessor.estimateAircraftAltitude(radarSignal); cnt = cnt + 1; end % Store the last signal for later analysis lastLandingSignal = radarSignal;
Analyze Altimeter Signal Processor Ground Detection Algorithm
FMCW radars determine the range and velocity of a target by dechirping the received signal with the transmitted signal and measuring the beat frequency of the dechirped signal. For more information on using FMCW radars for target range estimation, see Automotive Adaptive Cruise Control Using FMCW Technology. This section discusses the altimeter signal processor functionality before looking at the overall performance of the altimeter during the landing scenario. The helperAltimeterSignalProcessor class contains the logic for the signal processor.
Use the data from the first step of the landing scenario simulation to analyze the ground detection algorithm.
trueAltitude = landingTruthAltitude(1); % Get the true altitude based on aircraft positionThe first step in the signal processing algorithm is dechirping the received signal by mixing it with the known transmitted signal. The frequency of the dechirped signal is directly proportional to the range of the target, in this case the ground. After the signal is dechirped, the signal processor converts frequency values to range values using the known sweep slope. This processing step is handled by the phased.RangeResponse object in the altimeterSignalProcesor. The following plot illustrates the range response for the upsweep and downsweep signals that is the magnitude of the received signal as a function of range.
plotStruct = struct('TrueAltitude',trueAltitude,'XLims',[trueAltitude-50 trueAltitude+50]); altimeterSignalProcessor.plotRangeResponse(firstLandingSignal,plotStruct);
The plots above show that the returned energy from an altimeter at this altitude are spread over about 30 meters. This is expected because the target is the ground, which has a large RCS and the altimeter uses a large beamwidth.
A CFAR algorithm is used to generate detections using the smallest of cell averaging technique to select for detections on the edge of the clutter return, as this is where the ground lies. The number of detections is limited by using a false alarm rate of 1e-7. The following plots show the range response of the signals in addition to the detections that the CFAR algorithm in the signal processor found.
plotStruct = struct('TrueAltitude',trueAltitude,'Detections',[],'XLims',[trueAltitude-50 trueAltitude+50]); altimeterSignalProcessor.plotRangeResponse(firstLandingSignal,plotStruct);
Once the detection algorithm has been applied, the detection with the highest power is used to estimate the aircraft altitude. The following plot shows the ranges that are selected to estimate the altitude from the upsweep and downsweep signals.
plotStruct = struct('TrueAltitude',trueAltitude,'Detections',[],'SelectedRange',[],'XLims',[trueAltitude-50 trueAltitude+50]); altimeterSignalProcessor.plotRangeResponse(firstLandingSignal,plotStruct);
The effect of Doppler shift on the altitude estimate is removed by averaging the upsweep and downsweep altitude. The effectiveness of the altimeter algorithm can be checked by comparing the measured altitude to the ground truth altitude. For the first simulation step, the estimate is less than 5 meters from the true distance to the ground, which is an error of about 1%.
altitudeMeasured = altimeterSignalProcessor.estimateAircraftAltitude(firstLandingSignal);
fprintf('Measured Altitude = %0.1f m, True Altitude = %0.1f m, Altitude Error = %0.1f m',altitudeMeasured,trueAltitude,trueAltitude - altitudeMeasured);Measured Altitude = 435.6 m, True Altitude = 429.0 m, Altitude Error = -6.6 m
Landing Scenario Performance
The landing scenario measures the ability of the altimeter model to continue to produce good results as the range spread and power level of the clutter return vary. At higher altitudes, the return occurs over a much wider range and lower power level than at lower altitudes.
The plot below shows the range response of the highest and lowest altitude data point. Notice the difference in power as well as the width of the range response. The landing scenario tests the ability of the altimeter to meet performance requirements as the returned signal varies dramatically in shape and power.
plotStruct = struct('XLims',[0 600],'LegendName', ['Altitude = ',num2str(landingTruthAltitude(1),'%.1f')],'PlotDb',[]); [upaxis, downaxis] = altimeterSignalProcessor.plotRangeResponse(firstLandingSignal,plotStruct); plotStruct = struct('XLims',[0 600],'LegendName', ['Altitude = ',num2str(landingTruthAltitude(end),'%.1f')],'UpDownAxis',[upaxis,downaxis],'PlotDb',[]); altimeterSignalProcessor.plotRangeResponse(lastLandingSignal,plotStruct);
The following two plots visualize altimeter performance as the altitude of the aircraft changes. The first plot shows the altimeter altitude and true altitude throughout the duration of the landing scenario superimposed on top of the accuracy requirements derived from [1]. The second plot shows the altitude error throughout the landing scenario compared to the accuracy requirements derived from [1].
simTime = 0:tUpdate:(tStop-tStart); helperPlotAltimeterAccuracy(landingTruthAltitude, landingMeasuredAltitude, simTime, 'Time (s)', 'Measured vs. True Altitude: Landing Trajectory')
helperPlotAltimeterError(landingTruthAltitude,landingMeasuredAltitude,landingTruthAltitude,'True Altitude (m)','Measurement Error: Landing Trajectory',true)
In general, the altimeter as modeled meets the performance requirements laid out in [1] for the landing scenario. The error appears to reduce as the aircraft altitude drops. This makes sense, because the beam footprint gets smaller as the aircraft altitude drops, so the clutter returns happen over a smaller area. Note that the performance indicates that the algorithm could be improved.
Banking Scenario Simulation
The radar altimeter must continue to meet performance requirements as the aircraft roll and pitch angles change. A scenario in which the aircraft roll angle varies is simulated to evaluate this capability.
Generate Flight Trajectory and Update Radar Scenario
A helper function generates the flight trajectory for the aircraft banking scenario by modifying the landing trajectory to have a constant altitude and a roll angle that varies between positive and negative 1/2 beamwidth of the antenna. In this scenario, the flight direction remains constant despite the variation in aircraft roll for the sake of simplicity.
bankingFlightTrajectory = helperGetFlightTrajectory(surface,flightTrajectoryFile,tStart,tStop,tUpdate,antBeamwidth / 2);
With the banking trajectory created, the scenario is restarted, and the flight trajectory of the platform is updated to use the banking trajectory instead of the landing trajectory.
altimeterPlatform.Trajectory = bankingFlightTrajectory;
The following animations visualize the banking scenario. Execute the following line of code to generate the animations.
helperAltimeterGlobeAnimation(altimeterTransceiver,tStop,tStart,10,boundaryLatLim,boundaryLonLim,surfaceHeight,bankingFlightTrajectory,'Banking')
Run Simulation
Run the simulation now that the platform in the scenario follows a banking flight trajectory. In this case a helper function runs the simulation because the exact same block of code was used to simulate the landing scenario.
signalToSave = [1 ceil(nSimulationPoints/2)]; % Save the first radar signal and the middle radar signal for later analysis
[bankingMeasuredAltitude,bankingTruthAltitude,bankingRollAngle,bankingRadarSignal,maxErrorSignal,maxErrorIdx] = helperSimulateAltimeterScenario(scene,nSimulationPoints,altimeterSignalProcessor,surface,clutter,signalToSave);Banking Scenario Performance
The banking scenario measures the ability of the altimeter model to continue to produce good results as the shape of the clutter return changes due to aircraft roll angle. When the aircraft has a high roll angle, the nadir is at the edge of the beam footprint and the beam footprint covers a greater area. When the aircraft has a low roll angle, the nadir is in the center of the beam footprint and the beam footprint covers a smaller area.
The following plots show the range response when the roll angle of the aircraft is at its highest (20 degrees) vs. lowest (1 degree). Notice that when the roll angle is lower, the maximum range response power is greater, and the spread of ranges is more narrow compared to the highest roll angle case.
plotStruct = struct('XLims',[400 600],'LegendName', ['Roll Angle = ',num2str(bankingRollAngle(1),'%.1f')],'PlotDb',[]); [upaxis, downaxis] = altimeterSignalProcessor.plotRangeResponse(bankingRadarSignal{1},plotStruct); plotStruct = struct('XLims',[400 600],'LegendName', ['Roll Angle = ',num2str(bankingRollAngle(ceil(nSimulationPoints/2)),'%.1f')],'UpDownAxis',[upaxis,downaxis],'PlotDb',[]); altimeterSignalProcessor.plotRangeResponse(bankingRadarSignal{end},plotStruct);
The following two plots visualize altimeter performance as the roll angle of the aircraft changes. Once again plot the measured and true altitude on the same graph along with the acceptable error limits for altimeters derived from [1]. Additionally, plot the measurement error along with the acceptable error limits. Roll angle is the independent variable in this case.
helperPlotAltimeterAccuracy(bankingTruthAltitude, bankingMeasuredAltitude, bankingRollAngle, 'Roll Angle (deg)', 'Measured vs. True Altitude: Banking Trajectory');
helperPlotAltimeterError(bankingTruthAltitude,bankingMeasuredAltitude,bankingRollAngle,'Roll Angle (deg)','Measurement Error: Banking Trajectory',false);
The performance of the altimeter in the banking scenario meets the requirements described in [1] as aircraft rolls and the shape of the signal return changes. There are a couple of points that appear to be outliers in terms of error magnitude, which are discussed in the next section.
Results Discussion
In general, the altimeter meets the performance requirements described in [1] in the simulations of the landing and banking scenarios. In both cases, there are points with particularly high measurement errors, suggesting that improvements could still be made to this altimeter design.
As discussed, the ground is detected using a smallest of cell averaging CFAR algorithm to select for points on the edge of the clutter return. The highest power detection is used to estimate the aircraft altitude. The larger errors occur when the selected detection is not made on the leading edge of the clutter return. The following plot shows the clutter return with detections in addition to the CFAR threshold from the point with the largest error in the banking scenario.
maxErrorAltitude = bankingTruthAltitude(maxErrorIdx); % Return signal for the maximum measurement error point % Plot the detections for the maximum measurement error point in the % landing scenario plotStruct = struct('TrueAltitude',maxErrorAltitude,'SelectedRange',[],'Detections',[],'Threshold',[],'XLims',[maxErrorAltitude-10 maxErrorAltitude+30]); altimeterSignalProcessor.plotRangeResponse(maxErrorSignal,plotStruct);
In the plot above, there is a higher power detection that occurs on the trailing edge of the clutter return. This leads to a relatively high measurement error. A similar pattern is seen in other points with large errors This suggests that some further improvements could be made to the altimeter design by fine-tuning the ground detection algorithm.
Summary
This example models a simple radar altimeter and uses surface clutter as a target to determine aircraft altitude. The example starts by building the radar altimeter mode from a transceiver that generates IQ data and a signal processor that generates the ground range estimate from the IQ data. Then, two radar scenarios are simulated and performance of the altimeter is measured in both scenarios.
References
[1] Operational and Technical Characteristics and Protection Criteria of Radio Altimeters Utilizing the Band 4200-4400 MHz. Recommendation, M.2059-0, ITU-R, Feb. 2014, https://www.itu.int/dms_pubrec/itu-r/rec/m/R-REC-M.2059-0-201402-I!!PDF-E.pdf.
Supporting Functions
The helperAltimeterSignalProcessor class is shown below.
type('helperAltimeterSignalProcessor.m')classdef helperAltimeterSignalProcessor
% This class is used in the Altimeter Simulation example to perform the
% signal processing for the radar altimeter.
properties
UpsweepReference
DownsweepReference
Detector
RangeResponse
end
methods
function obj = helperAltimeterSignalProcessor(refSigUp,refSigDown,sweepSlope,fsWaveform,fsReceiver)
obj.UpsweepReference = refSigUp;
obj.DownsweepReference = refSigDown;
decimationFactor = floor(fsWaveform / fsReceiver);
obj.Detector = phased.CFARDetector('Method','SOCA','NumTrainingCells',50,'NumGuardCells',2,'ProbabilityFalseAlarm',1e-7,"ThresholdOutputPort",true);
obj.RangeResponse = phased.RangeResponse("RangeMethod","FFT",...
"SampleRate",fsWaveform,"SweepSlope",sweepSlope,"DechirpInput",true,"DecimationFactor",decimationFactor);
end
function [altitude,rngup,rngdn,respup,respdn] = estimateAircraftAltitude(obj,radarSignal)
% Apply CFAR detector to return signal
[grndrangesup,grndrangesdn,grndmagrespup,grndmagrespdn] = detectGround(obj,radarSignal);
% Find the maximum power detection for upsweep and downsweep
[rngup,respup] = findMaxMagResp(grndmagrespup,grndrangesup);
[rngdn,respdn] = findMaxMagResp(grndmagrespdn,grndrangesdn);
% Estimate altitude as average of up and down sweeps
if isempty(rngup)
rngup = -rngdn;
end
if isempty(rngdn)
rngdn = -rngup;
end
altitude = (-rngdn + rngup)/2;
function [rng,rngrspmax] = findMaxMagResp(rngrsp,fullrng)
[rngrspmax,rngidx] = max(abs(rngrsp));
rng = fullrng(rngidx);
end
end
function [grndrangesup,grndrangesdn,grndmagrespup,grndmagrespdn,range,upthreshold,downthreshold] = detectGround(obj,radarSignal)
% Get the range response of the radar signal
[rngrespup,rngrespdn,range] = obj.getRangeResponse(radarSignal);
% Detect ground ranges above threshold from upsweep
[upidx,upthreshold] = obj.Detector(abs(rngrespup).^2,1:length(rngrespup));
grndrangesup = range(upidx);
grndmagrespup = rngrespup(upidx);
% Detect ground ranges above threshold from downsweep
[dnidx,downthreshold] = obj.Detector(abs(rngrespdn).^2,1:length(rngrespdn));
grndrangesdn = range(dnidx);
grndmagrespdn = rngrespdn(dnidx);
% Take square root to get thresholds because thresholding done
% on squares
upthreshold = sqrt(upthreshold);
downthreshold = sqrt(downthreshold);
end
function [rngrespup,rngrespdn,range] = getRangeResponse(obj,radarSignal)
% Get the range response of the radar signal
[rngrespup,~] = obj.RangeResponse(radarSignal(:,1),obj.UpsweepReference);
[rngrespdn,range] = obj.RangeResponse(radarSignal(:,2),obj.DownsweepReference);
end
%------------------------------------------------------------------
% Visualization Functions
%------------------------------------------------------------------
function [upsweepSubplotAxis,downsweepSubplotAxis] = plotRangeResponse(obj,radarSignal,plotStruct)
% Plot the dechirped signals in the frequency domain.
[rangeResponseUpsweep,rangeResponseDownsweep,range] = obj.getRangeResponse(radarSignal);
if isfield(plotStruct,'UpDownAxis')
upsweepSubplotAxis = plotStruct.UpDownAxis(1);
downsweepSubplotAxis = plotStruct.UpDownAxis(2);
else
figure();
upsweepSubplotAxis = subplot(2,1,1);
downsweepSubplotAxis = subplot(2,1,2);
end
if isfield(plotStruct,'LegendName')
rangeResponseName = plotStruct.LegendName;
showLegend = true;
else
rangeResponseName = 'Range Response';
showLegend = false;
end
if isfield(plotStruct,'PlotDb')
rangeResponseUpsweepDb = 20*log10(abs(rangeResponseUpsweep));
rangeResponseDownsweepDb = 20*log10(abs(rangeResponseDownsweep));
plotRangeResponse(upsweepSubplotAxis,range,rangeResponseUpsweepDb,'Upsweep Range Response','Power (dB)',rangeResponseName,showLegend,[0 max(range)]);
plotRangeResponse(downsweepSubplotAxis,-range,rangeResponseDownsweepDb,'Downsweep Range Response','Power (dB)',rangeResponseName,showLegend,[0 max(-range)]);
else
plotRangeResponse(upsweepSubplotAxis,range,abs(rangeResponseUpsweep),'Upsweep Range Response','Magnitude',rangeResponseName,showLegend,[0 max(range)]);
plotRangeResponse(downsweepSubplotAxis,-range,abs(rangeResponseDownsweep),'Downsweep Range Response','Magnitude',rangeResponseName,showLegend,[0 max(-range)]);
end
if isfield(plotStruct,'TrueAltitude')
truthAltitude = plotStruct.TrueAltitude;
obj.plotTruthAltitude(truthAltitude,upsweepSubplotAxis,downsweepSubplotAxis);
end
if isfield(plotStruct,'Detections')
obj.plotDetections(radarSignal,upsweepSubplotAxis,downsweepSubplotAxis);
end
if isfield(plotStruct,'SelectedRange')
obj.plotSelectedRange(radarSignal,upsweepSubplotAxis,downsweepSubplotAxis);
end
if isfield(plotStruct,'Threshold')
obj.plotThreshold(radarSignal,upsweepSubplotAxis,downsweepSubplotAxis);
end
if isfield(plotStruct,'XLims')
xlim(upsweepSubplotAxis,plotStruct.XLims)
xlim(downsweepSubplotAxis,plotStruct.XLims)
end
function plotRangeResponse(axis,range,mResponse,plotTitle,yLabel,rangeResponseName,showLegend,plotLims)
hold(axis,"on")
plot(axis,range,mResponse,'DisplayName',rangeResponseName,'LineWidth',1);
title(axis,plotTitle);
ylabel(axis,yLabel);
xlabel(axis,'Range (m)');
xlim(axis,plotLims);
if showLegend
legend(axis,'Location','northeast');
end
hold(axis,"off")
end
end
function [upsweepSubplotAxis,downsweepSubplotAxis] = plotTruthAltitude(~,truthAltitude,upsweepSubplotAxis,downsweepSubplotAxis)
% This function plots the truth altitude on top of the
% range response
plotTruthAltitude(upsweepSubplotAxis,truthAltitude);
plotTruthAltitude(downsweepSubplotAxis,truthAltitude);
function plotTruthAltitude(axis,truthFrequency)
hold(axis,"on")
ylimit = axis.YLim;
x = [truthFrequency truthFrequency];
y = ylimit;
plot(axis,x,y,':','DisplayName','True Altitude','LineWidth',2);
legend(axis,'Location','northeast');
hold(axis,"off")
end
end
function [upsweepSubplotAxis,downsweepSubplotAxis] = plotDetections(obj,radarSignal,upsweepSubplotAxis,downsweepSubplotAxis)
% Plot the CFAR detections on the axis provided
[rDetectionsUpsweep,rDetectionsDownsweep,mDetectionsUpsweep,mDetectionsDownsweep] = obj.detectGround(radarSignal);
plotDetections(obj,upsweepSubplotAxis,rDetectionsUpsweep,mDetectionsUpsweep);
plotDetections(obj,downsweepSubplotAxis,-rDetectionsDownsweep,mDetectionsDownsweep);
function plotDetections(obj,axis,rDetection,mDetection)
obj.plotScatterOnAxis(axis,rDetection,abs(mDetection),'Detections');
end
end
function [upsweepSubplotAxis,downsweepSubplotAxis] = plotThreshold(obj,radarSignal,upsweepSubplotAxis,downsweepSubplotAxis)
% Plot the CFAR threshold on the axis provided
[~,~,~,~,range,upthreshold,downthreshold] = obj.detectGround(radarSignal);
plotThresholdOnAxis(upsweepSubplotAxis,range,upthreshold);
plotThresholdOnAxis(downsweepSubplotAxis,-range,downthreshold);
function plotThresholdOnAxis(axis,range,threshold)
hold(axis,"on");
plot(axis,range,threshold,'DisplayName','CFAR Threshold');
legend(axis,'Location','northeast');
hold(axis,"off");
end
end
function [upsweepSubplotAxis,downsweepSubplotAxis] = plotSelectedRange(obj,radarSignal,upsweepSubplotAxis,downsweepSubplotAxis)
% Plot the selected range on the axis provided
[~,rSelectedUpsweep,rSelectedDownsweep,mSelectedUpsweep,mSelectedDownsweep] = obj.estimateAircraftAltitude(radarSignal);
plotSelectedRange(obj,upsweepSubplotAxis,rSelectedUpsweep,mSelectedUpsweep);
plotSelectedRange(obj,downsweepSubplotAxis,-rSelectedDownsweep,mSelectedDownsweep);
function plotSelectedRange(obj,axis,fSelected,mSelected)
obj.plotScatterOnAxis(axis,fSelected,abs(mSelected),'Selected Range');
end
end
function plotScatterOnAxis(~,axis,xdata,ydata,name)
% Plot a scatter plot on an existing axis without modifying the
% axis.
hold(axis,"on");
scatter(axis,xdata,ydata,'filled','DisplayName',name);
legend(axis,'Location','northeast');
hold(axis,'off')
end
end
end
Various other functions defined below are used throughout this example.
function [TruthAltitude, RollAngle] = helperCalculateGroundTruth(surface,platform) % Calculate the actual altitude and aircraft roll angle based on current aircraft position % Terrain elevation height at the current aircraft position ElevationHeight = surface.height([platform.Position(1);platform.Position(2)]); % Ground truth altitude TruthAltitude = platform.Position(3) - ElevationHeight; % Roll angle RollAngle = platform.Orientation(2); end function flightTrajectory = helperGetFlightTrajectory(surface,FlightTrajectoryFile,n1,n2,updateInterval,varargin) % Return a flight trajectory with elevation from trajectory data file % Load flight trajectory load(FlightTrajectoryFile,'flightTrajectory'); % Init arrival time ArrivalTime = 0:size(flightTrajectory,1)-1; % Add the ground elevation to flight elevation flightTrajectory = flightTrajectory(n1:updateInterval:n2,:); groundElevation = surface.height(flightTrajectory(:,1:2)')'; flightTrajectory(:,3) = flightTrajectory(:,3) + groundElevation; % Generate the flight trajectory ArrivalTime = ArrivalTime(n1:updateInterval:n2)-ArrivalTime(n1); stepNumber = numel(ArrivalTime); orientationEuler = zeros(stepNumber,3); % If provided, generate a fixed altitude trajectory with bank angle if nargin > 5 bankAngle = varargin{1}; allBankAngles = -bankAngle:2*bankAngle/(stepNumber-1):bankAngle; orientationEuler(:,2) = allBankAngles; orientationQuat = quaternion(orientationEuler,'rotvecd'); flightTrajectory(:,3) = flightTrajectory(:,3) - groundElevation; flightTrajectory(:,3) = flightTrajectory(1,3); flightTrajectory(:,3) = flightTrajectory(:,3) + groundElevation; flightTrajectory = geoTrajectory(flightTrajectory,ArrivalTime,'Orientation',orientationQuat); else orientationQuat = quaternion(orientationEuler,'rotvecd'); flightTrajectory = geoTrajectory(flightTrajectory,ArrivalTime,'Orientation',orientationQuat); end end function [MeasuredAltitude,TruthAltitude,RollAngle,RadarSignal,MaxErrorSignal,MaxErrorIdx] = helperSimulateAltimeterScenario(scene,nSimulationPoints,altimeterSignalProcessor,surface,clutter,radarSignalPoints) % Simulate the scenario and return key metrics % Setup scene rng(1); restart(scene); % Initialize variables RadarSignal = cell(1,length(radarSignalPoints)); MeasuredAltitude = zeros(1,nSimulationPoints); TruthAltitude = zeros(1,nSimulationPoints); RollAngle = zeros(1,nSimulationPoints); maxMeasurementError = -Inf; cnt = 1; % Run the simulation while advance(scene) % Calculate truth altitude [TruthAltitude(cnt), RollAngle(cnt)] = helperCalculateGroundTruth(surface,scene.Platforms{1}); % Calculate required clutter resolution and set accordingly clutter.Resolution = helperDetermineClutterResolution(TruthAltitude(cnt)); % Receive the iq and track aircraft position and simulation time iqsig = receive(scene); radarSignal = iqsig{:}; % Calculate measured altitude MeasuredAltitude(cnt) = altimeterSignalProcessor.estimateAircraftAltitude(radarSignal); measurementError = abs(MeasuredAltitude(cnt)-TruthAltitude(cnt)); if (measurementError > maxMeasurementError) maxMeasurementError = measurementError; MaxErrorSignal = radarSignal; MaxErrorIdx = cnt; end if any(radarSignalPoints == cnt) RadarSignal{radarSignalPoints == cnt} = radarSignal; end cnt = cnt + 1; end end function helperPlotAltimeterAccuracy(TruthAltitude,MeasuredAltitude,xValues,xAxisName,chartTitle) % Plot the actual altimeter accuracy compared to the required accuracy % Get the accuracy limits for each truth altitude accuracyLimits = helperGetAccuracyLimits(TruthAltitude); % Determine whether any of the measured altitudes were too high or too low measureErrorLow = MeasuredAltitude < accuracyLimits.lowerLimit; measureErrorHigh = MeasuredAltitude > accuracyLimits.upperLimit; measureErrorIdx = measureErrorLow | measureErrorHigh; figure; hold on; % Plot the accuracy limits plotAllAccuracyLimits(xValues,accuracyLimits) plot(xValues,MeasuredAltitude,'-o','DisplayName','Altimeter Altitude', ... 'color',[0 0.45 0.74]) plot(xValues,TruthAltitude,'-s','DisplayName','Truth Altitude', ... 'color',[0.85 0.33 0.1]) plot(xValues(measureErrorIdx),MeasuredAltitude(measureErrorIdx),'x','MarkerSize',10, ... 'linewidth',2,'Color','red','DisplayName','Out of Accuracy Lim.'); ylabel('Altitude Above Ground (m)'); xlabel(xAxisName) title(chartTitle) legend('Location','best'); grid on; box on; xlim([min(xValues) max(xValues)]) ylim([0 max(TruthAltitude)*1.25]) function plotAllAccuracyLimits(xValues,accuracyLimits) % From 0 to 30 m, 1 m of accuracy curve1 = zeros(size(xValues)); curve2 = zeros(size(xValues)) + 30; helperPlotAccuracyLimitsSurface(xValues, curve1, curve2, [1 0.84 0],'1m Accuracy') % From 30 to 152 m, 3% accuracy curve1 = zeros(size(xValues)) + 30; curve2 = zeros(size(xValues)) + 152; helperPlotAccuracyLimitsSurface(xValues, curve1, curve2, [1 0.88 .47],'3% Accuracy') % From 152 to 762 m, 5% accuracy curve1 = zeros(size(xValues)) + 152; curve2 = zeros(size(xValues)) + 762; helperPlotAccuracyLimitsSurface(xValues, curve1, curve2, [1 0.92 0.8],'5% Accuracy') % Accuracy limit for the flight curve1 = accuracyLimits.lowerLimit(1:end); curve2 = accuracyLimits.upperLimit(1:end); helperPlotAccuracyLimitsSurface(xValues,curve1, curve2,[0.84 0.91 0.85],'Accuracy Limit') end end function helperPlotAltimeterError(TruthAltitude,MeasuredAltitude,xValues,xAxisName,chartTitle,reverseAxis) % Plot the altimeter error % Get the accuracy limits for each truth altitude accuracyLimits = helperGetAccuracyLimits(TruthAltitude); normalizedLimits.lowerLimit = accuracyLimits.lowerLimit - TruthAltitude; normalizedLimits.upperLimit = accuracyLimits.upperLimit - TruthAltitude; % Get the measurement error measurementError = MeasuredAltitude - TruthAltitude; figure; hold on; % Accuracy limit for the flight curve1 = normalizedLimits.lowerLimit; curve2 = normalizedLimits.upperLimit; lLower = plot(xValues,curve1,'--','DisplayName','Accuracy Limits','Color','k'); plot(xValues,curve2,'--','DisplayName','Accuracy Limits','Color','k'); lPlot = plot(xValues,measurementError,'-o','DisplayName','Measurement Error', 'color',[0 0.45 0.74]); ylabel('Measurement Error (m)'); xlabel(xAxisName) title(chartTitle) legend([lLower,lPlot],'Location','best'); if reverseAxis set(gca, 'XDir','reverse') end end function accuracy = helperGetAccuracyLimits(altitude) %% Get the accuracy requirements for each altitude in the altitude % vector based on the ITU document accuracy.lowerLimit = zeros(size(altitude)); accuracy.upperLimit = zeros(size(altitude)); % 0 to 100 ft. <= +- 1m idx = altitude < 30; accuracy.lowerLimit(idx) = altitude(idx) - 1; accuracy.upperLimit(idx) = altitude(idx) + 1; % 100 to 500 ft. <= 3% idx = and((altitude >= 30) , (altitude < 152)); accuracy.lowerLimit(idx) = altitude(idx) - 0.03*altitude(idx); accuracy.upperLimit(idx) = altitude(idx) + 0.03*altitude(idx); % 500 to 2500 ft. <= 5% idx = and((altitude >= 152), (altitude < 762)); accuracy.lowerLimit(idx>0) = altitude(idx) - 0.05*altitude(idx); accuracy.upperLimit(idx>0) = altitude(idx) + 0.05*altitude(idx); % above 2500 --> No limit (let's assume 10%) idx = ((altitude >= 762)); accuracy.lowerLimit(idx>0) = altitude(idx) - 0.1*altitude(idx); accuracy.upperLimit(idx>0) = altitude(idx) + 0.1*altitude(idx); end function error = helperGetAllowedError(altitude) limits = helperGetAccuracyLimits(altitude); error = limits.upperLimit - mean([limits.upperLimit,limits.lowerLimit]); end function helperPlotAccuracyLimitsSurface(x,curve1, curve2, color,name) % Plot an accuracy limit surface x2 = [x, fliplr(x)]; inBetween = [curve1, fliplr(curve2)]; fill(x2, inBetween,color,'DisplayName',name); end function clutterResolution = helperDetermineClutterResolution(altitude) % Set the resolution so that the density of scatters is high enough to % obtain sufficient results. This is set such that ~4 scatters fall within % the required system resolution. requiredResolution = helperGetAllowedError(altitude); longerRange = altitude + requiredResolution; clutterCircleRadiusSq = (longerRange)^2-altitude^2; clutterCircleArea = clutterCircleRadiusSq*pi; numScatterers = 4; clutterResolution = sqrt(clutterCircleArea/numScatterers); end
You can also select a web site from the following list
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)