From 2c0a509982a8fe8c87f8608293cfafe30620d732 Mon Sep 17 00:00:00 2001 From: yeokwanggoo-stack Date: Sun, 17 Aug 2025 20:40:07 +0900 Subject: [PATCH 1/2] [25.08.17] Fix angle error: 1) Fix snapshot generating 2) DDMA resolving --- System_Analysis_Main.asv | 669 +++++++++++++++++++++++++++++++++++++++ System_Analysis_Main.m | 64 ++-- 2 files changed, 708 insertions(+), 25 deletions(-) create mode 100644 System_Analysis_Main.asv diff --git a/System_Analysis_Main.asv b/System_Analysis_Main.asv new file mode 100644 index 0000000..b9e45c7 --- /dev/null +++ b/System_Analysis_Main.asv @@ -0,0 +1,669 @@ +clear variables; +close all; +clc + +addpath(genpath(fullfile(pwd, 'Functions'))); +addpath(genpath(fullfile(pwd, 'Antenna_Pattern'))); + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% 좌표계 : +% 1) 중심(0,0,0) : 레이더의 중심 +% 2) 레이더 boresight(레이더 안테나면의 수직방향) : +x axis +% 3) 레이더 boresight 기준 왼쪽 : +y axis, 오른쪽 : -y axis +% 4) X-Y 평면에서 위로 수직 : +z axis, 아래로 수직 : -z axis +% 5) Elevation : x-y 평면 0 deg 기준. +z 방향으로 (+), -z 방향으로 (-) +% 6) Azimuth : +x axis 를 0 deg 기준. x-y 평면상에서 반시계 : (+), 시계 : (-) +% +% +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +%% 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 = 1; +flag_plot_virtual_array = 1; +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 = [100]; +tgt_azi = [-40]; +tgt_elv = [-5]; +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 = [7 0 0].'; +tgt_acc = [0 0 0].'; +tgt_rcs = [10]; +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) = sum((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]', 'radi 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 = 1e-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 함. +temp_Ref_wave = zeros(2*NumSamples, NumChirps); +for m = 1 : NumChirps + Tx_wave = gen_fastramp_fmcw(timing_struct, 2*fs, fc-TxBw/2, TxBw, NumChirps, 0, tx_phase_initial, opt_waveform); + temp_Ref_wave(:,m) = Tx_wave.waveform; +end +Ref_wave = repmat(temp_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_fmcw(timing_struct, 2*fs, fc-TxBw/2, TxBw, NumChirps, target_prop_time, tx_phase_initial, opt_waveform); + delayed_Tx_sig(:,m,rxidx) = txidx * 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_fftshift = 0.5*dop2spd(gen_freqgrid(velFFTlen, 1/(Tx_wave.T_chirp), 1), lambda_c); +spd_grid_no_fftshift = fftshift(spd_grid_fftshift); + +% 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; + +PS_sig_single_ch = abs(sig_rng_dop_fft2).^2/NumChirps/NumSamples/rngFFTlen/velFFTlen; +PS_noise_single_ch = abs(noise_rng_dop_fft2).^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_fftshift, rng_grid(1:maxrngFFTidx), pow2db((squeeze(PS_sig_single_ch(1:maxrngFFTidx, :, ch_idx)))) + 30); + xlabel('Vel [m/s]'); + ylabel('Rng [m]'); + zlabel('Pow [dBm]'); + hold on + mesh(spd_grid_fftshift, rng_grid(1:maxrngFFTidx), noise_floor_dBm*ones(maxrngFFTidx, velFFTlen), 'EdgeColor', 'r'); + mesh(spd_grid_fftshift, 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_fftshift, rng_grid(1:maxrngFFTidx), pow2db((PS_sig_Rx_NCI(1:maxrngFFTidx,:))) + 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 = [192 0 64 128]; +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_no_fftshift(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.9); + figure + plot(test_ang, pow2db(azi_spectrum)); + 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'); + set(gca, 'XDir', 'reverse'); + 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 + + + diff --git a/System_Analysis_Main.m b/System_Analysis_Main.m index d2b2bd6..6b0a5f1 100644 --- a/System_Analysis_Main.m +++ b/System_Analysis_Main.m @@ -7,8 +7,13 @@ addpath(genpath(fullfile(pwd, 'Antenna_Pattern'))); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % -% 좌표계 : ego car 기준 앞뒤 : x 축, 좌우 : y 축, 위아래 : z 축 -% +% 좌표계 : +% 1) 중심(0,0,0) : 레이더의 중심 +% 2) 레이더 boresight(레이더 안테나면의 수직방향) : +x axis +% 3) 레이더 boresight 기준 왼쪽 : +y axis, 오른쪽 : -y axis +% 4) X-Y 평면에서 위로 수직 : +z axis, 아래로 수직 : -z axis +% 5) Elevation : x-y 평면 0 deg 기준. +z 방향으로 (+), -z 방향으로 (-) +% 6) Azimuth : +x axis 를 0 deg 기준. x-y 평면상에서 반시계 : (+), 시계 : (-) % % % @@ -90,12 +95,16 @@ if flag_plot_physical_array == 1 end if flag_plot_virtual_array == 1 - figure - scatter(mimoarray(2,:) / lambda_c, mimoarray(3,:) / lambda_c, 'o') + figure + scatter(mimoarray(2,1:Nr) / lambda_c, mimoarray(3,1:Nr) / lambda_c, 'bo') + hold on + scatter(mimoarray(2,Nr+1:2*Nr) / lambda_c, mimoarray(3,Nr+1:2*Nr) / lambda_c, 'ro') + scatter(mimoarray(2,2*Nr+1:3*Nr) / lambda_c, mimoarray(3,2*Nr+1:3*Nr) / lambda_c, 'go') grid on xlabel('y-axis [\lambda]'); ylabel('z-axis [\lambda]'); title('Virtual Array Positions'); + legend('Tx1','Tx2','Tx3'); end %% Environments @@ -105,13 +114,13 @@ velR = [0 0 0].'; % radar velocity accR = [0 0 0].'; % radar acceleration % 01. Target -tgt_rng = [10]; +tgt_rng = [100]; 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_vel = [0 0 0].'; tgt_acc = [0 0 0].'; tgt_rcs = [10]; tgt_class = 'TCR'; @@ -121,11 +130,11 @@ 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))); + tgt_elv_deg(tgt_idx, 1) = rad2deg(tgt_elv_rad(tgt_idx)); + tgt_radi_vel(tgt_idx, 1) = sum((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]'}); +tgt_table_truth = table([tgt_rng tgt_radi_vel tgt_azi_deg tgt_elv_deg]', 'VariableNames', {'Tgt.'}, 'RowNames', {'Rng[m]', 'radi Vel.[m/s]', 'Azi.[deg]','Elv.[deg]'}); disp(tgt_table_truth); % 02. Clutter @@ -237,7 +246,7 @@ for tgtidx = 1 : tgt.num target_prop_time = target_prop_time_dod + target_prop_time_doa; % [24.08.22] Lowpass filter 구현 때문에 2배 oversampling 함. 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)); + delayed_Tx_sig(:,m,rxidx) = txidx * 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 @@ -443,11 +452,14 @@ 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); +spd_grid_fftshift = 0.5*dop2spd(gen_freqgrid(velFFTlen, 1/(Tx_wave.T_chirp), 1), lambda_c); +spd_grid_no_fftshift = fftshift(spd_grid_fftshift); -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; +% 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; + +PS_sig_single_ch = abs(sig_rng_dop_fft2).^2/NumChirps/NumSamples/rngFFTlen/velFFTlen; +PS_noise_single_ch = abs(noise_rng_dop_fft2).^2/NumChirps/NumSamples/rngFFTlen/velFFTlen; Winloss_dB = 0; RBW_noise = fs/velFFTlen/rngFFTlen; @@ -470,13 +482,13 @@ 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]'); + mesh(spd_grid_fftshift, rng_grid(1:maxrngFFTidx), pow2db((squeeze(PS_sig_single_ch(1:maxrngFFTidx, :, ch_idx)))) + 30); + xlabel('Vel [m/s]'); + ylabel('Rng [m]'); 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'); + mesh(spd_grid_fftshift, rng_grid(1:maxrngFFTidx), noise_floor_dBm*ones(maxrngFFTidx, velFFTlen), 'EdgeColor', 'r'); + mesh(spd_grid_fftshift, 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)]); @@ -499,7 +511,7 @@ 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); + mesh(spd_grid_fftshift, rng_grid(1:maxrngFFTidx), pow2db((PS_sig_Rx_NCI(1:maxrngFFTidx,:))) + 30); axis([-vel_max vel_max 0 rng_max]); xlabel('Speed [m/s]'); ylabel('Range [m]'); @@ -565,7 +577,7 @@ end % 4. DOA estimation % 4-1. Generate Snapshot and Resolving Doppler ambiguity by DDMA -adddopidx = [64 128 192 0]; +adddopidx = [192 0 64 128]; 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) @@ -588,7 +600,7 @@ 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)); + temp_DET_set{det_idx}.v_amb_mps = spd_grid_no_fftshift(dopidx_set(det_idx)); end % 4-2. DOA estimation @@ -603,13 +615,14 @@ 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); + elv_phase_diff = conj(snapshot_data(9, det_idx)) * snapshot_data(4, 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); - + [pks, esti_azi_deg] = findpeaks(azi_spectrum/max(azi_spectrum), test_ang, 'MinPeakHeight', 0.9); + figure + plot(test_ang, pow2db(azi_spectrum)); esti_elv_deg = temp_esti_elv_deg * ones(1, length(esti_azi_deg)); @@ -637,6 +650,7 @@ if flag_plot_birdview == 1 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'); + set(gca, 'XDir', 'reverse'); 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); -- 2.54.0 From 73cd7e1a736b321b26e26042b5bf3708598b047a Mon Sep 17 00:00:00 2001 From: yeokwanggoo-stack Date: Sun, 17 Aug 2025 20:42:30 +0900 Subject: [PATCH 2/2] [25.08.17] Deleting backup file --- System_Analysis_Main.asv | 669 --------------------------------------- 1 file changed, 669 deletions(-) delete mode 100644 System_Analysis_Main.asv diff --git a/System_Analysis_Main.asv b/System_Analysis_Main.asv deleted file mode 100644 index b9e45c7..0000000 --- a/System_Analysis_Main.asv +++ /dev/null @@ -1,669 +0,0 @@ -clear variables; -close all; -clc - -addpath(genpath(fullfile(pwd, 'Functions'))); -addpath(genpath(fullfile(pwd, 'Antenna_Pattern'))); - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% -% 좌표계 : -% 1) 중심(0,0,0) : 레이더의 중심 -% 2) 레이더 boresight(레이더 안테나면의 수직방향) : +x axis -% 3) 레이더 boresight 기준 왼쪽 : +y axis, 오른쪽 : -y axis -% 4) X-Y 평면에서 위로 수직 : +z axis, 아래로 수직 : -z axis -% 5) Elevation : x-y 평면 0 deg 기준. +z 방향으로 (+), -z 방향으로 (-) -% 6) Azimuth : +x axis 를 0 deg 기준. x-y 평면상에서 반시계 : (+), 시계 : (-) -% -% -% -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -%% 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 = 1; -flag_plot_virtual_array = 1; -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 = [100]; -tgt_azi = [-40]; -tgt_elv = [-5]; -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 = [7 0 0].'; -tgt_acc = [0 0 0].'; -tgt_rcs = [10]; -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) = sum((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]', 'radi 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 = 1e-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 함. -temp_Ref_wave = zeros(2*NumSamples, NumChirps); -for m = 1 : NumChirps - Tx_wave = gen_fastramp_fmcw(timing_struct, 2*fs, fc-TxBw/2, TxBw, NumChirps, 0, tx_phase_initial, opt_waveform); - temp_Ref_wave(:,m) = Tx_wave.waveform; -end -Ref_wave = repmat(temp_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_fmcw(timing_struct, 2*fs, fc-TxBw/2, TxBw, NumChirps, target_prop_time, tx_phase_initial, opt_waveform); - delayed_Tx_sig(:,m,rxidx) = txidx * 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_fftshift = 0.5*dop2spd(gen_freqgrid(velFFTlen, 1/(Tx_wave.T_chirp), 1), lambda_c); -spd_grid_no_fftshift = fftshift(spd_grid_fftshift); - -% 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; - -PS_sig_single_ch = abs(sig_rng_dop_fft2).^2/NumChirps/NumSamples/rngFFTlen/velFFTlen; -PS_noise_single_ch = abs(noise_rng_dop_fft2).^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_fftshift, rng_grid(1:maxrngFFTidx), pow2db((squeeze(PS_sig_single_ch(1:maxrngFFTidx, :, ch_idx)))) + 30); - xlabel('Vel [m/s]'); - ylabel('Rng [m]'); - zlabel('Pow [dBm]'); - hold on - mesh(spd_grid_fftshift, rng_grid(1:maxrngFFTidx), noise_floor_dBm*ones(maxrngFFTidx, velFFTlen), 'EdgeColor', 'r'); - mesh(spd_grid_fftshift, 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_fftshift, rng_grid(1:maxrngFFTidx), pow2db((PS_sig_Rx_NCI(1:maxrngFFTidx,:))) + 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 = [192 0 64 128]; -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_no_fftshift(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.9); - figure - plot(test_ang, pow2db(azi_spectrum)); - 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'); - set(gca, 'XDir', 'reverse'); - 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 - - - -- 2.54.0