radar_target_generation_and_detection.m

clear all
clc;

%% Radar Specifications 
%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Frequency of operation = 77GHz
% Max Range = 200m
% Range Resolution = 1 m
% Max Velocity = 70 m/s
%%%%%%%%%%%%%%%%%%%%%%%%%%%

%speed of light = 3e8
%% User Defined Range and Velocity of target

% define the target's initial position and velocity. Note : Velocity
% remains contant
INITIAL_RANGE = 110; %meters initial range
INITIAL_VELOCITY = -20; %m/s initial velocity/speed of the target


%% FMCW Waveform Generation

%Design the FMCW waveform by giving the specs of each of its parameters.
% Calculate the Bandwidth (B), Chirp Time (Tchirp) and Slope (slope) of the FMCW
% chirp using the requirements above.
%Operating carrier frequency of Radar 
Fc= 77e9;             %carrier freq
MAX_RANGE = 200;
MAX_VELOCITY = 70; %m/s
RANGE_RESOLUTION = 1; %1meter
C = 3e8;

B_SWEEP = C/2 * RANGE_RESOLUTION;
T_CHIRP = 5.5 * 2 * MAX_RANGE / C;
SLOPE = B_SWEEP / T_CHIRP; 

%The number of chirps in one sequence. Its ideal to have 2^ value for the ease of running the FFT
%for Doppler Estimation. 
Nd=128;                   % #of doppler cells OR #of sent periods % number of chirps

%The number of samples on each chirp. 
Nr=1024;                  %for length of time OR # of range cells

% Timestamp for running the displacement scenario for every sample on each
% chirp
time_vector = linspace(0, Nd * T_CHIRP, Nr * Nd); %total time for samples

%Creating the vectors for Tx, Rx and Mix based on the total samples input.
Tx = zeros(1, length(time_vector)); %transmitted signal
Rx = zeros(1, length(time_vector)); %received signal
Mix = zeros(1, length(time_vector)); %beat signal

%Similar vectors for range_covered and time delay.
range_covered = zeros(1, length(time_vector));
time_delay = zeros(1, length(time_vector));


%% Signal generation and Moving Target simulation
% Running the radar scenario over the time. 

for i=1:length(time_vector)         
    
    %For each time stamp update the Range of the Target for constant velocity. 
    range_covered(i) = INITIAL_RANGE + INITIAL_VELOCITY * time_vector(i);
    time_delay(i) = 2 * range_covered(i) / C;
 
    %For each time sample we need update the transmitted and
    %received signal.
    % FM of CW takes place in Tx and Rx
    Tx(i) = cos(2 * pi * (Fc * time_vector(i) + SLOPE * time_vector(i)^2/2));
    Rx(i) = cos(2 * pi * (Fc * (time_vector(i)-time_delay(i)) + SLOPE * (time_vector(i)-time_delay(i))^2/2));
    
    %Now by mixing the Transmit and Receive generate the beat signal
    %This is done by element wise matrix multiplication of Transmit and
    %Receiver Signal
    Mix(i) = Tx(i) * Rx(i);
    
end

%% RANGE MEASUREMENT

%reshape the vector into Nr*Nd array. Nr and Nd here would also define the size of
%Range and Doppler FFT respectively.

Mix_reshape = reshape(Mix, Nr, Nd);

%run the FFT on the beat signal along the range bins dimension (Nr) and
fft_range = fft(Mix_reshape, Nr);

% Take the absolute value of FFT output
fft_range_abs = abs(fft_range);

% Output of FFT is double sided signal, but we are interested in only one side of the spectrum.
% Hence we throw out half of the samples.
fft_range_normalized = fft_range_abs ./ max(fft_range_abs) ; %elemewise normalization 
%fft_range_half_spectrum = fft_range_normalized(1:Nr/2, 1:2);
fft_range_half_spectrum = fft_range_normalized(1:Nr/2); %beat signal in a single half col. 

%plotting the range
figure ('Name', 'Range from First FFT(Range that is extracted from single col.)')
subplot(2,1,1)

 % plot FFT output 
 plot(fft_range_half_spectrum); 
