0
$\begingroup$

I'm trying to learn radar simulation using Matlab. I used this example in Matlab https://www.mathworks.com/help/driving/ug/radar-signal-simulation-and-processing-for-automated-driving.html

But when I tried to adapt to a new carrier frequency at 4.7MHz and bandwidth of 100kHz. The numeric range calculation shows the correct target position and velocity. Here is my code

close all;
clear;

%% System parameters

fc = 4.7e6;                     % Carrier frequency (Hz)
c = physconst('LightSpeed');    % Propagation speed (m/s)
lambda = c/fc;                  % Wavelength (m)

nm2m = 1852;                   % Coefficient to exchange from nautical miles to m
rangeMaxNm = 200;               % Maximum range in Nm
rangeMax = rangeMaxNm*nm2m;    % Max range in m

% Set up chirp duration.
chirp_dur = 5.5;                  % Setup value
tm_t = range2time(rangeMax,c);
tm = chirp_dur*range2time(rangeMax,c);    % Chirp duration. 

% Calculate range solution from bandwidth. Could be other way around
bw = 100e3;                     % Bandwidth (Hz)
range_res = bw2range(bw,c);     % Range resolution (m)

% Set the sampling rate to satisfy both the range and velocity requirements
% for the radar
sweepSlope = bw/tm;             % FMCW sweep slope (m)
fbeatMax = range2beat(rangeMax, sweepSlope, c); %Maximum beat frequency

vMaxKn = 136;                   % Maximum velocity (knots)
vMax = vMaxKn * nm2m /3600;    % Maximum velocity (m/s)
fdopMax = speed2dop(2*vMax,lambda); %Maximum Doppler shift (Hz)

fifMax = fbeatMax + fdopMax;    %Maximum received IF (Hz)
fs= max(2*fifMax,be);
test = round(fs * tm);
tm = test/fs;
%Configure FMCW waveform
waveform = phased.FMCWWaveform(...
    'SweepTime',tm,...
    'SweepBandwidth',be,...
    'SampleRate', fs,...
    'SweepDirection','Up');
if strcmp(waveform.SweepDirection, 'Down')
    sweepSlope = -sweepSlope;
end

Nsweep = 256;
sig = waveform();
figure(1);
subplot(211); plot(0:1/fs:tm-1/fs,real(sig));
xlabel('Time (s)'); ylabel('Amplitude (v)');
title('FMCW signal'); axis tight;
subplot(212); spectrogram(sig,32,16,32,fs,'yaxis');
title('FMCW signal spectrogram');
%%

%% Antenna model
%Model antenna element
antElement = phased.IsotropicAntennaElement('BackBaffled',true);

%Construct the receive array
Nrow = 2;                       % Number of rows
Ncolumn = 16;                   % Number of elements in a row
drow = lambda/2;                % row spacing
dcolumn = lambda/2;             % element spacing

rxArray = phased.URA('Size', [Nrow Ncolumn],...
    'Lattice', 'Rectangular',...
    'ElementSpacing', [drow dcolumn]);

% Form forward-facing beam to detect objects in front of the array
beamformer = phased.PhaseShiftBeamformer('SensorArray',rxArray,...
   'PropagationSpeed',c,'OperatingFrequency',fc,'Direction',[90;0]);

% Half-power beamwidth of the receive array - Needed to verify
hpbw = helperAutoDrivingRadarSigProc('Array Beamwidth',rxArray,c,fc);

%% Estimate direction of arrival
% Direction-of-arrival estimator for phased array signals
%estimator = phased.MUSICEstimator2D('SensorArray',rxArray,...
%    'OperatingFrequency',fc,...
%    'NumSignalsSource','Property',...
%    'DOAOutputPort',true,'NumSignals',2,...
%    'AzimuthScanAngles',-50:.5:50,...
%    'ElevationScanAngles',-30:.5:30);
%[~,doas] = estimator(x + noise)
%Scan beam in front
%angscan = -80:80;
%beamscan = phased.PhaseShiftBeamformer('Direction', [angscan;0*angscan],...
%    'SensorArray',rxArray,...
%    'OperatingFrequency',fc);
%%
%% Transmitter/receiver model
txPkPower = 0.4e3;       % Tx peak power (W)
txGain = 6;              % Tx antenna gain (dB)

antAperture = 6.06e-4;
antGain = 6;             % Antenna gain (dB)
rxNF = 2.9;              % Receiver noise figure (dB)
rxGain = antGain;
%Waveform transmitter
transmitter = phased.Transmitter('PeakPower',txPkPower, 'Gain',txGain,'InUseOutputPort',true);

%Radiator for single transmit element
radiator = phased.Radiator('Sensor', antElement, 'OperatingFrequency', fc);

%Collector for receive array
collector = phased.Collector('Sensor', rxArray, 'OperatingFrequency', fc);

