# 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 waveform`SamplingRate`

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 `receiverSamplingRa`

te 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 position`

The 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 = 433.6 m, True Altitude = 429.0 m, Altitude Error = -4.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