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 values

averaged_r = [];

averaged_v = [];

% Initialize variables to keep track of current group

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

averaged_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