PRE2024 3 Group18/MatLab code
Here is the MatLab code of the simulation:
Mfsk.mlx
close all;
clear;
clc;
rng(2025);
[fmcwwaveform,target,tgtmotion,channel,transmitter,receiver, ...
sensormotion,c,fc,lambda,fs,maxbeatfreq] = helperMFSKSystemSetup;
mfskwaveform = phased.MFSKWaveform( ...'SampleRate',45e6, ...
'SweepBandwidth',44e6, ...
'StepTime', 2e-6, ...
'StepsPerSweep',700,...
'FrequencyOffset',-294e3/4, ...
'OutputFormat','Sweeps', ...
'NumSweeps',1);
numsamp_step = round(mfskwaveform.SampleRate*mfskwaveform.StepTime);
sig_display = mfskwaveform();
spectrogram(sig_display(1:8192),kaiser(3*numsamp_step,100), ...
ceil(2*numsamp_step),linspace(0,4e6,2048),mfskwaveform.SampleRate, ...
'yaxis','reassigned','minthreshold',-60)
Nsweep = 1;
release(channel);
channel.SampleRate = mfskwaveform.SampleRate;
release(receiver);
receiver.SampleRate = mfskwaveform.SampleRate;
xr = helperFMCWSimulate(Nsweep,mfskwaveform,sensormotion,tgtmotion, ...
transmitter,channel,target,receiver);
x_dechirp = reshape(xr(numsamp_step:numsamp_step:end),2,[]).';
fs_dechirp = 1/(2*mfskwaveform.StepTime);
xf_dechirp = fft(x_dechirp);
num_xf_samp = size(xf_dechirp,1);
beatfreq_vec = (0:num_xf_samp-1).'/num_xf_samp*fs_dechirp;
clf
subplot(211),plot(beatfreq_vec/1e3,abs(xf_dechirp(:,1)))
grid on
ylabel('Magnitude')
title('Frequency spectrum for sweep 1')
subplot(212),plot(beatfreq_vec/1e3,abs(xf_dechirp(:,2)))
grid on
ylabel('Magnitude')
title('Frequency spectrum for sweep 2')
xlabel('Frequency (kHz)')
cfar = phased.CFARDetector('ProbabilityFalseAlarm',10e-2, ...
'NumTrainingCells',8);
peakidx = cfar(abs(xf_dechirp(:,1)),1:num_xf_samp);
Fbeat = beatfreq_vec(peakidx);
phi = angle(xf_dechirp(peakidx,2))-angle(xf_dechirp(peakidx,1));
sweep_slope = mfskwaveform.SweepBandwidth/ ...
(mfskwaveform.StepsPerSweep*mfskwaveform.StepTime);
temp = ...
[-2/lambda 2*sweep_slope/c;-2/lambda*mfskwaveform.StepTime 2*mfskwaveform.FrequencyOffset/c]\ ...
[Fbeat phi/(2*pi)].';
r_est = temp(2,:)
v_est = temp(1,:)*3.6
function [filtered_r, filtered_v] = filter_values(r_est, v_est)
% Filter out values where speed is larger than 100 km/h or distance is larger than 200 meters
% Filters out the vehicles that passed the person already.
valid_indices = (abs(v_est) <= 100) & (r_est <= 200) & (v_est >= 0);
filtered_r = r_est(valid_indices);
filtered_v = v_est(valid_indices);
end
function [averaged_r, averaged_v] = average_similar_values(r_est, v_est)
% Combine r_est and v_est into a single matrix
combined = [r_est(:), v_est(:)];
% Sort the matrix based on the first column (r_est values)combined = sortrows(combined, 1);
% Initialize arrays to store averaged valuesaveraged_r = [];
averaged_v = [];
% Initialize variables to keep track of current groupcurrent_r = combined(1, 1);
current_v = combined(1, 2);
count = 1;
for i = 2:length(combined)if abs(combined(i, 1) - current_r) <= 4.5 & abs(combined(i, 2) - current_v) <= 4.5 % If the value is similar (within 4 meters and 4 km/h)
current_r = [current_r, combined(i, 1)];
current_v = [current_v, combined(i, 2)];
count = count + 1;
else
averaged_r = [averaged_r, mean(current_r)];
averaged_v = [averaged_v, mean(current_v)];
current_r = combined(i, 1);
current_v = combined(i, 2);
count = 1;
end
end
% Append the last groupaveraged_r = [averaged_r, mean(current_r)];
averaged_v = [averaged_v, mean(current_v)];
end
function arrival_times = calculate_arrival_times(r_est, v_est)
% Calculate the time it takes for each vehicle to arrive at the position
arrival_times = r_est ./ (v_est./3.6);
end
% Filter values
[filtered_r, filtered_v] = filter_values(r_est, v_est);
% Calculate averaged values
[averaged_r, averaged_v] = average_similar_values(filtered_r, filtered_v);
% Calculate arrival times
arrival_times = calculate_arrival_times(averaged_r, averaged_v);
disp('Averaged r values:');
disp(averaged_r);
disp('Averaged v values:');
disp(averaged_v);
disp('Arrival times (in seconds):');
disp(arrival_times);
function result = check_time_below_threshold(arrival_times, threshold)
% Check if there is any time below the threshold
if any(arrival_times < threshold)
result = 1;
disp('Stop');
else
result = 0;
disp('Go')
end
end
result = check_time_below_threshold(arrival_times, 8);
helperMFSKSystemSetup.m
function [fmcwwaveform,target,tgtmotion,channel,transmitter,receiver,...
sensormotion,c,fc,lambda,fs,fr_max] = helperMFSKSystemSetup
% This function helperMFSKSystemSetup is only in support of MFSKExample.
% System parameter
fc = 78e9; % operating frequency
c = 3e8; % propagation speed
lambda = c/fc; % wavelength
tm = 2222e-11; % sweep time
%range_res = 1; % range resolution
bw = 5e9;%rangeres2bw(range_res,c); % bandwidth
sweep_slope = bw/tm; % sweep slope
range_max = 250;
fr_max = range2beat(range_max,sweep_slope,c);
v_max = 210*1000/3600;
fd_max = speed2dop(2*v_max,lambda);
fb_max = fr_max+fd_max;
fs = max(2*fb_max,bw);
fs = 7.5e11;
fmcwwaveform = phased.FMCWWaveform('SweepTime',tm,'SweepBandwidth',bw,...
'SampleRate',fs,'SweepDirection','Triangle');
car_dist = 60;
car_speed = 55*1000/3600;
car_rcs = db2pow(min(10*log10(car_dist)+5,20));
truck_dist = 40;
truck_speed = -50*1000/3600;
truck_rcs = db2pow(min(10*log10(truck_dist)+5,20));
bike_dist = 35;
bike_speed = -20*1000/3600;
bike_rcs = db2pow(min(10*log10(bike_dist)+5,20));
tgtpos = [[car_dist;0;0],[truck_dist;0;0],[bike_dist;0;0]];
tgtvel = [[car_speed;0;0],[truck_speed;0;0],[bike_speed;0;0]];
tgtmotion = phased.Platform('InitialPosition',tgtpos,'Velocity',tgtvel);
tgtrcs = [car_rcs,truck_rcs,bike_rcs];
target = phased.RadarTarget('MeanRCS',tgtrcs,'PropagationSpeed',c,...
'OperatingFrequency',fc);
channel = phased.FreeSpace('PropagationSpeed',c,...
'OperatingFrequency',fc,'SampleRate',fs,'TwoWayPropagation',true);
ant_aperture = 18.5*lambda*0.8*lambda+18.5/38*8*lambda*0.4*lambda; %6.06e-4;
ant_gain = aperture2gain(ant_aperture,lambda) ; % in dB
tx_ppower = db2pow(13)*1e-3; % in watts
tx_gain = 9+ant_gain; % in dB
rx_gain = 44; % in dB
rx_nf = 11.5; % in dB
transmitter = phased.Transmitter('PeakPower',tx_ppower,'Gain',tx_gain);
receiver = phased.ReceiverPreamp('Gain',rx_gain,'NoiseFigure',rx_nf,...
'SampleRate',fs);
radar_speed = 0*1000/3600;
sensormotion = phased.Platform('Velocity',[radar_speed;0;0]);
HelperFMCWSimulate.m
function [xr,xr_unmixed] = helperFMCWSimulate(Nsweep,waveform,...
radarmotion,carmotion,transmitter,channel,cartarget,receiver)
% This function helperFMCWSimulate is only in support of FMCWExample. It
% may be removed in a future release.
% RSWEEP =
% helperFMCWSimulate(NSWEEP,WAVEFORM,RADARMOTION,CARMOTION,TRANSMITTER,
% CHANNEL,CARTARGET,RECEIVER) returns the simulated sweep train RSWEEP.
%
% The input parameters are:
% NSWEEP: number of sweeps
% WAVEFORM: waveform object
% RADARMOTION: platform object for the radar
% CARMOTION: platform object for target car
% TRANSMITTER: transmitter object
% CHANNEL: propagation channel object
% CARTARGET: target car object
% RECEIVER: receiver object
%
% The rows of RSWEEP represent fast time and its columns represent slow
% time (pulses). When the pulse transmitter uses staggered PRFs, the
% length of the fast time sequences is determined by the highest PRF.
% Copyright 2010-2016 The MathWorks, Inc.
release(waveform);
release(radarmotion);
release(transmitter);
release(receiver);
release(carmotion);
release(channel);
release(cartarget);
if isa(waveform,'phased.MFSKWaveform')
sweeptime = waveform.StepTime*waveform.StepsPerSweep;
else
sweeptime = waveform.SweepTime;
end
Nsamp = round(waveform.SampleRate*sweeptime);
xr = complex(zeros(Nsamp,Nsweep));
xr_unmixed = xr;
Ntgt = numel(cartarget.MeanRCS);
for m = 1:Nsweep
% Update radar and target positions
[radar_pos,radar_vel] = radarmotion(sweeptime);
[tgt_pos,tgt_vel] = carmotion(sweeptime);
% Transmit FMCW waveform
sig = waveform();
txsig = transmitter(sig);
% Propagate the signal and reflect off the target
rxsig = complex(zeros(Nsamp,Ntgt));
for n = 1:Ntgt
rxsig(:,n) = channel(txsig,radar_pos,tgt_pos(:,n),radar_vel,tgt_vel(:,n));
end
rxsig = cartarget(rxsig);
% Dechirp the received radar return
rxsig = receiver(sum(rxsig,2));
xd = dechirp(rxsig,sig);
xr_unmixed(:,m) = rxsig;
xr(:,m) = xd;
end