How to modify this code from MATHWORKS website?

I am trying to simulate the following code from mathworks.com:

https://www.mathworks.com/help/radar/ug/fmcw-radar-altimeter-simulation.html

But not with fmcw waveform. I want to simulate this with the following waveform as altimeter waveform.

h = phased.PhaseCodedWaveform;
h.SampleRate = 100000000;
h.Code = 'Barker';
h.ChipWidth = 2e-08;
h.NumChips = 13;
h.PRF = 50000;
h.NumPulses = 1;

But i don't understand how to modify the whole code and other parameters according to this waveform

Answer 1

[Satwik]

Hi Basmah!

To simulate the radar altimeter example with a phase-coded waveform instead of an FMCW waveform, key modifications are required. These involve replacing the waveform generation, adjusting the signal processing steps, and ensuring the matched filter is properly configured for the new waveform. Here is an example of how you can do this:

% Define the phase-coded waveform
h = phased.PhaseCodedWaveform('SampleRate', 100000000, 'Code', 'Barker', ...
'ChipWidth', 2e-08, 'NumChips', 13, 'PRF', 50000, 'NumPulses', 1);

% Generate the waveform
waveform = h();

% Create a matched filter
mf = phased.MatchedFilter('Coefficients', getMatchedFilter(h), ...
'GainOutputPort', true);

% Radar parameters
fc = 10e9; % Carrier frequency
transmitter = phased.Transmitter('PeakPower', 1e3, 'Gain', 20);
receiver = phased.ReceiverPreamp('Gain', 20, 'NoiseFigure', 10);

% Simulation setup
range = 1000; % Target range in meters
target = phased.RadarTarget('MeanRCS', 1, 'OperatingFrequency', fc);
targetPlatform = phased.Platform('InitialPosition', [0; 0; range]);
radarPlatform = phased.Platform('InitialPosition', [0; 0; 0]);
channel = phased.FreeSpace('SampleRate', h.SampleRate, ...    
'OperatingFrequency', fc);



% Simulation loop
numPulses = 10;
rxSignal = zeros(size(waveform, 1), numPulses);

for n = 1:numPulses    
    [radarPos, radarVel] = radarPlatform(1/h.PRF);
    [tgtPos, tgtVel] = targetPlatform(1/h.PRF);
    
    txsig = transmitter(waveform);
    txsig = channel(txsig, radarPos, tgtPos, radarVel, tgtVel);
    rxsig = target(txsig);
    rxsig = channel(rxsig, tgtPos, radarPos, tgtVel, radarVel);
    rxSignal(:, n) = receiver(rxsig);
end

% Apply matched filter
mfOut = mf(rxSignal);

% Display results
plot(real(mfOut));
xlabel('Samples');
ylabel('Amplitude');
title('Matched Filter Output');

Hope this helps!

Answer 2

[George]

Hello Basmah,

There are a couple of things that you will need to change in order to get this example working with the phase coded waveform instead of the FMCW.

First of all, the radarTransceiver is used to generate IQ data from a radarScenario. You will need to replace the waveform in the radarTransceiver with your desired waveform.

So this line:

altimeterWaveform = phased.FMCWWaveform('SweepBandwidth',bandwidth,'SampleRate',waveformSamplingRate,'SweepDirection',sweepDir,'SweepTime',tSweep,'OutputFormat','Sweeps');

Should be replaced with this line:

altimeterWaveform = phased.PhaseCodedWaveform(SampleRate=100e6,Code='Barker',ChipWidth=2e-8,NumChips=13,PRF=50000,NumPulses=1);

Also, the radarTransceiver expects that the sample rate setting on the waveform matches the setting on the receiver, so you can update this line:

altimeterReceiver = phased.ReceiverPreamp('SampleRate', waveformSamlpingRate,'NoiseFigure', noiseFigure);

To this:

altimeterReceiver = phased.ReceiverPreamp('SampleRate', 100e6,'NoiseFigure', noiseFigure);

Finally, you will need to implement your own signal processing routine. This example uses a helper class (helperAltimeterSignalProcessor) to generate detections. However, this class has the key underlying assumption that an FMCW is being transmitted.

So you have two options. You could implement your own class that has the same methods as the helperAltimeterSignalProcessor, or you could remove references to this object and write your own signal processing pipeline.

For example, here is a script which cuts out all of the signal processing steps in the example, and just generates the IQ data. You can try running this from the folder where the example is located. In this case you will need to write your own signal processing routine, you can look at the helperAltimeterSignalProcessor for inspiration, but this script will generate IQ data using your waveform:

%% Setup parameters

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)

%% Create antenna

% 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);

altimeterTxAntenna = phased.Radiator('OperatingFrequency',fCenter,'Sensor',altimeterAntennaElement, 'SensorGainMeasure','dBi');
altimeterRxAntenna = phased.Collector('OperatingFrequency',fCenter,'Sensor',altimeterAntennaElement,'SensorGainMeasure','dBi');

%% Create waveform

altimeterWaveform = phased.PhaseCodedWaveform(SampleRate=100e6,Code='Barker',ChipWidth=2e-8,NumChips=13,PRF=50000,NumPulses=1);

%% Create Transmitter and Receiver

altimeterTransmitter = phased.Transmitter('PeakPower',pTransmitter);
altimeterReceiver = phased.ReceiverPreamp('SampleRate', 100e6,'NoiseFigure', noiseFigure);

%% Create radarTransceiver

altimeterTransceiver = radarTransceiver('Waveform',altimeterWaveform,...
    'Transmitter',altimeterTransmitter, ...
    'TransmitAntenna',altimeterTxAntenna, ...
    'ReceiveAntenna',altimeterRxAntenna, ...
    'Receiver',altimeterReceiver, ...
    'MountingAngles', [0 -90 0], ...
    'NumRepetitions', 1);

%% Create scenario

% Create scenario
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);

% Creatge height map
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);

% Create platform
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]));

%% Run Simulation

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

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

    % This is where you would do your signal processing


    cnt = cnt + 1;
end

%% Helpers

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

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