669 lines
24 KiB
Matlab
669 lines
24 KiB
Matlab
clear variables;
|
||
close all;
|
||
clc
|
||
|
||
addpath(genpath(fullfile(pwd, 'Functions')));
|
||
addpath(genpath(fullfile(pwd, 'Antenna_Pattern')));
|
||
|
||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||
%
|
||
% 좌표계 : ego car 기준 앞뒤 : x 축, 좌우 : y 축, 위아래 : z 축
|
||
%
|
||
%
|
||
%
|
||
%
|
||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||
|
||
%% Control Panel
|
||
opt_waveform = 1; % 0 : cos waveform, 1 : exp waveform
|
||
opt_interference = 0;
|
||
opt_clutter = 0;
|
||
|
||
opt_tgt = 0; % 0 : Point targets, 1 : Realistic targerts
|
||
opt_noise = 1;
|
||
opt_filter = 1;
|
||
opt_window = 1;
|
||
|
||
flag_plot_physical_array = 0;
|
||
flag_plot_virtual_array = 0;
|
||
flag_plot_RVM = 1;
|
||
flag_plot_noise_esti = 0;
|
||
flag_plot_birdview = 1;
|
||
|
||
%% Basic Parameters
|
||
c0 = physconst('Lightspeed'); % Light speed [m/s]
|
||
kb = physconst('Boltzmann'); % Boltzmann Constant
|
||
T0 = 290; % Room Temperature [K]
|
||
|
||
%% System Parameters
|
||
fc = 76.75e9; % Center frequency [Hz]
|
||
lambda_c = c0 / fc; % wavelength for center frequency [m]
|
||
fs = 20e6; % ADC sampling rate [Hz]
|
||
NumChirps = 256; % Number of Chirps (= # of samples in slow-time)
|
||
NumSamples = 1024; % Number of samples in fast-time
|
||
rngFFTlen = 2^(nextpow2(NumSamples)); % FFT length in fast-time
|
||
velFFTlen = 2^(nextpow2(NumChirps)); % FFT length in slow-time
|
||
|
||
maxrngFFTidx = 13/16 * rngFFTlen / 2; % Max Range index in FFT
|
||
|
||
%% Array Design
|
||
Nt = 3; % Number of Tx
|
||
Nr = 4; % Number of Rx
|
||
|
||
u_azi = 1.96e-3; % [m]
|
||
u_elv = 5.488e-3; % [m]
|
||
|
||
txarray_loc = [ 0 0 0 ;
|
||
0 -8*u_azi 0 ;
|
||
0 -15*u_azi u_elv;]';
|
||
rxarray_loc = [ 0 0 10*u_azi;
|
||
0 -5*u_azi 10*u_azi;
|
||
0 -9*u_azi 10*u_azi;
|
||
0 -15*u_azi 10*u_azi;]';
|
||
|
||
[mimoarray] = gen_virarray(txarray_loc, rxarray_loc);
|
||
|
||
if flag_plot_physical_array == 1
|
||
ant_width = 2.55e-3;
|
||
ant_height = 10.77e-3;
|
||
|
||
figure
|
||
for eidx = 1 : Nt
|
||
scatter((txarray_loc(2, eidx)+ant_width/2)*1e3, (txarray_loc(3, eidx)+ant_height/2)*1e3, 'ro');
|
||
hold on
|
||
rectangle('Position', 1e3*[(txarray_loc(2, eidx)), (txarray_loc(3, eidx)), ant_width, ant_height]);
|
||
hold on
|
||
text(txarray_loc(2,eidx)*1e3, txarray_loc(3,eidx)*1e3, ['#', num2str(eidx)], 'Color', 'red', 'HorizontalAlignment', 'center');
|
||
end
|
||
grid on
|
||
|
||
for eidx = 1 : Nr
|
||
scatter((rxarray_loc(2, eidx)+ant_width/2)*1e3, (rxarray_loc(3, eidx)+ant_height/2)*1e3, 'bo');
|
||
hold on
|
||
rectangle('Position', 1e3*[(rxarray_loc(2, eidx)), (rxarray_loc(3, eidx)), ant_width, ant_height]);
|
||
hold on
|
||
text(rxarray_loc(2,eidx)*1e3, rxarray_loc(3, eidx)*1e3, ['#', num2str(eidx)], 'Color', 'blue', 'HorizontalAlignment', 'center');
|
||
end
|
||
xlabel('y-axis [mm]');
|
||
ylabel('y-axis [mm]');
|
||
title('Physical array locations');
|
||
end
|
||
|
||
if flag_plot_virtual_array == 1
|
||
figure
|
||
scatter(mimoarray(2,:) / lambda_c, mimoarray(3,:) / lambda_c, 'o')
|
||
grid on
|
||
xlabel('y-axis [\lambda]');
|
||
ylabel('z-axis [\lambda]');
|
||
title('Virtual Array Positions');
|
||
end
|
||
|
||
%% Environments
|
||
% 00. Radar platform
|
||
posR = [0 0 0].'; % radar position
|
||
velR = [0 0 0].'; % radar velocity
|
||
accR = [0 0 0].'; % radar acceleration
|
||
|
||
% 01. Target
|
||
tgt_rng = [10];
|
||
tgt_azi = [0];
|
||
tgt_elv = [0];
|
||
tgt_pos = [tgt_rng.*cosd(tgt_azi).*cosd(tgt_elv) tgt_rng.*sind(tgt_azi).*cosd(tgt_elv) tgt_rng.*sind(tgt_elv)].';
|
||
|
||
%tgt_pos = [150 0 0].';
|
||
tgt_vel = [5 0 0].';
|
||
tgt_acc = [0 0 0].';
|
||
tgt_rcs = [20];
|
||
tgt_class = 'TCR';
|
||
tgt = gen_tgt(tgt_pos, tgt_vel, tgt_acc, tgt_rcs, 'TCR');
|
||
tgt.num = size(tgt.rcs, 1);
|
||
|
||
for tgt_idx = 1 : tgt.num
|
||
[tgt_azi_rad(tgt_idx), tgt_elv_rad(tgt_idx), ~] = cart2sph(tgt_pos(1, tgt_idx), tgt_pos(2, tgt_idx), tgt_pos(3, tgt_idx));
|
||
tgt_azi_deg(tgt_idx, 1) = rad2deg(tgt_azi_rad(tgt_idx));
|
||
tgt_elv_deg(tgt_idx, 1) = rad2deg(tgt_elv_rad(tgt_idx));
|
||
tgt_radi_vel(tgt_idx, 1) = norm((tgt_vel(:, tgt_idx) - velR) .* (tgt_pos(:, tgt_idx) - posR)/norm((tgt_pos(:, tgt_idx) - posR)));
|
||
end
|
||
|
||
tgt_table_truth = table([tgt_rng tgt_radi_vel tgt_azi_deg tgt_elv_deg]', 'VariableNames', {'Tgt.'}, 'RowNames', {'Rng[m]', 'Vel,[m/s]', 'Azi.[deg]','Elv.[deg]'});
|
||
disp(tgt_table_truth);
|
||
|
||
% 02. Clutter
|
||
|
||
% 03. Interference
|
||
|
||
%% Transmitter
|
||
% 00. Sampling interval
|
||
ts = 1/fs;
|
||
|
||
% 00. Transmitter Parameters
|
||
Ptx_dBm = 11;
|
||
|
||
% 01. Waveform generation
|
||
TxBw = 344e6;
|
||
timing_struct.T_dwell = 0e-6;
|
||
timing_struct.T_settle = 3e-6;
|
||
timing_struct.T_jumpback = 0.5e-6;
|
||
timing_struct.T_reset = 1e-6;
|
||
timing_struct.T_acq = 51.2e-6;
|
||
tx_phase_initial = 0;
|
||
timing_struct.T_idle = 9.5e-6; %[1e-6 9.5e-6 18e-6];
|
||
|
||
% timing_struct.T_dwell = 0e-6;
|
||
% timing_struct.T_settle = 0e-6;
|
||
% timing_struct.T_jumpback = 0e-6;
|
||
% timing_struct.T_reset = 0e-6;
|
||
% timing_struct.T_acq = 51.2e-6;
|
||
% tx_phase_initial = 0;
|
||
% timing_struct.T_idle = 0e-6; %[1e-6 9.5e-6 18e-6];
|
||
timing_struct.PRI = timing_struct.T_dwell + timing_struct.T_settle + timing_struct.T_jumpback + timing_struct.T_reset + timing_struct.T_acq + timing_struct.T_idle;
|
||
%[56.7e-6 65.2e-6 73.7e-6];
|
||
|
||
|
||
% [24.08.22] Lowpass filter 구현 때문에 2배 oversampling 함.
|
||
for m = 1 : NumChirps
|
||
%Tx_wave = gen_fastramp_fmcw2(timing_struct, 2*fs, fc-TxBw/2, TxBw, m, NumChirps, 0, tx_phase_initial, opt_waveform);
|
||
Tx_wave = gen_fastramp_fmcw(timing_struct, 2*fs, fc-TxBw/2, TxBw, NumChirps, 0, tx_phase_initial, opt_waveform);
|
||
Ref_wave(:,m) = Tx_wave.waveform;
|
||
end
|
||
Ref_wave = repmat(Ref_wave, 1, 1, Nr);
|
||
|
||
% 02. DDMA
|
||
DDMA_freq = [0 1/4 2/4] * 1/Tx_wave.T_chirp;
|
||
DDMA_idx = DDMA_freq * Tx_wave.T_chirp * velFFTlen;
|
||
|
||
%% System Parameters (Analysis)
|
||
rng_max = beat2rng(fs/2*13/16, Tx_wave.f_slope, c0);
|
||
rng_res = beat2rng(fs/rngFFTlen, Tx_wave.f_slope, c0);
|
||
rng_sep = bw2rngsep(TxBw, 1, c0);
|
||
|
||
vel_max = 0.5*dop2spd(1/(Tx_wave.T_chirp)/2, lambda_c);
|
||
vel_res = dop2spd(1/Tx_wave.T_frame/velFFTlen, lambda_c);
|
||
vel_sep = dop2spd(1/Tx_wave.T_frame/velFFTlen, lambda_c);
|
||
|
||
% figure(1)
|
||
% subplot(211); plot(Tx_wave.timeline,real(Tx_wave.waveform));
|
||
% xlabel('Time (s)'); ylabel('Amplitude (v)');
|
||
% title('FMCW signal'); axis tight;
|
||
% subplot(212); spectrogram(Tx_wave.waveform,32,16,32,fs,'yaxis');
|
||
% title('FMCW signal spectrogram');
|
||
|
||
%% Antenna (Tx)
|
||
load('azi_ant_pat.mat');
|
||
load('elev_ant_pat.mat');
|
||
|
||
azi_ant_pat(:,2) = azi_ant_pat(:,2) - 4;
|
||
elev_ant_pat(:,2) = elev_ant_pat(:,2) - 4;
|
||
|
||
%% Propagation (From Tx ant to Rx ant)
|
||
delayed_Tx_sig = zeros(length(Tx_wave.waveform), NumChirps, Nr);
|
||
|
||
for tgtidx = 1 : tgt.num
|
||
|
||
for m = 1 : NumChirps
|
||
|
||
if m > 1
|
||
velR = velR + accR * timing_struct.PRI;
|
||
posR = posR + velR * timing_struct.PRI;
|
||
tgt.pos(:,tgtidx) = tgt.pos(:,tgtidx) + tgt.vel(:,tgtidx) * timing_struct.PRI;
|
||
tgt.vel(:,tgtidx) = tgt.vel(:,tgtidx) + tgt.acc(:,tgtidx) * timing_struct.PRI;
|
||
end
|
||
|
||
% RF out
|
||
powVar_dB = (Ptx_dBm - 30) * ones(Nt, Nr); % -30 : dBm -> dBW
|
||
for txidx = 1 : Nt
|
||
tgt_loc_vec = tgt.pos(:, tgtidx) - (posR + txarray_loc(:, txidx));
|
||
[tgt_azi_dod, tgt_elev_dod, tgt_rng_dod] = cart2sph(tgt_loc_vec(1),tgt_loc_vec(2),tgt_loc_vec(3));
|
||
target_prop_time_dod = rng2time(tgt_rng_dod, c0);
|
||
tx_ant_gain_azi = interp1(azi_ant_pat(:,1), azi_ant_pat(:,2), rad2deg(tgt_azi_dod));
|
||
tx_ant_gain_elev_reduction = max(elev_ant_pat(:,2)) - interp1(elev_ant_pat(:,1), elev_ant_pat(:,2), rad2deg(tgt_elev_dod));
|
||
tx_ant_gain_db = tx_ant_gain_azi - tx_ant_gain_elev_reduction;
|
||
|
||
% Tx ant gain
|
||
powVar_dB(txidx,:) = powVar_dB(txidx,:) + tx_ant_gain_db;
|
||
|
||
% Freespace loss (Radar to Target)
|
||
powVar_dB(txidx,:) = powVar_dB(txidx,:) + pow2db(1/(4*pi*tgt_rng_dod^2));
|
||
|
||
% RCS
|
||
powVar_dB(txidx,:) = powVar_dB(txidx,:) + tgt.rcs(tgtidx);
|
||
|
||
for rxidx = 1 : Nr
|
||
tgt_loc_vec = tgt.pos(:, tgtidx) - (posR + rxarray_loc(:, rxidx));
|
||
[tgt_azi_doa, tgt_elev_doa, tgt_rng_doa] = cart2sph(tgt_loc_vec(1),tgt_loc_vec(2),tgt_loc_vec(3));
|
||
target_prop_time_doa = rng2time(tgt_rng_doa, c0);
|
||
rx_ant_gain_azi = interp1(azi_ant_pat(:,1), azi_ant_pat(:,2), rad2deg(tgt_azi_doa));
|
||
rx_ant_gain_elev_reduction = max(elev_ant_pat(:,2)) - interp1(elev_ant_pat(:,1), elev_ant_pat(:,2), rad2deg(tgt_elev_doa));
|
||
rx_ant_gain_db = rx_ant_gain_azi - rx_ant_gain_elev_reduction;
|
||
|
||
% Freespace loss (Target to Radar)
|
||
powVar_dB(txidx, rxidx) = powVar_dB(txidx, rxidx) + pow2db(1/(4*pi*tgt_rng_doa^2));
|
||
|
||
% Rx ant gain
|
||
powVar_dB(txidx, rxidx) = powVar_dB(txidx, rxidx) + rx_ant_gain_db + pow2db(lambda_c^2/(4*pi));
|
||
|
||
target_prop_time = target_prop_time_dod + target_prop_time_doa;
|
||
% [24.08.22] Lowpass filter 구현 때문에 2배 oversampling 함.
|
||
%delayed_wave = gen_fastramp_fmcw2(timing_struct, 2*fs, fc-TxBw/2, TxBw, m, NumChirps, target_prop_time, tx_phase_initial, opt_waveform);
|
||
delayed_wave = gen_fastramp_fmcw(timing_struct, 2*fs, fc-TxBw/2, TxBw, NumChirps, target_prop_time, tx_phase_initial, opt_waveform);
|
||
delayed_Tx_sig(:,m,rxidx) = delayed_Tx_sig(:,m,rxidx) + sqrt(db2pow(powVar_dB(txidx, rxidx))) * (delayed_wave.waveform).' .* exp(-1i * 2 * pi * DDMA_freq(txidx) * ((m-1) * delayed_wave.T_chirp));
|
||
end
|
||
end
|
||
end
|
||
end
|
||
|
||
|
||
%% Antenna (Rx)
|
||
rxsig_wo_n = delayed_Tx_sig;
|
||
|
||
%% Receiver
|
||
% Compute the cascaded noise figure and total gain of a receiver system. The system has seven stages, with these values:
|
||
% 01. LNA with a noise figure of 1.0 dB and a gain of 15.0 dB
|
||
% 02. Mixer with a noise figure of 5.0 dB and a gain of –7.0 dB
|
||
% 03. HPF1 with a noise figure of 0.5 dB and a gain of –0.5 dB
|
||
% 04. IF VGA1 with a noise figure of 0.6 dB and a gain of –15.0 dB
|
||
% 05. HPF2 with a noise figure of 0.5 dB and a gain of –0.5 dB
|
||
% 06. IF VGA2 with a noise figure of 0.6 dB and a gain of –15.0 dB
|
||
% 07. LPF1 with a noise figure of 0.6 dB and a gain of –15.0 dB
|
||
% 08. Buffer with a noise figure of 1.0 dB and a gain of –1.0 dB
|
||
% 09. LPF2 with a noise figure of 0.6 dB and a gain of –15.0 dB
|
||
% 10. ADC with Vpp : 1.2 [V], Bit : 12 [bit]
|
||
|
||
% 00. Receiver Parameters
|
||
NF_dB_LNA = 11.5;
|
||
NF_dB_Mixer = 12;
|
||
NF_dB_HPF1 = 0.5;
|
||
NF_dB_VGA1 = 20;
|
||
NF_dB_HPF2 = 0.5;
|
||
NF_dB_VGA2 = 20;
|
||
NF_dB_LPF1 = 0.5;
|
||
NF_dB_Buffer = 30;
|
||
NF_dB_LPF2 = 0.5;
|
||
|
||
G_dB_LNA = 20;
|
||
G_dB_Mixer = 0;
|
||
G_dB_HPF1 = -0.5;
|
||
G_dB_VGA1 = 13.5;
|
||
G_dB_HPF2 = -0.5;
|
||
G_dB_VGA2 = 13.5;
|
||
G_dB_LPF1 = -0.5;
|
||
G_dB_Buffer = 0;
|
||
G_dB_LPF2 = -0.5;
|
||
|
||
|
||
nf = [NF_dB_LNA NF_dB_Mixer NF_dB_HPF1 NF_dB_VGA1 NF_dB_HPF2, NF_dB_VGA2 NF_dB_LPF1 NF_dB_Buffer NF_dB_LPF2];
|
||
g = [G_dB_LNA G_dB_Mixer G_dB_HPF1 G_dB_VGA1 G_dB_HPF2, G_dB_VGA2 G_dB_LPF1 G_dB_Buffer G_dB_LPF2];
|
||
|
||
|
||
% nf = [4.0 0.5 5.0 1.0 0.6 1.0 6.0];
|
||
% g = [20.0 -0.5 -7.0 -1.0 21.75 21.75 -5.0];
|
||
|
||
[cnf,ng] = noisefigure(nf,g);
|
||
% NF_dB2 = cal_nf(nf, g, T0);
|
||
NF_dB = 12;
|
||
|
||
% 01. LNA
|
||
NF_dB_LNA = 12;
|
||
if opt_noise == 0
|
||
Pn = 0;
|
||
elseif opt_noise == 1
|
||
Pn = kb * T0 * fs * db2pow(NF_dB_LNA);
|
||
end
|
||
|
||
thermal_noise_LNA = sqrt(Pn) * randn([size(rxsig_wo_n)]);
|
||
|
||
% 02. Mixing
|
||
% 2-1. Product
|
||
mixed_sig = zeros(size(rxsig_wo_n));
|
||
mixed_noise = zeros(size(thermal_noise_LNA));
|
||
for rxidx = 1 : Nr
|
||
mixed_sig(:,:,rxidx) = dechirping(rxsig_wo_n(:,:,rxidx), Ref_wave(:,:,rxidx));
|
||
mixed_noise(:,:,rxidx) = dechirping(thermal_noise_LNA(:,:,rxidx), Ref_wave(:,:,rxidx));
|
||
end
|
||
|
||
% 2-2. High Pass Filter
|
||
% Second-order, -6 dB, (100 kHz ~ 3.2 MHz) : NXP (tef82xx)
|
||
f_cut_hpf = 1.6e6;
|
||
|
||
[zb1, pb1, kb1] = butter(1, 2*pi*f_cut_hpf, 'high', 's');
|
||
[zd1, pd1, kd1] = bilinear(zb1, pb1, kb1, 2*fs);
|
||
[hpf1_d_num, hpf1_d_den] = zp2tf(zd1, pd1, kd1);
|
||
|
||
temp_mixed_sig_hpf = zeros(size(mixed_sig));
|
||
mixed_sig_hpf = zeros(size(mixed_sig));
|
||
temp_mixed_noise_hpf = zeros(size(mixed_sig));
|
||
mixed_noise_hpf = zeros(size(mixed_sig));
|
||
for rxidx = 1 : Nr
|
||
if opt_filter == 1
|
||
temp_mixed_sig_hpf(:,:,rxidx) = filter(hpf1_d_num, hpf1_d_den, mixed_sig(:,:,rxidx));
|
||
mixed_sig_hpf(:,:,rxidx) = filter(hpf1_d_num, hpf1_d_den, temp_mixed_sig_hpf(:,:,rxidx));
|
||
|
||
temp_mixed_noise_hpf(:,:,rxidx) = filter(hpf1_d_num, hpf1_d_den, mixed_noise(:,:,rxidx));
|
||
mixed_noise_hpf(:,:,rxidx) = filter(hpf1_d_num, hpf1_d_den, temp_mixed_noise_hpf(:,:,rxidx));
|
||
else
|
||
mixed_sig_hpf(:,:,rxidx) = mixed_sig(:,:,rxidx);
|
||
|
||
mixed_noise_hpf(:,:,rxidx) = mixed_noise(:,:,rxidx);
|
||
end
|
||
end
|
||
[h_hpf1, ~] = freqz(hpf1_d_num,hpf1_d_den, rngFFTlen);
|
||
Loss_hpf1_dB = pow2db(sum(abs(h_hpf1).^2)/rngFFTlen);
|
||
Loss_hpf2_dB = pow2db(sum(abs(h_hpf1).^2)/rngFFTlen);
|
||
|
||
|
||
|
||
% 2-3. Low Pass Filter
|
||
% 1) Third-order, -6 dB, (12.5 MHz ~ 25 MHz) : NXP (tef82xx)
|
||
% 2) Third-order, -6 dB, > 40 MHz(Wide-bandwidth mode) : NXP (tef82xx)
|
||
f_cut_lpf = 12.5e6; % (20e6 * 416/512)
|
||
|
||
[zl1, pl1, kl1] = butter(1, 2*pi*f_cut_lpf, 'low', 's');
|
||
[zdl1, pdl1, kdl1] = bilinear(zl1, pl1, kl1, 2*fs);
|
||
[lpf1_d_num, lpf1_d_den] = zp2tf(zdl1, pdl1, kdl1);
|
||
|
||
[zl2, pl2, kl2] = butter(2, 2*pi*f_cut_lpf, 'low', 's');
|
||
[zdl2, pdl2, kdl2] = bilinear(zl2, pl2, kl2, 2*fs);
|
||
[lpf2_d_num, lpf2_d_den] = zp2tf(zdl2, pdl2, kdl2);
|
||
|
||
temp_mixed_sig_lpf = zeros(size(mixed_sig));
|
||
mixed_sig_lpf = zeros(size(mixed_sig));
|
||
temp_mixed_noise_lpf = zeros(size(mixed_sig));
|
||
mixed_noise_lpf = zeros(size(mixed_sig));
|
||
for rxidx = 1 : Nr
|
||
if opt_filter == 1
|
||
temp_mixed_sig_lpf(:,:,rxidx) = filter(lpf1_d_num, lpf1_d_den, mixed_sig_hpf(:,:,rxidx));
|
||
mixed_sig_lpf(:,:,rxidx) = filter(lpf2_d_num, lpf2_d_den, temp_mixed_sig_lpf(:,:,rxidx));
|
||
|
||
temp_mixed_noise_lpf(:,:,rxidx) = filter(lpf1_d_num, lpf1_d_den, mixed_noise_hpf(:,:,rxidx));
|
||
mixed_noise_lpf(:,:,rxidx) = filter(lpf2_d_num, lpf2_d_den, temp_mixed_noise_lpf(:,:,rxidx));
|
||
else
|
||
mixed_sig_lpf(:,:,rxidx) = mixed_sig_hpf(:,:,rxidx);
|
||
|
||
mixed_noise_lpf(:,:,rxidx) = mixed_noise_hpf(:,:,rxidx);
|
||
end
|
||
end
|
||
[h_lpf1, ~] = freqz(lpf1_d_num, lpf1_d_den, rngFFTlen);
|
||
Loss_lpf1_dB = pow2db(sum(abs(h_lpf1).^2)/rngFFTlen);
|
||
[h_lpf2, ~] = freqz(lpf2_d_num, lpf2_d_den, rngFFTlen);
|
||
Loss_lpf2_dB = pow2db(sum(abs(h_hpf1).^2)/rngFFTlen);
|
||
|
||
% figure()
|
||
% subplot(211); plot(Tx_wave.timeline, real(mixed_sig_lpf(:,10)));
|
||
% xlabel('Time (s)'); ylabel('Amplitude (v)');
|
||
% title('dechirped FMCW signal'); axis tight;
|
||
% subplot(212); spectrogram(mixed_sig_lpf(:,10), 320, 160, 320, fs,'yaxis');
|
||
% title('dechirped FMCW signal spectrogram');
|
||
|
||
% 2-4. Receiver gain
|
||
Rxgain_dB = 45;
|
||
mixed_sig_lpf_rxgain = mixed_sig_lpf(1:2:end,:,:) * db2pow((Rxgain_dB - (Loss_hpf1_dB + Loss_hpf1_dB + Loss_lpf1_dB + Loss_lpf2_dB))/2);
|
||
mixed_noise_lpf_rxgain = mixed_noise_lpf(1:2:end,:,:) * db2pow((Rxgain_dB - (Loss_hpf1_dB + Loss_hpf1_dB + Loss_lpf1_dB + Loss_lpf2_dB))/2);
|
||
|
||
% 2-5 Quantization Noise
|
||
Vpp_ADC = 1.2;
|
||
ADC_bit = 12;
|
||
ADC_impd = 50;
|
||
Vq_rms = Vpp_ADC/2^ADC_bit / sqrt(12);
|
||
qt_pow = Vq_rms^2 / ADC_impd;
|
||
qt_noise = sqrt(qt_pow) * randn(size(mixed_sig_lpf_rxgain));
|
||
|
||
sig_adc_out = mixed_sig_lpf_rxgain;
|
||
noise_adc_out = mixed_noise_lpf_rxgain;
|
||
|
||
%% 03. FFT
|
||
% 3-1 Range Windowing
|
||
if opt_window == 1
|
||
winRng = repmat(hann(rngFFTlen), [1,NumChirps, Nr]);
|
||
else
|
||
winRng = repmat(ones(NumSamples, 1), [1,NumChirps, Nr]);
|
||
end
|
||
|
||
scalwinRng = sum(winRng(:,1,1)) / length(winRng(:,1,1));
|
||
winR_mixed_sig_lpf = sig_adc_out .* winRng;
|
||
winR_mixed_noise_lpf = noise_adc_out .* winRng;
|
||
|
||
% 3-2. Range FFT
|
||
sig_rng_fft = fft(winR_mixed_sig_lpf, rngFFTlen, 1);
|
||
noise_rng_fft = fft(winR_mixed_noise_lpf, rngFFTlen, 1);
|
||
|
||
% 3-3 Doppler Windowing
|
||
if opt_window == 1
|
||
winDop = repmat(hann(velFFTlen).', [rngFFTlen, 1, Nr]);
|
||
else
|
||
winDop = repmat(ones(1, NumChirps), [rngFFTlen, 1, Nr]);
|
||
end
|
||
|
||
scalwinDop = sum(winDop(1,:,1)) / length(winDop(1,:,1));
|
||
winD_sig_rng_fft = sig_rng_fft .* winDop;
|
||
winD_noise_rng_fft = noise_rng_fft .* winDop;
|
||
|
||
% 3-4. Doppler FFT
|
||
sig_rng_dop_fft2_wo_n = fft(winD_sig_rng_fft, velFFTlen, 2);
|
||
noise_rng_dop_fft2 = fft(winD_noise_rng_fft, velFFTlen, 2);
|
||
|
||
%noise_att =[0.2120 0.2143 0.2206 0.2276 0.2342 0.2424 0.2468 ];
|
||
noise_att = ones(1, maxrngFFTidx);
|
||
noise_rng_dop_fft2(1:maxrngFFTidx, :, :) = noise_rng_dop_fft2(1:maxrngFFTidx, : ,:) .* repmat(noise_att', 1, 256, 4);
|
||
|
||
sig_rng_dop_fft2 = sig_rng_dop_fft2_wo_n + noise_rng_dop_fft2;
|
||
|
||
|
||
|
||
|
||
%% Signal Processing
|
||
rng_grid = beat2rng(gen_freqgrid(rngFFTlen, fs, 0), Tx_wave.f_slope, c0);
|
||
spd_grid = 0.5*dop2spd(gen_freqgrid(velFFTlen, 1/(Tx_wave.T_chirp), 1), lambda_c);
|
||
spd_grid_2 = fftshift(spd_grid);
|
||
|
||
PS_sig_single_ch = abs(fftshift(sig_rng_dop_fft2 ,2)).^2/NumChirps/NumSamples/rngFFTlen/velFFTlen;
|
||
PS_noise_single_ch = abs(fftshift(noise_rng_dop_fft2 ,2)).^2/NumChirps/NumSamples/rngFFTlen/velFFTlen;
|
||
|
||
Winloss_dB = 0;
|
||
RBW_noise = fs/velFFTlen/rngFFTlen;
|
||
%noise_floor_dBm = pow2db(kb * T0) + NF_dB + pow2db(RBW_noise) + Rxgain_dB + Winloss_dB + pow2db(sum(winDop(1,:,1).^2)/velFFTlen) + pow2db(sum(winRng(:,1,1).^2)/rngFFTlen) + 30;
|
||
%noise_floor_dBm = pow2db(kb * T0) + NF_dB + pow2db(RBW_noise) + Rxgain_dB + Winloss_dB + pow2db(1/4) - 1.76*2 + 30;
|
||
noise_floor_dBm = pow2db(kb * T0) + NF_dB + pow2db(RBW_noise) + Rxgain_dB + pow2db(sum(winDop(1,:,1).^2)/velFFTlen) + pow2db(sum(winRng(:,1,1).^2)/rngFFTlen) + 30;
|
||
|
||
% Mean noise power in single channel
|
||
np_true_single_ch = sum(PS_noise_single_ch,2) / NumChirps;
|
||
|
||
np_esti_single_ch = zeros(maxrngFFTidx, Nr);
|
||
for ch_idx = 1 : Nr
|
||
for rng_idx = 1 : maxrngFFTidx
|
||
[find_noise_index, l, u, np_esti_single_ch(rng_idx, ch_idx)] = isoutlier(PS_sig_single_ch(rng_idx, :, ch_idx), "percentiles", [10 80]);
|
||
end
|
||
end
|
||
|
||
|
||
if flag_plot_RVM == 1
|
||
figure
|
||
for ch_idx = 1 : Nr
|
||
subplot(sqrt(Nr), sqrt(Nr), ch_idx);
|
||
mesh(spd_grid, rng_grid(1:maxrngFFTidx), pow2db(squeeze(PS_sig_single_ch(1:maxrngFFTidx, :, ch_idx))) + 30);
|
||
xlabel('Rng [m]');
|
||
ylabel('Vel [m/s]');
|
||
zlabel('Pow [dBm]');
|
||
hold on
|
||
mesh(spd_grid, rng_grid(1:maxrngFFTidx), noise_floor_dBm*ones(maxrngFFTidx, velFFTlen), 'EdgeColor', 'r');
|
||
mesh(spd_grid, rng_grid(1:maxrngFFTidx), pow2db(np_esti_single_ch(:, ch_idx).*ones(1,velFFTlen)) + 30, 'EdgeColor', 'g');
|
||
legend('Received Data', 'Expected Noise Floor', 'Esti. Noise Floor', 'Location', 'best');
|
||
hold off
|
||
title(['Ch : ', num2str(ch_idx)]);
|
||
end
|
||
end
|
||
|
||
% figure
|
||
% imagesc(spd_grid, rng_grid(1:416), (pow2db(abs(sig_fft2(1:416,:,1)))))
|
||
% xlabel('Speed (m/s)'); ylabel('Range (m)'); title('Range Velocity Map');
|
||
% axis([-vel_max vel_max 0 rng_max])
|
||
% colorbar
|
||
|
||
% 1. Generate RV Quarter Matrix
|
||
% 1-1. NCI Rx
|
||
sig_Rx_NCI = sum(abs(sig_rng_dop_fft2).^2, 3) / Nr;
|
||
noise_Rx_NCI = sum(abs(noise_rng_dop_fft2).^2, 3) / Nr;
|
||
|
||
PS_sig_Rx_NCI = sig_Rx_NCI/NumChirps/NumSamples/rngFFTlen/velFFTlen;
|
||
PS_noise_Rx_NCI = noise_Rx_NCI/NumChirps/NumSamples/rngFFTlen/velFFTlen;
|
||
|
||
if flag_plot_RVM == 1
|
||
figure
|
||
mesh(spd_grid, rng_grid(1:maxrngFFTidx), pow2db(fftshift(PS_sig_Rx_NCI(1:maxrngFFTidx,:), 2)) + 30);
|
||
axis([-vel_max vel_max 0 rng_max]);
|
||
xlabel('Speed [m/s]');
|
||
ylabel('Range [m]');
|
||
zlabel('Power [dBm]');
|
||
title('RVM - NCI w.r.t Rx');
|
||
end
|
||
|
||
% NCI Tx
|
||
sig_TxRx_NCI = (1/4) * (sig_Rx_NCI(:,1:velFFTlen/4) + sig_Rx_NCI(:,velFFTlen/4 + 1 : 2*velFFTlen/4) + sig_Rx_NCI(:,2*velFFTlen/4 + 1 : 3*velFFTlen/4 ) + sig_Rx_NCI(:,3*velFFTlen/4+1 : end));
|
||
noise_TxRx_NCI = (1/4) * (noise_Rx_NCI(:,1:velFFTlen/4) + noise_Rx_NCI(:,velFFTlen/4 + 1 : 2*velFFTlen/4) + noise_Rx_NCI(:,2*velFFTlen/4 + 1 : 3*velFFTlen/4 ) + noise_Rx_NCI(:,3*velFFTlen/4+1 : end));
|
||
|
||
PS_sig_TxRx_NCI = sig_TxRx_NCI/NumChirps/NumSamples/rngFFTlen/velFFTlen;
|
||
PS_noise_TxRx_NCI = noise_TxRx_NCI/NumChirps/NumSamples/rngFFTlen/velFFTlen;
|
||
|
||
if flag_plot_RVM == 1
|
||
figure
|
||
mesh(pow2db(PS_sig_TxRx_NCI(1:maxrngFFTidx,:)) + 30);
|
||
xlabel('Dop idx');
|
||
ylabel('Rng idx');
|
||
zlabel('Power [dBm]');
|
||
title('RVM - NCI w.r.t Tx and Rx');
|
||
end
|
||
|
||
% 2. Noise Estimation
|
||
np_esti_TxRx_NCI = zeros(maxrngFFTidx, 1);
|
||
for rng_idx = 1 : maxrngFFTidx
|
||
[cal_noise_index, l, u, np_esti_TxRx_NCI(rng_idx)] = isoutlier(PS_sig_TxRx_NCI(rng_idx, :), "percentiles", [10 60]);
|
||
end
|
||
|
||
np_true_TxRx_NCI = sum(PS_noise_TxRx_NCI(1:maxrngFFTidx,:), 2) / (NumChirps/(length(DDMA_idx)+1));
|
||
|
||
if flag_plot_noise_esti == 1
|
||
figure
|
||
plot(pow2db(np_true_single_ch(1:maxrngFFTidx)) + 30);
|
||
hold on
|
||
plot(pow2db(np_esti_single_ch) + 30);
|
||
plot(pow2db(np_true_TxRx_NCI(1:maxrngFFTidx)) + 30);
|
||
plot(pow2db(np_esti_TxRx_NCI) + 30);
|
||
legend('NP ture in single ch.', 'NP esti in single ch.', 'NP true at TxRxNCI', 'NP esti at TxRxNCI');
|
||
zlabel('Power [dBm]');
|
||
end
|
||
|
||
% 3. Detection
|
||
% 3-1. Setting Threshold
|
||
snr_th = db2pow(13);
|
||
|
||
% 3-2. Find det index
|
||
[det_index_peak, b] = findpeak2D(PS_sig_TxRx_NCI(1:maxrngFFTidx,:), []);
|
||
|
||
det_index_th = [];
|
||
for rngidx = 1 : maxrngFFTidx
|
||
det_index_th = [det_index_th ; PS_sig_TxRx_NCI(rngidx, :) > np_esti_TxRx_NCI(rngidx) * snr_th];
|
||
end
|
||
|
||
det_index = det_index_th .* det_index_peak;
|
||
|
||
[rngidx_set, dopidx_set] = find(det_index>0);
|
||
|
||
for det_idx = 1 : length(rngidx_set)
|
||
temp_DET_set{det_idx}.r_m = rng_grid(rngidx_set(det_idx));
|
||
temp_DET_set{det_idx}.snr_db = pow2db((PS_sig_TxRx_NCI(rngidx_set(det_idx), dopidx_set(det_idx))) / np_esti_TxRx_NCI(rngidx_set(det_idx)));
|
||
end
|
||
|
||
% 4. DOA estimation
|
||
% 4-1. Generate Snapshot and Resolving Doppler ambiguity by DDMA
|
||
adddopidx = [64 128 192 0];
|
||
snapshot_idx_set = [4 1 2 3; 3 4 1 2; 2 3 4 1; 1 2 3 4];
|
||
temp_snapshot_data = zeros((Nt+1)*Nr,length(rngidx_set));
|
||
for idx = 1 : length(rngidx_set)
|
||
temp_snapshot = zeros((Nt+1)*Nr,1);
|
||
for txidx = 1 : Nt+1
|
||
temp_snapshot(Nr*(txidx-1)+1 : Nr*txidx) = squeeze(sig_rng_dop_fft2_wo_n(rngidx_set(idx), dopidx_set(idx)+64*(txidx-1), :));
|
||
end
|
||
|
||
% Find virtual channel by selecting min power channel
|
||
rxchpow = sum(abs(reshape(temp_snapshot, Nt+1, Nr)).^2);
|
||
[minval, minloc] = min(rxchpow);
|
||
|
||
dopidx_set(idx) = dopidx_set(idx) + adddopidx(minloc);
|
||
|
||
for txidx = 1 : Nt+1
|
||
temp_snapshot_data(Nr*(txidx-1)+1 : Nr*txidx, idx) = temp_snapshot(4*snapshot_idx_set(minloc, txidx)-3 : 4*snapshot_idx_set(minloc, txidx));
|
||
end
|
||
end
|
||
snapshot_data = temp_snapshot_data(1:Nt*Nr, :);
|
||
|
||
% Add resolved vel info to DET_set structure
|
||
for det_idx = 1 : length(rngidx_set)
|
||
temp_DET_set{det_idx}.v_amb_mps = spd_grid_2(dopidx_set(det_idx));
|
||
end
|
||
|
||
% 4-2. DOA estimation
|
||
az_array_pos = mimoarray(2, :) / u_azi;
|
||
|
||
elv_Unit = u_elv / lambda_c;
|
||
|
||
test_ang = -90 : 0.1 : 90;
|
||
svmat = exp(-1i * 2 * pi / lambda_c * az_array_pos.' * u_azi .* sind(test_ang));
|
||
|
||
det_jdx = 1;
|
||
for det_idx = 1 : size(snapshot_data, 2)
|
||
% DOA estimation
|
||
% Elv esti.
|
||
elv_phase_diff = conj(snapshot_data(4, det_idx)) * snapshot_data(9, det_idx);
|
||
temp_esti_elv_deg = asind(angle(elv_phase_diff)/2/pi/elv_Unit);
|
||
|
||
% Azi esti.
|
||
azi_spectrum = abs(svmat' * snapshot_data(:, det_idx)).^2;
|
||
[pks, esti_azi_deg] = findpeaks(azi_spectrum/max(azi_spectrum), test_ang, 'MinPeakHeight', 0.5);
|
||
|
||
esti_elv_deg = temp_esti_elv_deg * ones(1, length(esti_azi_deg));
|
||
|
||
|
||
% Construct DET set structure
|
||
for angidx = 1 : length(esti_azi_deg)
|
||
DET_set{det_jdx} = temp_DET_set{det_idx};
|
||
DET_set{det_jdx}.az_deg = esti_azi_deg(angidx);
|
||
DET_set{det_jdx}.el_deg = esti_elv_deg(angidx);
|
||
DET_set{det_jdx}.x_m = DET_set{det_jdx}.r_m * cosd(DET_set{det_jdx}.az_deg);
|
||
DET_set{det_jdx}.y_m = DET_set{det_jdx}.r_m * sind(DET_set{det_jdx}.az_deg);
|
||
DET_set{det_jdx}.z_m = DET_set{det_jdx}.r_m * sind(DET_set{det_jdx}.el_deg);
|
||
det_jdx = det_jdx + 1;
|
||
end
|
||
end
|
||
|
||
%% Results
|
||
if flag_plot_birdview == 1
|
||
figure;
|
||
x_m_set = cellfun(@(x) x.x_m, DET_set);
|
||
y_m_set = cellfun(@(x) x.y_m, DET_set);
|
||
z_m_set = cellfun(@(x) x.z_m, DET_set);
|
||
r_m_set = cellfun(@(x) x.r_m, DET_set);
|
||
v_amb_mps_set = cellfun(@(x) x.v_amb_mps, DET_set);
|
||
az_deg_set = cellfun(@(x) x.az_deg, DET_set);
|
||
el_deg_set = cellfun(@(x) x.el_deg, DET_set);
|
||
snr_db_set = cellfun(@(x) x.snr_db, DET_set);
|
||
plotdet = plot3(y_m_set, x_m_set, z_m_set, 'bo');
|
||
fn_add_data_tip(plotdet, 'Rng [m]:', r_m_set, 1);
|
||
fn_add_data_tip(plotdet, 'aVel [m/s]:', v_amb_mps_set, 2);
|
||
fn_add_data_tip(plotdet, 'Azi [deg]:', az_deg_set, 3);
|
||
fn_add_data_tip(plotdet, 'Elv [deg]:', el_deg_set, 4);
|
||
fn_add_data_tip(plotdet, 'SNR [dB]:', snr_db_set, 5);
|
||
fn_add_data_tip(plotdet, 'X [m]:', x_m_set, 6);
|
||
fn_add_data_tip(plotdet, 'Y [m]:', y_m_set, 7);
|
||
fn_add_data_tip(plotdet, 'Z [m]:', z_m_set, 8);
|
||
view([0, 90]);
|
||
grid on
|
||
xlabel('Y [m]')
|
||
ylabel('X [m]')
|
||
zlabel('Z [m]')
|
||
xlim([-rng_max, rng_max])
|
||
ylim([0 rng_max])
|
||
title('BirdView')
|
||
end
|
||
|
||
|
||
|