axis ([0 200 0 1]);


%% RANGE DOPPLER RESPONSE
% The 2D FFT implementation is already provided here. This will run a 2DFFT
% on the mixed signal (beat signal) output and generate a range doppler
% map.You will implement CFAR on the generated RDM


% Range Doppler Map Generation.

% The output of the 2D FFT is an image that has reponse in the range and
% doppler FFT bins. So, it is important to convert the axis from bin sizes
% to range and doppler based on their Max values.

Mix = reshape(Mix, [Nr,Nd]);

% 2D FFT using the FFT size for both dimensions.
sig_fft2 = fft2(Mix, Nr, Nd);

% Taking just one side of signal from Range dimension.
sig_fft2 = sig_fft2(1:Nr/2, 1:Nd);
sig_fft2 = fftshift (sig_fft2);
RDM = abs(sig_fft2);
RDM = 10 * log10(RDM) ;

%use the surf function to plot the output of 2DFFT and to show axis in both
%dimensions
doppler_axis = linspace(-100, 100, Nd);
range_axis = linspace(-200,200,Nr/2) * ((Nr/2)/400);
figure('Name', 'Range-Doppler Map');
surf(doppler_axis, range_axis, RDM);
xlabel('Velocity');
ylabel('Range');
zlabel('RDM Amplitude')

%% CFAR implementation

%Slide Window through the complete Range Doppler Map

%Select the number of Training Cells in both the dimensions.
Tr = 12; %Training Cells in range dimension meaning along changing cols dir.
Td = 10; %Training Cells in doppler dimension meaning along changing rows dir.

%Select the number of Guard Cells in both dimensions around the Cell under 
%test (CUT) for accurate estimation
Gr = 4;
Gd = 4;

% offset the threshold by SNR value in dB
Offset = 1.2; 

%Create a vector to store noise_level for each iteration on training cells
noise_level = zeros(1,1);

%design a loop such that it slides the CUT across range doppler map by
%giving margins at the edges for Training and Guard Cells.
%For every iteration sum the signal level within all the training
%cells. To sum convert the value from logarithmic to linear using db2pow
%function. Average the summed values for all of the training
%cells used. After averaging convert it back to logarithimic using pow2db.
%Further add the offset to it to determine the threshold. Next, compare the
%signal under CUT with this threshold. If the CUT level > threshold assign
%it a value of 1, else equate it to 0.

% Use RDM[x,y] as the matrix from the output of 2D FFT for implementing
  % CFAR
RDM = RDM / max(RDM(:)); %normalizing RDM

for i = Tr+Gr+1:(Nr/2)-(Gr+Tr)  %Nr/2 because of taking one side of the spectrum see line 127 and 128
    for j = Td+Gd+1:Nd-(Gd+Td)
        noise_level = zeros(1,1);
        
        for p = i-(Tr+Gr) : i+(Tr+Gr)
            for q = j-(Td+Gd) : j+(Td+Gd)
                if (abs(i-p) > Gr || abs(j-q) > Gd)
                    noise_level = noise_level + db2pow(RDM(p,q));
                end
            end
        end
        
        threshold = pow2db(noise_level / (2 * (Tr+Gr+1) * 2 * (Td+Gd+1) - (Gr * Gd)-1)); %pow2db on normalized noise level 
        threshold = threshold + Offset;
        CUT = RDM(i, j);
        
        if (CUT < threshold)
            RDM(i, j) = 0;
        else
            RDM(i, j) = 1;
        end
        
    end
end

% The process above will generate a thresholded block, which is smaller 
%than the Range Doppler Map as the CUT cannot be located at the edges of
%matrix. Hence,few cells will not be thresholded. To keep the map size same
% set those values to 0. 
 
RDM(RDM ~= 0 & RDM ~= 1) = 0;

%display the CFAR output using the Surf function like we did for Range
%Doppler Response output.
figure('Name', '2D-CFAR');
surf(doppler_axis, range_axis, RDM);
colorbar;
xlabel('Velocity');
ylabel('Range');
zlabel('Normalized RDM');