%Receiver preamplifier
receiver = phased.ReceiverPreamp('Gain', rxGain, 'NoiseFigure', rxNF, 'SampleRate', fs);

radarStation = phased.Platform(...
    'InitialPosition',[0;0;0],...
    'Velocity',[0;0;0]);

%% Target model
target1_disNm = 100;                % Target 1 distance (Nautical miles)
target1_dis = target1_disNm * nm2m;   %Target distcane (m)
target2_disNm = 200;                % Target 1 distance (Nautical miles)
target2_dis = target2_disNm * nm2m;   %Target distcane (m)
target1_speed = vMax;
target2_speed = vMax/2;
%target_rcs = db2pow(min(10*log10(target_dis)+5,20));
target_rcs = 1000;
target = phased.RadarTarget(...
    'Model', 'Nonfluctuating',...
    'MeanRCS', target_rcs,...
    'PropagationSpeed',c,...
    'OperatingFrequency',fc);
%tgtpos = [[target1_dis;0;0],[target2_dis;0;0],[50000;0;0]];
tgtpos = [target1_dis;0;0];
%tgtvel = [[target1_speed;0;0],[target2_speed;0;0],[0;0;0]];
tgtvel = [target1_speed;0;0];
targetMotion = phased.Platform(...
    'InitialPosition', [target1_dis;0;0],...
    'Velocity',[30;0;0]);
[tgtrng,tgtang] = rangeangle(targetMotion.InitialPosition,...
    radarStation.InitialPosition);

%% Propagation Environment
channel = phased.FreeSpace(...
    'PropagationSpeed',c,...
    'SampleRate',fs,...
    'TwoWayPropagation',true,...
    'OperatingFrequency',fc);
%% Signal processing chain
Nft = waveform.SweepTime * waveform.SampleRate; % Number of fast-time samples
Nst = 56;                                   % Number of slow time samples
Nr = 2^nextpow2(Nft);                           % Number of range samples after processing
Nd = 512;                                       % Number of Doppler samples after processing

xr = complex(zeros(waveform.SampleRate*waveform.SweepTime,Nsweep));
%specanalyzer = dsp.SpectrumAnalyzer('SampleRate',fs,...
%   'PlotAsTwoSidedSpectrum',true,...
%   'Title','Spectrum for received and dechirped signal',...
%   'ShowLegend',true);

for m = 1:Nsweep
   %Update target position
    [radar_pos, radar_vel] = radarStation(waveform.SweepTime);
    [tgt_pos, tgt_vel] = targetMotion(waveform.SweepTime);

    %Transmit FMCW waveform
    sig = waveform();
    [txsig,txstatus] = transmitter(sig);

    % Propagate the signal and reflect off the target
    txsig = radiator(txsig,tgtang);
    txsig = channel(txsig, radar_pos, tgt_pos,radar_vel,tgt_vel);
    txsig = target(txsig);

    %Dechirp the received radar return
    %txsig = collector(txsig,tgtang);
    txsig = receiver(txsig);
    dechirpsig = dechirp(txsig,sig);

    %Visualize the spectrum
    %specanalyzer([txsig dechirpsig]);

    xr(:,m) = dechirpsig;
end

rngDopResp = phased.RangeDopplerResponse(...
    'RangeMethod', 'FFT',...
    'DopplerOutput', 'Speed',...
    'RangeWindow', 'Hann',...
    'RangeFFTLengthSource', 'Property',...
    'RangeFFTLength', Nr,...
    'DopplerFFTLengthSource', 'Property',...
    'DopplerFFTLength', Nd,...
    'DopplerWindow','Hann',...
    'PropagationSpeed',c,...
    'OperatingFrequency',fc,...
    'SampleRate',fs);
clf;
figure(2)
plotResponse(rngDopResp, xr);
%axis([-vMax vMax 0 rangeMax])
clim = caxis;

Dn = fix(fs/(2*fbeatMax));
for m = size(xr,2):-1:1
     xr_d(:,m) = decimate(xr(:,m),Dn,'FIR');
end
fs_d = fs/Dn;
fb_rng = rootmusic(pulsint(xr_d,'coherent'),1,fs_d);
rng_est = beat2range(fb_rng,sweepSlope,c)
rng_estNm = rng_est/nm2m;

peak_loc = val2ind(rng_est,c/(fs_d*2));
fd = -rootmusic(xr_d(peak_loc,:),1,1/tm);
v_est = dop2speed(fd,lambda)/2

However, in the Range Doppler response plot, only the velocity is correct, the range display is incorrect. It should be at 184.5km and -30 m/s. I tried to change the range but the range doppler stayed the same. If I change the velocity then it changed correctly. enter image description here

Do I need to do any scaling with the data of range doppler response? Any suggestion? Thank you for your time!

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.