10 Commits
22 changed files with 2172 additions and 0 deletions
+31
View File
@@ -0,0 +1,31 @@
% chirpIdx , t와 fs
function tx_matrix = apply_mimo_mode(x_amp, t, Timing, NumTx, mimoMode, fs)
N_samples = length(x_amp);
tx_matrix = zeros(NumTx, N_samples);
T_chirp = Timing.IdleTime + Timing.RampEndTime;
SamplesPerChirp = round(T_chirp * fs);
% []: (1, 2, 3...)
chirp_indices = floor((0:N_samples-1) / SamplesPerChirp) + 1;
switch upper(mimoMode)
case 'TDM'
% 1,2,3,1,2,3
active_tx_seq = mod(chirp_indices - 1, NumTx) + 1;
for tx = 1:NumTx
%
tx_matrix(tx, :) = x_amp .* (active_tx_seq == tx);
end
case 'DDMA'
%
for tx = 1:NumTx
phase_shift_seq = 2 * pi * (tx - 1) * (chirp_indices - 1) / NumTx;
tx_matrix(tx, :) = x_amp .* exp(1j * phase_shift_seq);
end
otherwise
error(' MIMO . TDM DDMA를 .');
end
end
+26
View File
@@ -0,0 +1,26 @@
function [t, tx_mask] = generate_waveform_timing(RadarParams)
% RadarParams TX
% TX .
%
% :
% - RadarParams:
%
% :
% - t: ()
% - tx_mask: TX (0 1)
Timing = RadarParams.Waveform.Timing;
NumChirps = RadarParams.Waveform.NumChirps;
T_chirp = Timing.IdleTime + Timing.RampEndTime;
T_frame = T_chirp * NumChirps; %
fs = RadarParams.Waveform.fs_waveform;
N_samples = round(T_frame * fs);
% ( )
t = linspace(0, T_frame, N_samples);
% (t_mod) ( )
t_mod = mod(t, T_chirp);
% TX
tx_mask = zeros(1, N_samples);
tx_start_abs = Timing.IdleTime + Timing.TxStartTime;
tx_mask(t_mod >= tx_start_abs) = 1;
end
+25
View File
@@ -0,0 +1,25 @@
function TxOut = radiate_antenna(Target, TxPatternArray)
% : , TX
% : TxOut.G_tx_amp (, TX별 [NumTargets x NumTx])
NumTx = length(TxPatternArray);
NumTargets = Target.NumTargets;
TxOut.G_tx_amp = zeros(NumTargets, NumTx);
for k = 1:NumTargets
for tx = 1:NumTx
pat = TxPatternArray(tx);
%
if strcmpi(pat.Type, '2D')
g_dBi = interp2(pat.az_angles, pat.el_angles, pat.gain_dBi, Target.az(k), Target.el(k), 'linear', -20);
else
g_az = interp1(pat.az_angles, pat.gain_az_dBi, Target.az(k), 'linear', -20);
g_el = interp1(pat.el_angles, pat.gain_el_dBi, Target.el(k), 'linear', -20);
g_dBi = g_az + g_el - pat.max_gain_dBi;
end
% (sqrt(10^(G_dBi/10)))
TxOut.G_tx_amp(k, tx) = sqrt(10^(g_dBi/10));
end
end
end
+145
View File
@@ -0,0 +1,145 @@
function fig = visualize_multi_tx_waveform(t, Timing, fc, f_start, Slope, tx_mask, TotalNumChirps, mimoMode, NumTx)
fig = figure('Name', 'Multi-Chirp & MIMO Modulation', 'Position', [150, 150, 1100, 750]);
T_chirp = Timing.IdleTime + Timing.RampEndTime;
NumChirpsToPlot = min(TotalNumChirps, 4);
plot_idx = (t <= NumChirpsToPlot * T_chirp);
t_plot = t(plot_idx);
tx_mask_plot = tx_mask(plot_idx);
t_mod = mod(t_plot, T_chirp);
t_ramp = t_mod - Timing.IdleTime;
inst_freq_theoretical = f_start + Slope * t_ramp;
inst_freq_masked = inst_freq_theoretical;
inst_freq_masked(tx_mask_plot == 0) = NaN;
% =========================================================
% [Subplot 1] :
% =========================================================
subplot(2, 1, 1);
y_min = (f_start / 1e9) - 0.25;
y_max = (f_start / 1e9) + (Slope * Timing.RampEndTime / 1e9) + 0.1;
plot(t_plot * 1e6, inst_freq_masked / 1e9, 'b', 'LineWidth', 2);
grid on; hold on;
title(sprintf('Continuous Multi-Chirp Sequence (Showing %d of %d Chirps)', NumChirpsToPlot, TotalNumChirps), 'FontSize', 12);
ylabel('Absolute Frequency (GHz)', 'FontSize', 11);
ylim([y_min, y_max]);
if ~isempty(t_plot)
xlim([0, max(t_plot)*1e6]);
else
xlim([0, 1]); % fallback range when no data
end
f_valid_start = f_start + (Slope * Timing.AdcStartTime);
f_valid_end = f_start + (Slope * (Timing.AdcStartTime + Timing.AdcSampTime));
valid_bandwidth = f_valid_end - f_valid_start;
info_str = {
sprintf(' Valid Start Freq : %.4f GHz', f_valid_start / 1e9), ...
sprintf(' Center Freq : %.4f GHz', fc / 1e9), ...
sprintf(' Valid End Freq : %.4f GHz', f_valid_end / 1e9), ...
sprintf(' Transmit Bandwidth : %.4f GHz', valid_bandwidth / 1e9), ...
sprintf(' Total Gen Chirps : %d', TotalNumChirps)
};
text(0.02, 0.96, info_str, 'Units', 'normalized', 'FontSize', 10, 'FontWeight', 'bold', 'BackgroundColor', [1 1 1 0.85], 'EdgeColor', 'k', 'VerticalAlignment', 'top', 'Margin', 5);
guide_line_args = {'Color', [0 0 0.5], 'LineStyle', '--', 'LineWidth', 1};
for i = 0:NumChirpsToPlot
plot([i * T_chirp * 1e6, i * T_chirp * 1e6], [y_min, y_max], guide_line_args{:});
end
draw_dim_arrow = @(x1, x2, y, label_str) ...
[plot([x1, x2], [y, y], 'k-', 'LineWidth', 1.2), ...
fill([x1, x1 + min(0.6, (x2-x1)*0.35), x1 + min(0.6, (x2-x1)*0.35)], [y, y + 0.015, y - 0.015], 'k', 'EdgeColor', 'none'), ...
fill([x2, x2 - min(0.6, (x2-x1)*0.35), x2 - min(0.6, (x2-x1)*0.35)], [y, y + 0.015, y - 0.015], 'k', 'EdgeColor', 'none'), ...
text((x1+x2)/2, y + 0.03, label_str, 'HorizontalAlignment', 'center', 'VerticalAlignment', 'bottom', 'FontSize', 9, 'FontWeight', 'bold', 'BackgroundColor', 'w', 'EdgeColor', 'k')];
y_pri = y_min + 0.08;
pri_us = T_chirp * 1e6;
prf_khz = (1 / T_chirp) / 1e3;
pri_prf_label = sprintf('PRI: %.1f \\mus\nPRF: %.1f kHz', pri_us, prf_khz);
draw_dim_arrow(0, pri_us, y_pri, pri_prf_label);
hold off;
% =========================================================
% [Subplot 2] : MIMO Active TX & Phase
% =========================================================
subplot(2, 1, 2);
hold on; grid on;
% --- [ ]: TX ---
% lines() .
tx_colors = lines(NumTx);
for tx = 1:NumTx
plot([0, NumChirpsToPlot * T_chirp * 1e6], [tx, tx], ':', 'Color', [0.8 0.8 0.8], 'HandleVisibility', 'off');
end
for i = 0 : NumChirpsToPlot - 1
t_center_us = (i + 0.5) * T_chirp * 1e6;
if strcmpi(mimoMode, 'TDM')
active_tx_list = mod(i, NumTx) + 1;
elseif strcmpi(mimoMode, 'DDMA')
active_tx_list = 1:NumTx;
else
active_tx_list = 1;
end
for tx = 1:NumTx
% TX
current_color = tx_colors(tx, :);
if ismember(tx, active_tx_list)
if strcmpi(mimoMode, 'DDMA')
phase_rad = 2 * pi * (tx - 1) * i / NumTx;
phase_deg = mod(rad2deg(phase_rad), 360);
else
phase_deg = 0;
end
% :
plot(t_center_us, tx, 'o', 'MarkerSize', 12, 'MarkerFaceColor', current_color, 'MarkerEdgeColor', current_color, 'HandleVisibility', 'off');
text_offset_us = T_chirp * 1e6 * 0.08;
text(t_center_us + text_offset_us, tx, sprintf('%d^\\circ', round(phase_deg)), ...
'VerticalAlignment', 'middle', 'HorizontalAlignment', 'left', ...
'FontSize', 9, 'FontWeight', 'bold', 'Color', current_color);
else
% :
plot(t_center_us, tx, 'o', 'MarkerSize', 12, 'MarkerFaceColor', 'w', 'MarkerEdgeColor', current_color, 'HandleVisibility', 'off');
end
end
end
title(sprintf('%s Active TX Antenna & Phase Map', upper(mimoMode)), 'FontSize', 12);
ylabel('TX Antenna', 'FontSize', 11);
yticks(1:NumTx);
% Y축 ( )
% yticklabels(arrayfun(@(x) sprintf('\\color[rgb]{%f,%f,%f}TX %d', tx_colors(x,1), tx_colors(x,2), tx_colors(x,3), x), 1:NumTx, 'UniformOutput', false));
yticklabels(arrayfun(@(x) sprintf('TX %d', x), 1:NumTx, 'UniformOutput', false));
ylim([0.5, NumTx + 0.5]);
% --- [ ]: ( ) ---
neutral_color = [0.3 0.3 0.3];
h_on = plot(NaN, NaN, 'o', 'MarkerSize', 10, 'MarkerFaceColor', neutral_color, 'MarkerEdgeColor', neutral_color, 'DisplayName', 'Active');
h_off = plot(NaN, NaN, 'o', 'MarkerSize', 10, 'MarkerFaceColor', 'w', 'MarkerEdgeColor', neutral_color, 'DisplayName', 'Inactive');
legend([h_on, h_off], 'Location', 'best', 'FontSize', 10, 'Box', 'on');
xlabel('Time (\mus)', 'FontSize', 11);
xlim([0, max(t_plot)*1e6]);
y_lims_bottom = ylim;
for i = 0:NumChirpsToPlot
plot([i * T_chirp * 1e6, i * T_chirp * 1e6], y_lims_bottom, guide_line_args{:}, 'HandleVisibility', 'off');
end
hold off;
end
+78
View File
@@ -0,0 +1,78 @@
function fig = visualize_tx_waveform(t, Timing, fc, f_start, Slope, tx_mask, peak_phase_error, f_ripple)
fig = figure('Name', 'FMCW Radar Timing & True RF Frequency', 'Position', [100, 100, 1100, 700]);
t_ramp = t - Timing.IdleTime;
%
inst_freq_theoretical = f_start + Slope * t_ramp;
freq_nonlin = peak_phase_error * f_ripple * cos(2 * pi * f_ripple * t_ramp);
inst_freq_theoretical = inst_freq_theoretical + freq_nonlin;
inst_freq_masked = inst_freq_theoretical;
inst_freq_masked(tx_mask == 0) = NaN;
% Y축
y_min = (f_start / 1e9) - 0.2;
y_max = (f_start / 1e9) + (Slope * Timing.RampEndTime / 1e9) + 0.1;
% --- [ ] ---
plot(t * 1e6, inst_freq_masked / 1e9, 'b', 'LineWidth', 2.5);
title(sprintf('Mathematical Instantaneous RF Frequency & HW Timing (Center fc = %.2f GHz)', fc/1e9), 'FontSize', 12);
xlabel('Time (\mus)', 'FontSize', 11); ylabel('Absolute Frequency (GHz)', 'FontSize', 11);
ylim([y_min, y_max]);
grid on; hold on;
% --- [ : ] ---
% ADC /
f_valid_start = f_start + (Slope * Timing.AdcStartTime);
f_valid_end = f_start + (Slope * (Timing.AdcStartTime + Timing.AdcSampTime));
valid_bandwidth = f_valid_end - f_valid_start; %
info_str = {
sprintf(' Valid Start Freq : %.4f GHz', f_valid_start / 1e9), ...
sprintf(' Center Freq : %.4f GHz', fc / 1e9), ...
sprintf(' Valid End Freq : %.4f GHz', f_valid_end / 1e9), ...
sprintf(' Valid Bandwidth : %.4f GHz', valid_bandwidth / 1e9) %
};
text(0.02, 0.96, info_str, 'Units', 'normalized', ...
'FontSize', 10, 'FontWeight', 'bold', 'BackgroundColor', [1 1 1 0.85], ...
'EdgeColor', 'k', 'VerticalAlignment', 'top', 'Margin', 5);
% --- [ ] ---
y_lims = ylim;
guide_line_args = {'Color', [0 0 0.5], 'LineStyle', '--', 'LineWidth', 1};
x_coords = [0, Timing.IdleTime, ...
(Timing.IdleTime + Timing.TxStartTime), ...
(Timing.IdleTime + Timing.AdcStartTime), ...
(Timing.IdleTime + Timing.AdcStartTime + Timing.AdcSampTime), ...
(Timing.IdleTime + Timing.RampEndTime)];
x_coords_us = x_coords * 1e6;
for i = 1:length(x_coords_us)
plot([x_coords_us(i), x_coords_us(i)], y_lims, guide_line_args{:});
end
% --- [ ] ---
draw_dim_arrow = @(x1, x2, y, name, val_us) ...
[plot([x1, x2], [y, y], 'k-', 'LineWidth', 1.2), ...
fill([x1, x1 + min(0.6, (x2-x1)*0.35), x1 + min(0.6, (x2-x1)*0.35)], [y, y + 0.015, y - 0.015], 'k', 'EdgeColor', 'none'), ...
fill([x2, x2 - min(0.6, (x2-x1)*0.35), x2 - min(0.6, (x2-x1)*0.35)], [y, y + 0.015, y - 0.015], 'k', 'EdgeColor', 'none'), ...
text((x1+x2)/2, y + 0.03, sprintf('%s\n(%.1f \\mus)', name, val_us), 'HorizontalAlignment', 'center', 'VerticalAlignment', 'bottom', 'FontSize', 9, 'FontWeight', 'bold', 'BackgroundColor', 'w', 'EdgeColor', 'k')];
%
y1 = y_min + 0.05;
y2 = y_min + 0.18;
y3 = y_min + 0.31;
draw_dim_arrow(x_coords_us(1), x_coords_us(2), y1, 'Idle Time', Timing.IdleTime * 1e6);
draw_dim_arrow(x_coords_us(2), x_coords_us(3), y2, 'TX Start', Timing.TxStartTime * 1e6);
draw_dim_arrow(x_coords_us(2), x_coords_us(4), y3, 'ADC Delay', Timing.AdcStartTime * 1e6);
draw_dim_arrow(x_coords_us(4), x_coords_us(5), y2, 'ADC Sampling (Valid)', Timing.AdcSampTime * 1e6);
draw_dim_arrow(x_coords_us(5), x_coords_us(6), y1, 'Excess', Timing.ExcessTime * 1e6);
hold off;
end
+59
View File
@@ -0,0 +1,59 @@
function ChannelOut = apply_channel_effects(RadarParams)
Target = RadarParams.Target;
TxPos = RadarParams.Antenna.TxPos;
RxPos = RadarParams.Antenna.RxPos;
NumChirps = RadarParams.Waveform.NumChirps;
Timing = RadarParams.Waveform.Timing;
fs_adc = RadarParams.Waveform.fs_adc;
f_start = RadarParams.Waveform.f_start;
Slope = RadarParams.Waveform.Slope;
c = RadarParams.Basic.c;
NumTx = size(TxPos, 2);
NumRx = size(RxPos, 2);
NumTargets = Target.NumTargets;
T_pri = RadarParams.Waveform.PRI;
t_adc = 0 : 1/fs_adc : (Timing.AdcSampTime - 1/fs_adc);
N_adc = length(t_adc);
% []
% t_adc에 (lambda)
f_inst = f_start + Slope * t_adc;
lambda_t = c ./ f_inst;
ChannelOut.tau = zeros(NumTargets, NumRx, NumTx, NumChirps, N_adc);
ChannelOut.space_loss_amp = zeros(NumTargets, NumRx, NumTx, NumChirps, N_adc);
ref_point = (TxPos(:,1) + RxPos(:,1)) / 2;
for k = 1:NumTargets
az_rad = deg2rad(Target.az(k));
el_rad = deg2rad(Target.el(k));
rel_pos0 = Target.R(k) * [cos(el_rad)*sin(az_rad); cos(el_rad)*cos(az_rad); sin(el_rad)];
target_world_pos0 = ref_point + rel_pos0;
v_vec = Target.v(k) * (rel_pos0 / norm(rel_pos0));
sigma = 10^(Target.rcs(k)/10);
for m = 0:NumChirps-1
%
target_pos_at_chirp = target_world_pos0 + v_vec * (m * T_pri);
for tx = 1:NumTx
for rx = 1:NumRx
%
curr_target_pos = target_pos_at_chirp + v_vec * t_adc;
d_tx = sqrt(sum((curr_target_pos - TxPos(:,tx)).^2, 1));
d_rx = sqrt(sum((curr_target_pos - RxPos(:,rx)).^2, 1));
% TTD
ChannelOut.tau(k, rx, tx, m+1, :) = (d_tx + d_rx) / c;
% [ ]
denom = (4*pi)^3 * (d_tx.^2 .* d_rx.^2);
ChannelOut.space_loss_amp(k, rx, tx, m+1, :) = sqrt((lambda_t.^2 * sigma) ./ denom);
end
end
end
end
end
+39
View File
@@ -0,0 +1,39 @@
function quantized_data = apply_adc_quantization(adc_raw_data, n_bits, v_full_scale, mode)
% :
% - adc_raw_data: HPF/LPF가 (Voltage Scale)
% - n_bits: ADC
% - v_full_scale: ADC Peak-to-Peak (e.g., 2.0V)
% - mode: 'IQ' (Complex ) 'Real' (Real )
v_max = v_full_scale / 2;
lsb = v_full_scale / (2^n_bits);
if strcmpi(mode, 'IQ')
% --- Case 1: I/Q Demodulation (Complex) ---
% I채널과 Q채널 Clipping Quantization
r_part = real(adc_raw_data);
i_part = imag(adc_raw_data);
r_part = max(min(r_part, v_max), -v_max);
i_part = max(min(i_part, v_max), -v_max);
q_r = round(r_part / lsb) * lsb;
q_i = round(i_part / lsb) * lsb;
quantized_data = q_r + 1j * q_i;
elseif strcmpi(mode, 'Real')
% --- Case 2: Real Demodulation ---
% (I채널) ADC
%
r_part = real(adc_raw_data);
r_part = max(min(r_part, v_max), -v_max);
q_r = round(r_part / lsb) * lsb;
% (Real)
quantized_data = q_r;
else
error(' IQ Real .');
end
end
+39
View File
@@ -0,0 +1,39 @@
function filtered_data = apply_analog_hpf(adc_raw_data, fs_adc, fc_hpf_Hz)
% :
% - adc_raw_data: [NumRx, NumTx, NumChirps, N_samples]
% - fs_adc: ADC (Hz)
% - fc_hpf_Hz: HPF (Cut-off Frequency, Hz)
[NumRx, NumTx, NumChirps, N_samples] = size(adc_raw_data);
filtered_data = zeros(size(adc_raw_data));
% 1. 1 HPF (S-plane -> Z-plane )
% Analog Transfer Function: H(s) = s / (s + omega_c)
omega_c = 2 * pi * fc_hpf_Hz;
% (Bilinear Transform)
% [b, a] = butter(1, fc_hpf_Hz / (fs_adc/2), 'high'); %
T = 1 / fs_adc;
alpha = 2 / T;
b0 = alpha / (alpha + omega_c);
b1 = -alpha / (alpha + omega_c);
a1 = (omega_c - alpha) / (alpha + omega_c);
b = [b0, b1];
a = [1, a1];
% 2.
for rx = 1:NumRx
for tx = 1:NumTx
for m = 1:NumChirps
% 퀀
raw_sig = squeeze(adc_raw_data(rx, tx, m, :));
% ( 'filtfilt' 'filter' )
% (Causal) filter
filtered_data(rx, tx, m, :) = filter(b, a, raw_sig);
end
end
end
end
+206
View File
@@ -0,0 +1,206 @@
function [adc_raw_data, adc_combined] = apply_lna_and_mixer(RxOut, TxOut, RadarParams)
% : adc_raw_data [NumRx, NumTx, NumChirps, N_adc_samples] -
% adc_combined [NumRx, NumChirps, N_adc_samples] - ADC ( TX )
% : RxOut, TxOut, RadarParams
% MIMO TX .
Timing = RadarParams.Waveform.Timing;
fs_adc = RadarParams.Waveform.fs_adc;
f_start = RadarParams.Waveform.f_start;
Slope = RadarParams.Waveform.Slope;
pn_level = RadarParams.Waveform.nonideal.pn_level;
f_ripple = RadarParams.Waveform.nonideal.f_ripple;
peak_phase_error = RadarParams.Waveform.nonideal.peak_phase_error;
enable_phase_noise = true;
enable_nonlinearity = true;
if isfield(RadarParams.Waveform.nonideal, 'enable_phase_noise')
enable_phase_noise = logical(RadarParams.Waveform.nonideal.enable_phase_noise);
end
if isfield(RadarParams.Waveform.nonideal, 'enable_nonlinearity')
enable_nonlinearity = logical(RadarParams.Waveform.nonideal.enable_nonlinearity);
end
use_datasheet_phase_noise = false;
phase_noise_cfg = struct();
if isfield(RadarParams.Waveform.nonideal, 'phaseNoise')
phase_noise_cfg = RadarParams.Waveform.nonideal.phaseNoise;
has_offset = isfield(phase_noise_cfg,'offset_Hz') && ~isempty(phase_noise_cfg.offset_Hz);
has_level = isfield(phase_noise_cfg,'level_dBc_Hz') && ~isempty(phase_noise_cfg.level_dBc_Hz);
same_length = has_offset && has_level && (numel(phase_noise_cfg.offset_Hz) == numel(phase_noise_cfg.level_dBc_Hz));
if same_length
use_datasheet_phase_noise = true;
end
end
mimoMode = RadarParams.Waveform.mimoMode;
NumTx = RadarParams.Antenna.NumTx;
rxPathGain_dB = RadarParams.Rxpath.rxPathGain_dB;
system_NF_dB = RadarParams.Rxpath.system_NF_dB;
impedance_ohm = 50; % ( 50)
PA_Profile = RadarParams.RFOutput.PA_Profile;
SpurParams = RadarParams.SpurParams;
enable_spur = true;
if isstruct(SpurParams) && isfield(SpurParams, 'enabled')
enable_spur = logical(SpurParams.enabled);
end
[NumTargets, NumRx, ~, NumChirps, N_adc_samples] = size(RxOut.space_loss_amp);
% TDM : TX별
if strcmpi(mimoMode, 'TDM') && mod(NumChirps, NumTx) ~= 0
error('TDM mode requires NumChirps to be a multiple of NumTx. Current NumChirps=%d, NumTx=%d.', NumChirps, NumTx);
end
% t=0 ( ) , ADC
t_adc = Timing.AdcStartTime : 1/fs_adc : (Timing.AdcStartTime + Timing.AdcSampTime - 1/fs_adc);
f_inst = f_start + Slope * t_adc;
% --- 1. PA_Profile을 (Vtx) ---
% PA_Profile의 ( PA_Profile(:,1), PA_Profile(:,2) )
if isstruct(PA_Profile)
P_inst_dBm = interp1(PA_Profile.freqs, PA_Profile.power_dBm, f_inst, 'linear', 'extrap');
else
P_inst_dBm = interp1(PA_Profile(:,1), PA_Profile(:,2), f_inst, 'linear', 'extrap');
end
% dBm -> Watt -> (V) ( )
Vtx_inst = sqrt(10.^((P_inst_dBm - 30) / 10) * impedance_ohm);
% --- 2. ---
T_ref = RadarParams.Basic.T0;
P_noise_floor_W = RadarParams.Basic.kb * T_ref * fs_adc;
rxPathGain_lin = 10^(rxPathGain_dB/10);
system_NF_lin = 10^(system_NF_dB/10);
% ADC
P_noise_total_W = P_noise_floor_W * system_NF_lin * rxPathGain_lin;
sigma_n = sqrt(P_noise_total_W * impedance_ohm);
% --- MIMO TX ---
% TDM: TX만
% DDMA: TX +
% legacy phase-noise model ( )
window_size = 50;
ma_filter = ones(1,window_size)/window_size;
% --- 3. ADC Raw Data ---
adc_raw_data = zeros(NumRx, NumTx, NumChirps, N_adc_samples);
for m = 1:NumChirps
% generate phase noise & nonlinearity for this chirp
t_rel = t_adc; % time within chirp
if enable_phase_noise && use_datasheet_phase_noise
phi_noise = generate_phase_noise_from_datasheet(t_rel, phase_noise_cfg);
elseif enable_phase_noise
phi_noise = pn_level * filter(ma_filter, 1, randn(1, N_adc_samples));
else
phi_noise = zeros(1, N_adc_samples);
end
phi_nonlin = zeros(1, N_adc_samples);
if enable_nonlinearity
ramp_idx2 = (t_rel >= Timing.IdleTime);
phi_nonlin(ramp_idx2) = peak_phase_error * sin(2 * pi * f_ripple * (t_rel(ramp_idx2) - Timing.IdleTime));
end
% MIMO TX
if strcmpi(mimoMode, 'TDM')
% TDM: TX (1,2,3,1,2,3,...)
active_tx_list = mod(m - 1, NumTx) + 1; % : m=1TX1, m=2TX2, m=3TX1...
elseif strcmpi(mimoMode, 'DDMA')
% DDMA: TX
active_tx_list = 1:NumTx;
else
% : TX
active_tx_list = 1:NumTx;
end
for tx = 1:NumTx
% TDM TX는
if strcmpi(mimoMode, 'TDM') && ~ismember(tx, active_tx_list)
adc_raw_data(:, tx, m, :) = 0; % TX는 0
continue;
end
for rx = 1:NumRx
signal_mix = zeros(1, N_adc_samples);
for kt = 1:NumTargets
% (1) : Vtx(f) * G_tx * Space_Loss(f) * G_rx * rxPathGain
A_path = squeeze(RxOut.space_loss_amp(kt, rx, tx, m, :))';
% (Vtx_inst)
A_total = Vtx_inst .* (TxOut.G_tx_amp(kt, tx) * A_path * sqrt(rxPathGain_lin));
% (2) TTD Beat Phase ( t_adc )
tau = squeeze(RxOut.tau(kt, rx, tx, m, :))';
beat_phase = 2*pi * (f_start * tau + Slope * tau .* t_adc - 0.5 * Slope * tau.^2);
% add transmit impairments (phase noise + nonlinearity)
beat_phase = beat_phase + phi_noise + phi_nonlin;
% (3) DDMA
if strcmpi(mimoMode, 'DDMA')
phase_shift = 2 * pi * (tx - 1) * (m - 1) / NumTx;
beat_phase = beat_phase + phase_shift;
end
% (4)
base_sig = A_total .* exp(1j * beat_phase);
if enable_spur && exist('SpurParams','var') && ~isempty(SpurParams)
spur_vec = generate_spur(SpurParams, t_adc, base_sig);
signal_mix = signal_mix + base_sig + spur_vec;
else
signal_mix = signal_mix + base_sig;
end
end
% (5)
noise = (sigma_n/sqrt(2)) * (randn(1, N_adc_samples) + 1j*randn(1, N_adc_samples));
adc_raw_data(rx, tx, m, :) = signal_mix + noise;
end
end
end
% --- 4. ADC ( TX RX별로 ) ---
% RX가 TX
adc_combined = zeros(NumRx, NumChirps, N_adc_samples);
for m = 1:NumChirps
for rx = 1:NumRx
% TX의
combined_signal = squeeze(sum(adc_raw_data(rx, :, m, :), 2));
% ADC spur이 (orientation )
if enable_spur && exist('SpurParams','var') && isfield(SpurParams,'adc')
spur_adc = generate_spur(SpurParams, t_adc);
% ensure same shape as combined_signal (Nx1 vs 1xN)
spur_adc = reshape(spur_adc, size(combined_signal));
combined_signal = combined_signal + spur_adc;
end
adc_combined(rx, m, :) = combined_signal;
end
end
% --- 5. : DC offset/LO leakage clipping ---
if enable_spur && exist('SpurParams','var')
% LO leakage: /DC
if isfield(SpurParams,'lo_leak')
dc_amp = SpurParams.lo_leak.amp;
adc_combined = adc_combined + dc_amp;
end
% Clipping spur:
if isfield(SpurParams,'clip')
target_range = SpurParams.clip.range;
beat_freq = 2 * Slope * target_range / RadarParams.Basic.c;
clip_amp = SpurParams.clip.amp;
for m = 1:NumChirps
% ()
clip_tone = clip_amp * cos(2*pi*beat_freq*t_adc);
for rx = 1:NumRx
adc_combined(rx, m, :) = adc_combined(rx, m, :) + reshape(clip_tone, [1, 1, length(clip_tone)]);
end
end
end
end
end
@@ -0,0 +1,58 @@
function phi_noise = generate_phase_noise_from_datasheet(t, phaseNoiseCfg)
% DATASHEET SSB phase noise(dBc/Hz) (rad)
% :
% t : [1xN] [Nx1] ()
% phaseNoiseCfg.offset_Hz : (Hz)
% phaseNoiseCfg.level_dBc_Hz : SSB phase noise (dBc/Hz)
% :
% phi_noise : [1xN] (rad)
t = t(:);
N = numel(t);
phi_noise = zeros(1, N);
if N < 2
return;
end
dt = mean(diff(t));
fs = 1 / dt;
df = fs / N;
offsets = phaseNoiseCfg.offset_Hz(:);
levels_dBc = phaseNoiseCfg.level_dBc_Hz(:);
valid = isfinite(offsets) & isfinite(levels_dBc) & offsets > 0;
offsets = offsets(valid);
levels_dBc = levels_dBc(valid);
if isempty(offsets)
return;
end
[offsets, order] = sort(offsets, 'ascend');
levels_dBc = levels_dBc(order);
f_pos = (1:floor(N/2))' * df;
if isempty(f_pos)
return;
end
if numel(offsets) == 1
L_dBc = levels_dBc(1) * ones(size(f_pos));
else
L_dBc = interp1(log10(offsets), levels_dBc, log10(f_pos), 'linear', 'extrap');
end
% SSB phase noise L(f) PSD : L(f) ~= 0.5 * S_phi(f)
% => S_phi(f) ~= 2 * 10^(L(f)/10) [rad^2/Hz]
S_phi = 2 * 10.^(L_dBc / 10);
%
amp = sqrt(2 * S_phi * df); % (rad)
rand_phase = 2 * pi * rand(size(f_pos));
t_rel = t - t(1);
% NxK
phase_matrix = 2*pi*(t_rel * f_pos.') + rand_phase.';
phi_noise = (cos(phase_matrix) * amp).';
phi_noise = phi_noise(:).';
end
+83
View File
@@ -0,0 +1,83 @@
function spur = generate_spur(SpurParams, t, baseSignal)
%GENERATE_SPUR create spur signals based on a parameter structure
% spur = GENERATE_SPUR(SpurParams, t, baseSignal)
% t : time vector (row or column)
% baseSignal : (optional) complex baseband signal used for some spur
% mechanisms such as mixer nonlinearity. If omitted the
% function will only generate independent tones.
%
% SpurParams is a struct that may contain any of the following fields
% .lo : struct with fields amp, freq, phase (optional) -
% LO-related spur added to phase of signal
% .mixer : struct with fields alpha2, alpha3 - coefficients for
% 2nd/3rd order nonlinearity. Requires baseSignal input.
% .adc : struct with fields amp, freq, phase - tone added after
% ADC (complex). Does not require baseSignal.
% .switch : struct with fields amp, freq - square/pulse spur
% .pulse : struct with fields amp, freq, duty (01), phase -
% periodic gating/pulse train (lock) added to output
%
% Example:
% SpurParams = struct();
% SpurParams.lo = struct('amp',0.01,'freq',1e6,'phase',0);
% SpurParams.mixer= struct('alpha2',1e-4,'alpha3',1e-6);
% SpurParams.adc = struct('amp',1e-3,'freq',2e6,'phase',0);
% SpurParams.switch = struct('amp',5e-4,'freq',5e5);
% SpurParams.pulse = struct('amp',0.02,'freq',500e3,'duty',0.1,'phase',0);
% LO , , ADC ,
% 10% duty pulse train .
%
% The returned spur vector has the same dimensions as t. If baseSignal
% is supplied the mixer nonlinearity is computed elementwise using that
% signal; otherwise only independent spur terms are returned.
if nargin < 3
baseSignal = [];
end
% ensure t is row for consistent operations
t = t(:)';
spur = zeros(size(t));
if isfield(SpurParams, 'lo')
p = SpurParams.lo;
phase = 0;
if isfield(p, 'phase'); phase = p.phase; end
spur = spur + p.amp .* sin(2*pi*p.freq .* t + phase);
end
if isfield(SpurParams, 'mixer') && ~isempty(baseSignal)
m = SpurParams.mixer;
% apply polynomial nonlinearity to the provided base signal
if isfield(m, 'alpha2')
spur = spur + m.alpha2 .* (baseSignal.^2);
end
if isfield(m, 'alpha3')
spur = spur + m.alpha3 .* (baseSignal.^3);
end
end
if isfield(SpurParams, 'adc')
a = SpurParams.adc;
phase = 0;
if isfield(a, 'phase'); phase = a.phase; end
spur = spur + a.amp .* exp(1j*(2*pi*a.freq .* t + phase));
end
if isfield(SpurParams, 'switch')
s = SpurParams.switch;
spur = spur + s.amp .* square(2*pi*s.freq .* t);
end
if isfield(SpurParams, 'pulse')
p = SpurParams.pulse;
% duty default 50%
duty = 0.5;
if isfield(p,'duty'); duty = p.duty; end
phase = 0;
if isfield(p,'phase'); phase = p.phase; end
% generate pulse train (0/1) scaled by amp
% square with duty cycle multiplied then shifted to [0,1]
pulse_wave = (square(2*pi*p.freq .* t + phase, duty*100)+1)/2;
spur = spur + p.amp .* pulse_wave;
end
end
+21
View File
@@ -0,0 +1,21 @@
function RxOut = receive_antenna(ChannelOut, Target, RxPattern)
[NumTargets, NumRx, NumTx, NumChirps, N_adc] = size(ChannelOut.space_loss_amp);
RxOut = ChannelOut;
for k = 1:NumTargets
for rx = 1:NumRx
pat = RxPattern(rx);
if strcmpi(pat.Type, '2D')
g_dBi = interp2(pat.az_angles, pat.el_angles, pat.gain_dBi, Target.az(k), Target.el(k), 'linear', -20);
else
g_az = interp1(pat.az_angles, pat.gain_az_dBi, Target.az(k), 'linear', -20);
g_el = interp1(pat.el_angles, pat.gain_el_dBi, Target.el(k), 'linear', -20);
g_dBi = g_az + g_el - pat.max_gain_dBi;
end
%
rx_gain_amp = sqrt(10^(g_dBi/10));
RxOut.space_loss_amp(k, rx, :, :, :) = RxOut.space_loss_amp(k, rx, :, :, :) * rx_gain_amp;
end
end
end
@@ -0,0 +1,118 @@
function [windowVector, metrics] = create_window_with_metrics(windowType, windowLength, varargin)
% CREATE_WINDOW_WITH_METRICS -
%
% :
% windowType - : 'none', 'hann', 'hamming', 'blackman', 'chebwin'
% windowLength - ( )
% varargin - (: chebwin의 sidelobe level)
%
% :
% windowVector - (1 x windowLength)
% metrics -
% .type :
% .length :
% .snr_loss_dB : SNR (dB)
% .scalloping_loss_dB: Scalloping (dB)
% .coherent_gain :
% .enbw : Equivalent Noise Bandwidth (bins)
%
if windowLength <= 0
error('Window length must be positive.');
end
%
switch lower(windowType)
case 'none'
windowVector = ones(1, windowLength);
case 'hann'
windowVector = hann(windowLength)';
case 'hamming'
windowVector = hamming(windowLength)';
case 'blackman'
windowVector = blackman(windowLength)';
case 'chebwin'
% Chebyshev sidelobe level (: 60dB)
if ~isempty(varargin)
sidelobe_dB = varargin{1};
else
sidelobe_dB = 60;
end
windowVector = chebwin(windowLength, sidelobe_dB)';
otherwise
warning('Unknown window type "%s". Using Hann window as default.', windowType);
windowVector = hann(windowLength)';
windowType = 'hann';
end
%
metrics = calculate_window_metrics(windowVector, windowType);
end
function metrics = calculate_window_metrics(windowVector, windowType)
% CALCULATE_WINDOW_METRICS -
%
% :
% 1. SNR Loss (dB) : SNR
% 2. Scalloping Loss (dB): FFT bin (0.5 bin offset)
% 3. Coherent Gain :
% 4. ENBW (bins) : Equivalent Noise Bandwidth
w = windowVector(:).'; % Row vector로
N = numel(w);
if N == 0
metrics = struct('type', windowType, 'length', 0, ...
'snr_loss_dB', NaN, 'scalloping_loss_dB', NaN, ...
'coherent_gain', NaN, 'enbw', NaN);
return;
end
% 1. Coherent Gain ( )
coherent_gain = mean(w);
% 2. Noise Power Gain ( )
noise_power_gain = mean(abs(w).^2);
% 3. SNR Loss (dB)
% SNR_loss = (Noise Power Gain) / (Coherent Gain)^2
%
if abs(coherent_gain) > eps
snr_loss_linear = noise_power_gain / (abs(coherent_gain)^2);
snr_loss_dB = 10 * log10(snr_loss_linear);
else
snr_loss_dB = NaN;
end
% 4. Equivalent Noise Bandwidth (ENBW)
% ENBW는 bin의
enbw = N * noise_power_gain / (sum(w)^2);
% 5. Scalloping Loss (dB)
% FFT bin (0.5 bin offset)
%
sample_index = 0:(N-1);
half_bin_response = abs(sum(w .* exp(-1j * 2 * pi * 0.5 * sample_index / N)));
dc_response = abs(sum(w));
if dc_response > eps
scalloping_loss_dB = -20 * log10(half_bin_response / dc_response);
else
scalloping_loss_dB = NaN;
end
%
metrics = struct(...
'type', lower(windowType), ...
'length', N, ...
'snr_loss_dB', snr_loss_dB, ...
'scalloping_loss_dB', scalloping_loss_dB, ...
'coherent_gain', coherent_gain, ...
'enbw', enbw);
end
+140
View File
@@ -0,0 +1,140 @@
function [det_mask, threshold_map, detections] = detect_targets_cfar(rd_map, RadarParams)
% CFAR target detection for range-doppler map
% - method: 'CA' or 'OS'
% - dimension: '1D' or '2D'
% - axis (for 1D): 'range' or 'doppler'
%
% Input:
% rd_map : [Ndoppler x Nrange] complex or real map
% RadarParams : struct containing SP.CFAR options
%
% Output:
% det_mask : logical detection mask [Ndoppler x Nrange]
% threshold_map : threshold map [Ndoppler x Nrange]
% detections : [Ndet x 2] = [doppler_bin, range_bin]
cfg = RadarParams.SP.CFAR;
method = upper(string(cfg.method));
dim_mode = upper(string(cfg.dimension));
axis_mode = lower(string(cfg.axis));
pfa = cfg.pfa;
train = cfg.train;
guard = cfg.guard;
if numel(train) == 1
train = [train, train];
end
if numel(guard) == 1
guard = [guard, guard];
end
os_rank_ratio = cfg.rank;
os_scale = cfg.os_scale;
rd_power = abs(rd_map).^2;
[n_dop, n_rng] = size(rd_power);
det_mask = false(n_dop, n_rng);
threshold_map = nan(n_dop, n_rng);
switch dim_mode
case "2D"
td = train(1); tr = train(2);
gd = guard(1); gr = guard(2);
for d = (td+gd+1):(n_dop-(td+gd))
for r = (tr+gr+1):(n_rng-(tr+gr))
d_idx = (d-(td+gd)):(d+(td+gd));
r_idx = (r-(tr+gr)):(r+(tr+gr));
win = rd_power(d_idx, r_idx);
cut_d = td+gd+1;
cut_r = tr+gr+1;
guard_mask = false(size(win));
guard_mask((cut_d-gd):(cut_d+gd), (cut_r-gr):(cut_r+gr)) = true;
train_cells = win(~guard_mask);
th = local_cfar_threshold(train_cells, method, pfa, os_rank_ratio, os_scale);
threshold_map(d, r) = th;
det_mask(d, r) = rd_power(d, r) > th;
end
end
case "1D"
switch axis_mode
case "range"
tr = train(2); gr = guard(2);
for d = 1:n_dop
[det_row, th_row] = cfar_1d_line(rd_power(d, :), tr, gr, method, pfa, os_rank_ratio, os_scale);
det_mask(d, :) = det_row;
threshold_map(d, :) = th_row;
end
case "doppler"
td = train(1); gd = guard(1);
for r = 1:n_rng
[det_col, th_col] = cfar_1d_line(rd_power(:, r).', td, gd, method, pfa, os_rank_ratio, os_scale);
det_mask(:, r) = det_col.';
threshold_map(:, r) = th_col.';
end
otherwise
error('CFAR axis must be ''range'' or ''doppler'' when dimension is 1D.');
end
otherwise
error('CFAR dimension must be ''1D'' or ''2D''.');
end
[d_idx, r_idx] = find(det_mask);
detections = [d_idx, r_idx];
end
function [det_line, th_line] = cfar_1d_line(x, t, g, method, pfa, rank_ratio, os_scale)
n = numel(x);
det_line = false(1, n);
th_line = nan(1, n);
left = t + g;
right = t + g;
for i = (left+1):(n-right)
l_train = x((i-g-t):(i-g-1));
r_train = x((i+g+1):(i+g+t));
train_cells = [l_train, r_train];
th = local_cfar_threshold(train_cells, method, pfa, rank_ratio, os_scale);
th_line(i) = th;
det_line(i) = x(i) > th;
end
end
function th = local_cfar_threshold(train_cells, method, pfa, rank_ratio, os_scale)
train_cells = train_cells(:);
n_train = numel(train_cells);
if n_train == 0
th = inf;
return;
end
switch method
case "CA"
noise_hat = mean(train_cells);
alpha = n_train * (pfa^(-1/n_train) - 1);
th = alpha * noise_hat;
case "OS"
sorted_cells = sort(train_cells, 'ascend');
k = max(1, min(n_train, round(rank_ratio * n_train)));
noise_hat = sorted_cells(k);
th = os_scale * noise_hat;
otherwise
error('CFAR method must be ''CA'' or ''OS''.');
end
end
+18
View File
@@ -0,0 +1,18 @@
function target_rd_map = integrate_nci_rdm(rd_cube, RadarParams)
% NCI (Noncoherent Integration) 2D RDM
% :
% rd_cube: [RX, TX, Doppler, Range]
% RadarParams.Waveform.mimoMode: 'TDM'
% :
% target_rd_map: [Doppler, Range]
%
% :
% - TDM : TX-RX (|.|^2)
% - : TX=1 RX
if strcmpi(RadarParams.Waveform.mimoMode, 'TDM')
target_rd_map = squeeze(sum(sum(abs(rd_cube).^2, 1), 2));
else
target_rd_map = squeeze(sum(abs(rd_cube(:, 1, :, :)).^2, 1));
end
end
@@ -0,0 +1,36 @@
function [rd_map, doppler_axis] = process_doppler_fft(range_profile, RadarParams)
% :
% - range_profile: [NumRx, NumTx, NumChirps, NumRangeBins]
% - RadarParams:
NumChirps = RadarParams.Waveform.NumChirps;
[~, ~, ~, ~] = size(range_profile);
%
lambda = RadarParams.Waveform.lambda_c; % = c/fc
% 1. (PRI, Pulse Repetition Interval)
T_pri = RadarParams.Waveform.PRI;
% 2. Main.m에서 Doppler-FFT용
win_doppler = RadarParams.SP.RDM.window_doppler;
if length(win_doppler) ~= NumChirps
error('Doppler window length mismatch: expected %d, got %d', NumChirps, length(win_doppler));
end
win_data = range_profile .* reshape(win_doppler, [1, 1, NumChirps, 1]);
% 3. Doppler-FFT (3 : Chirp)
% NFFT를 NumChirps보다 .
rd_fft = fft(win_data, NumChirps, 3);
% 4. fftshift ( 0 )
% [-, 0, +] .
rd_map = fftshift(rd_fft, 3);
% 5. (Velocity Axis)
% Vmax = lambda / (4 * T_pri)
v_max = lambda / (4 * T_pri);
% dv = lambda / (2 * NumChirps * T_pri)
doppler_axis = linspace(-v_max, v_max, NumChirps);
end
@@ -0,0 +1,40 @@
function [range_profile, range_axis] = process_range_fft_lpf(adc_raw_data, RadarParams)
% :
% - adc_raw_data: [NumRx, NumTx, NumChirps, N_samples]
% - RadarParams:
fs_adc = RadarParams.Waveform.fs_adc;
Slope = RadarParams.Waveform.Slope;
fc_lpf_Hz = RadarParams.Rxpath.fc_lpf;
[~, ~, ~, N_samples] = size(adc_raw_data);
c = RadarParams.Basic.c;
% 1. Main.m에서
win = RadarParams.SP.RDM.window_range;
if length(win) ~= N_samples
error('Window length mismatch: expected %d, got %d', N_samples, length(win));
end
win_data = adc_raw_data .* reshape(win, [1, 1, 1, N_samples]);
% 2. Range-FFT
range_fft = fft(win_data, N_samples, 4);
% 3.
df = fs_adc / N_samples;
freq_axis = (0 : N_samples-1) * df;
full_range_axis = (freq_axis * c) / (2 * Slope);
% 4. Ideal LPF
cutoff_idx = floor(fc_lpf_Hz / df);
%
half_idx = floor(N_samples/2);
final_idx = min(cutoff_idx, half_idx);
% 5. [ ] (Truncation)
% 1 final_idx까지만
range_profile = range_fft(:, :, :, 1:final_idx);
range_axis = full_range_axis(1:final_idx);
end
@@ -0,0 +1,175 @@
function fig = visualize_rd_map_with_spurs(target_rd_map, r_axis, v_axis, Target, SpurParams, RadarParams, Slope)
% 2D Range-Doppler .
%
% :
% - target_rd_map: 2D RD [Doppler x Range]
% - r_axis:
% - v_axis:
% - Target: (R, v, NumTargets )
% - SpurParams: (lo, adc, mixer, pulse, lo_leak, clip )
% - RadarParams:
% - Slope: (Hz/s)
%
% :
% - fig: figure
if ~isfield(RadarParams, 'SP') || ~isfield(RadarParams.SP, 'RDM') || ...
~isfield(RadarParams.SP.RDM, 'window_range') || ~isfield(RadarParams.SP.RDM, 'window_doppler')
error('RadarParams.SP.RDM.window_range / window_doppler .');
end
win_r = RadarParams.SP.RDM.window_range;
win_d = RadarParams.SP.RDM.window_doppler;
normal_factor = RadarParams.Rxpath.fc_lpf * RadarParams.Waveform.PRF * sum(win_r.^2) * sum(win_d.^2);
if ~isfield(RadarParams, 'Rxpath') || ~isfield(RadarParams.Rxpath, 'load_impedance_ohm')
error('RadarParams.Rxpath.load_impedance_ohm .');
end
load_impedance_ohm = RadarParams.Rxpath.load_impedance_ohm;
psd2d_w_hz2 = (abs(target_rd_map).^2) / (load_impedance_ohm * max(normal_factor, eps));
rbw_range_hz = RadarParams.Waveform.fs_adc * sum(win_r.^2) / max(eps, (sum(win_r)^2));
rbw_doppler_hz = RadarParams.Waveform.PRF * sum(win_d.^2) / max(eps, (sum(win_d)^2));
rd_map_dbm_bin = 10*log10(max(psd2d_w_hz2, realmin) / 1e-3) + 10*log10(max(rbw_range_hz * rbw_doppler_hz, eps));
fig = figure('Name', '2D Range-Doppler Map');
imagesc(r_axis, v_axis, rd_map_dbm_bin);
axis xy; % y축 ()
cb = colorbar;
cb.Label.String = 'Power (dBm/bin)';
xlabel('Range (m)');
ylabel('Velocity (m/s)');
title('Range-Doppler Map (dBm/bin)');
colormap(jet);
%
hLines = gobjects(0);
hLabels = {};
% --- 1. ---
hold on;
hTarget = gobjects(0);
if exist('Target','var') && isstruct(Target)
for k = 1:Target.NumTargets
% /
if Target.R(k) < min(r_axis) || Target.R(k) > max(r_axis) || ...
Target.v(k) < min(v_axis) || Target.v(k) > max(v_axis)
continue; %
end
[~, ix] = min(abs(r_axis - Target.R(k)));
[~, iy] = min(abs(v_axis - Target.v(k)));
hTarget(end+1) = plot(r_axis(ix), v_axis(iy), 'ro', 'MarkerSize',8, 'LineWidth',1.5);
text(r_axis(ix), v_axis(iy), sprintf(' T%d', k), 'Color','r','FontSize',8);
end
end
hold off;
%
if exist('hTarget','var') && ~isempty(hTarget)
hLines(end+1) = hTarget(1);
hLabels{end+1} = 'actual target';
end
% --- nonideal/spur ---
show_nonideal_overlay = true;
if exist('SpurParams','var') && isstruct(SpurParams) && isfield(SpurParams,'enabled')
show_nonideal_overlay = logical(SpurParams.enabled);
end
if exist('RadarParams','var') && isstruct(RadarParams) && ...
isfield(RadarParams,'Waveform') && isfield(RadarParams.Waveform,'nonideal')
ni = RadarParams.Waveform.nonideal;
enable_phase_noise = false;
enable_nonlinearity = false;
if isfield(ni,'enable_phase_noise')
enable_phase_noise = logical(ni.enable_phase_noise);
end
if isfield(ni,'enable_nonlinearity')
enable_nonlinearity = logical(ni.enable_nonlinearity);
end
show_nonideal_overlay = show_nonideal_overlay && ...
(enable_phase_noise || enable_nonlinearity);
end
% --- 2. LO/ADC ---
if show_nonideal_overlay && exist('SpurParams','var') && isstruct(SpurParams)
spur_freqs = [];
spur_names = {};
if isfield(SpurParams,'lo') && ~isempty(SpurParams.lo)
spur_freqs(end+1) = SpurParams.lo.freq;
spur_names{end+1} = 'LO spur';
end
if isfield(SpurParams,'adc') && ~isempty(SpurParams.adc)
spur_freqs(end+1) = SpurParams.adc.freq;
spur_names{end+1} = 'ADC spur';
end
hold on;
for idx = 1:length(spur_freqs)
f = spur_freqs(idx);
name = sprintf('%s @ %.2f MHz', spur_names{idx}, f/1e6);
r_spur = (RadarParams.Basic.c * f) / (2 * Slope);
hLines(end+1) = xline(r_spur, '--m', name);
hLabels{end+1} = name;
end
% --- 3. Pulse comb ---
if isfield(SpurParams,'pulse') && ~isempty(SpurParams.pulse)
p = SpurParams.pulse;
num_harm = 5;
for n = 1:num_harm
f_n = n * p.freq;
r_n = (RadarParams.Basic.c * f_n) / (2 * Slope);
if r_n <= max(r_axis)
lbl = sprintf('pulse comb %d', n);
hLines(end+1) = xline(r_n, '--g', lbl);
hLabels{end+1} = lbl;
end
end
end
% --- 4. Mixer ---
if isfield(SpurParams,'mixer') && ~isempty(SpurParams.mixer) && exist('Target','var')
base_ranges = Target.R;
harmonics = [2,3];
colors = {'--c','--y'};
for hi = 1:length(harmonics)
mul = harmonics(hi);
for r0 = base_ranges
r_spur2 = r0 * mul;
if r_spur2 <= max(r_axis)
lbl = sprintf('mixer x%d', mul);
hLines(end+1) = xline(r_spur2, colors{hi}, lbl);
hLabels{end+1} = lbl;
end
end
end
end
% --- 5. DC offset / LO leakage ---
if isfield(SpurParams,'lo_leak') && ~isempty(SpurParams.lo_leak)
hLines(end+1) = line([min(r_axis), max(r_axis)], [0, 0], 'Color','r','LineStyle',':', 'LineWidth',1.5);
hLabels{end+1} = 'DC offset/LO leak';
end
% --- 6. Clipping spur ---
if isfield(SpurParams,'clip') && ~isempty(SpurParams.clip)
r_clip = SpurParams.clip.range;
if r_clip >= min(r_axis) && r_clip <= max(r_axis)
hLines(end+1) = line([r_clip, r_clip], [min(v_axis), max(v_axis)], 'Color', 'r', 'LineStyle', ':', 'LineWidth', 1.2);
hLabels{end+1} = sprintf('clipping @ %.1f m', r_clip);
end
end
hold off;
end
% --- 7. ---
if ~isempty(hLines)
legend(hLines, hLabels, 'Location', 'northeastoutside');
end
end
+379
View File
@@ -0,0 +1,379 @@
function [snr_coverage_2d, azimuth_deg, elevation_deg, coverage_info, fig] = analyze_coverage(RadarParams, target_range_m, target_rcs_dBsm)
% Coverage : 2D SNR
%
% :
% 1) Azimuth 0, Elevation 0 SNR
% 2) SNR ()
% 3) RCS
%
% :
% - RadarParams:
% - target_range_m: (m), 100
% - target_rcs_dBsm: RCS (dBsm), 0
%
% :
% - snr_coverage_2d: 2D SNR [elevation x azimuth]
% - azimuth_deg:
% - elevation_deg:
% - coverage_info:
% - fig: figure
%%
if nargin < 2
target_range_m = 100; % : 100m
end
if nargin < 3
target_rcs_dBsm = 0; % : 0 dBsm = 1
end
%%
lambda = RadarParams.Waveform.lambda_c;
fc = RadarParams.Waveform.fc;
kb = RadarParams.Basic.kb;
T0 = RadarParams.Basic.T0;
%
TxPattern = RadarParams.Antenna.TxPattern;
RxPattern = RadarParams.Antenna.RxPattern;
NumTx = RadarParams.Antenna.NumTx;
NumRx = RadarParams.Antenna.NumRx;
%
rxPathGain_dB = RadarParams.Rxpath.rxPathGain_dB;
system_NF_dB = RadarParams.Rxpath.system_NF_dB;
% RCS ()
rcs_target = 10^(target_rcs_dBsm / 10);
%% 2D (1D -> 2D )
% ( )
az_common = -90:1:90;
el_common = -90:1:90;
[AZ_common_grid, EL_common_grid] = meshgrid(az_common, el_common);
% TX 2D로
avg_tx_gain_2d = zeros(length(el_common), length(az_common));
for tx = 1:NumTx
if isfield(TxPattern(tx), 'gain_az_dBi') % 1D
az_angles = TxPattern(tx).az_angles;
el_angles = TxPattern(tx).el_angles;
gain_az = TxPattern(tx).gain_az_dBi;
gain_el = TxPattern(tx).gain_el_dBi;
%
max_gain_az = max(gain_az);
max_gain_el = max(gain_el);
max_gain_ref = max(max_gain_az, max_gain_el);
% Normalize ( = 0 dB)
gain_az_norm = gain_az - max_gain_az;
gain_el_norm = gain_el - max_gain_el;
%
gain_az_interp = interp1(az_angles, gain_az_norm, AZ_common_grid, 'linear', 'extrap');
gain_el_interp = interp1(el_angles, gain_el_norm, EL_common_grid, 'linear', 'extrap');
% 2D로
tx_gain_2d = gain_az_interp + gain_el_interp + max_gain_ref;
else % 2D
az_angles = TxPattern(tx).az_angles;
el_angles = TxPattern(tx).el_angles;
gain_dBi = TxPattern(tx).gain_dBi;
[AZ_grid, EL_grid] = ndgrid(az_angles, el_angles);
F = scatteredInterpolant(AZ_grid(:), EL_grid(:), gain_dBi(:), 'linear', 'nearest');
tx_gain_2d = F(AZ_common_grid, EL_common_grid);
end
avg_tx_gain_2d = avg_tx_gain_2d + tx_gain_2d;
end
avg_tx_gain_2d = avg_tx_gain_2d / NumTx;
% RX 2D로 ( )
avg_rx_gain_2d = zeros(length(el_common), length(az_common));
for rx = 1:NumRx
if isfield(RxPattern(rx), 'gain_az_dBi') % 1D
az_angles = RxPattern(rx).az_angles;
el_angles = RxPattern(rx).el_angles;
gain_az = RxPattern(rx).gain_az_dBi;
gain_el = RxPattern(rx).gain_el_dBi;
%
max_gain_az = max(gain_az);
max_gain_el = max(gain_el);
max_gain_ref = max(max_gain_az, max_gain_el);
% Normalize ( = 0 dB)
gain_az_norm = gain_az - max_gain_az;
gain_el_norm = gain_el - max_gain_el;
%
gain_az_interp = interp1(az_angles, gain_az_norm, AZ_common_grid, 'linear', 'extrap');
gain_el_interp = interp1(el_angles, gain_el_norm, EL_common_grid, 'linear', 'extrap');
% 2D로
rx_gain_2d = gain_az_interp + gain_el_interp + max_gain_ref;
else % 2D
az_angles = RxPattern(rx).az_angles;
el_angles = RxPattern(rx).el_angles;
gain_dBi = RxPattern(rx).gain_dBi;
[AZ_grid, EL_grid] = ndgrid(az_angles, el_angles);
F = scatteredInterpolant(AZ_grid(:), EL_grid(:), gain_dBi(:), 'linear', 'nearest');
rx_gain_2d = F(AZ_common_grid, EL_common_grid);
end
avg_rx_gain_2d = avg_rx_gain_2d + rx_gain_2d;
end
avg_rx_gain_2d = avg_rx_gain_2d / NumRx;
%% TX
% TX
pa_profile_freq = RadarParams.RFOutput.PA_Profile.freqs;
pa_profile_power_dbm = RadarParams.RFOutput.PA_Profile.power_dBm;
ptx_dbm = interp1(pa_profile_freq, pa_profile_power_dbm, fc, 'linear', 'extrap');
ptx_w = 10^((ptx_dbm - 30) / 10);
% ( )
path_loss_factor = (4 * pi * target_range_m)^2;
rxGain_linear = 10^(rxPathGain_dB / 10);
noise_power_w = kb * T0 * RadarParams.Waveform.fs_adc;
system_NF_linear = 10^(system_NF_dB / 10);
total_noise_power_w = noise_power_w * system_NF_linear * rxGain_linear; % RX
%% SNR (2D - )
% X축: Azimuth, Y축: Elevation, Z축: SNR
azimuth_deg = -90:1:90; % X축: -90 ~ 90, 1
elevation_deg = -90:1:90; % Y축: -90 ~ 90, 1
% az_common과 elevation_deg가
% (dB)
g_tx_db_grid = avg_tx_gain_2d; % [num_el x num_az]
g_rx_db_grid = avg_rx_gain_2d; % [num_el x num_az]
%
g_tx_linear_grid = 10.^(g_tx_db_grid / 10);
g_rx_linear_grid = 10.^(g_rx_db_grid / 10);
%
SP_gain_rngFFT = RadarParams.Waveform.Timing.AdcSampTime * RadarParams.Waveform.fs_adc;
SP_gain_dopFFT = RadarParams.Waveform.NumChirps;
SP_loss_rngwin = 10^(RadarParams.SP.RDM.window_metrics_range.snr_loss_dB / 10);
SP_loss_dopwin = 10^(RadarParams.SP.RDM.window_metrics_doppler.snr_loss_dB / 10);
SP_loss_rng_straddle = 10^(RadarParams.SP.RDM.window_metrics_range.scalloping_loss_dB / 10);
SP_loss_dop_straddle = 10^(RadarParams.SP.RDM.window_metrics_doppler.scalloping_loss_dB / 10);
SP_total = SP_gain_rngFFT * SP_gain_dopFFT / (SP_loss_rngwin * SP_loss_dopwin * SP_loss_rng_straddle * SP_loss_dop_straddle);
% Secondary Surface Loss
secondary_loss = 10^(RadarParams.Antenna.SecondarySurfaceLoss_dB / 10);
fprintf('SP gain rng FFT: %.2f, SP gain doppler FFT: %.2f\n', SP_gain_rngFFT, SP_gain_dopFFT);
fprintf('SP loss range window: %.2f dB, SP loss doppler window: %.2f dB\n', RadarParams.SP.RDM.window_metrics_range.snr_loss_dB, RadarParams.SP.RDM.window_metrics_doppler.snr_loss_dB);
fprintf('SP loss range straddle: %.2f dB, SP loss doppler straddle: %.2f dB\n', RadarParams.SP.RDM.window_metrics_range.scalloping_loss_dB, RadarParams.SP.RDM.window_metrics_doppler.scalloping_loss_dB);
fprintf('SP loss doppler window: %.2f dB\n', RadarParams.SP.RDM.window_metrics_doppler.snr_loss_dB);
fprintf('SP loss doppler straddle: %.2f dB\n', RadarParams.SP.RDM.window_metrics_doppler.scalloping_loss_dB);
fprintf('Secondary surface loss: %.2f dB\n', RadarParams.Antenna.SecondarySurfaceLoss_dB);
% SNR : ( × 2D )
p_rx_w = (ptx_w * (lambda^2) * rcs_target) / ((path_loss_factor^2) * (4 * pi)) * SP_total / secondary_loss;
p_rx_w = p_rx_w * g_tx_linear_grid .* g_rx_linear_grid; %
p_rx_after_rxgain_w = p_rx_w * rxGain_linear;
snr_linear = p_rx_after_rxgain_w / total_noise_power_w;
snr_coverage_2d = 10 * log10(max(snr_linear, eps));
%% 1D CUT (elevation=0 Azimuth cut, azimuth=0 Elevation cut)
[~, el_idx_zero] = min(abs(elevation_deg - 0));
snr_az_cut = snr_coverage_2d(el_idx_zero, :); % elevation=0
[~, az_idx_zero] = min(abs(azimuth_deg - 0));
snr_el_cut = snr_coverage_2d(:, az_idx_zero); % azimuth=0
%% Cartesian (Range, Azimuth, Elevation -> x, y, z)
[AZ_deg_grid, EL_deg_grid] = meshgrid(azimuth_deg, elevation_deg);
x_grid_m = target_range_m .* cosd(EL_deg_grid) .* cosd(AZ_deg_grid);
y_grid_m = target_range_m .* cosd(EL_deg_grid) .* sind(AZ_deg_grid);
z_grid_m = target_range_m .* sind(EL_deg_grid);
%% Coverage
snr_2d_vec = snr_coverage_2d(:);
coverage_info.azimuth_deg = azimuth_deg;
coverage_info.elevation_deg = elevation_deg;
coverage_info.snr_coverage_2d = snr_coverage_2d;
coverage_info.snr_az_cut = snr_az_cut;
coverage_info.snr_el_cut = snr_el_cut;
coverage_info.target_range_m = target_range_m;
coverage_info.target_rcs_dBsm = target_rcs_dBsm;
coverage_info.x_grid_m = x_grid_m;
coverage_info.y_grid_m = y_grid_m;
coverage_info.z_grid_m = z_grid_m;
coverage_info.mean_snr = mean(snr_2d_vec);
coverage_info.max_snr = max(snr_2d_vec);
coverage_info.min_snr = min(snr_2d_vec);
[~, max_idx_2d] = max(snr_2d_vec);
[max_el_idx_2d, max_az_idx_2d] = ind2sub(size(snr_coverage_2d), max_idx_2d);
coverage_info.max_snr_azimuth = azimuth_deg(max_az_idx_2d);
coverage_info.max_snr_elevation = elevation_deg(max_el_idx_2d);
coverage_info.max_snr_x_m = x_grid_m(max_el_idx_2d, max_az_idx_2d);
coverage_info.max_snr_y_m = y_grid_m(max_el_idx_2d, max_az_idx_2d);
coverage_info.max_snr_z_m = z_grid_m(max_el_idx_2d, max_az_idx_2d);
%% 3D Coverage
fig = figure('Name', 'Coverage Analysis - 3D');
set(fig, 'Position', [100, 100, 1600, 900]);
% 1D cut SNR
[~, el_idx_zero] = min(abs(elevation_deg - 0));
[~, az_idx_zero] = min(abs(azimuth_deg - 0));
snr_az_cut = snr_coverage_2d(el_idx_zero, :); % elevation=0 cut
snr_el_cut = snr_coverage_2d(:, az_idx_zero); % azimuth=0 cut
% Subplot 1: x-y top view (color = SNR at elevation=0 cut)
ax_topview = subplot(2, 2, 3);
x_top_m = target_range_m .* cosd(azimuth_deg);
y_top_m = target_range_m .* sind(azimuth_deg);
scatter(x_top_m, y_top_m, 45, snr_az_cut, 'filled');
axis equal;
grid on;
cb = colorbar;
ylabel(cb, 'SNR (dB)');
xlabel('x (m)');
ylabel('y (m)');
title('Top View (Azimuth)');
% Subplot 2: 3D Surface (geometry = x,y,z, color = SNR) - spans full width at top
ax_3d = subplot(2, 2, 1:2);
surf(x_grid_m, y_grid_m, z_grid_m, snr_coverage_2d, 'EdgeColor', 'none');
colormap(jet);
cb_3d = colorbar;
ylabel(cb_3d, 'SNR (dB)');
xlabel('x (m)');
ylabel('y (m)');
zlabel('z (m)');
title('3D Coverage in Cartesian Coordinates');
view(45, 30);
grid on;
hold on;
plot3(coverage_info.max_snr_x_m, coverage_info.max_snr_y_m, coverage_info.max_snr_z_m, 'r*', 'MarkerSize', 20, 'LineWidth', 2);
hold off;
% Subplot 3: x-z side view (azimuth=0, color=SNR at azimuth=0 cut)
ax_sideview = subplot(2, 2, 4);
x_side_m = target_range_m .* cosd(elevation_deg);
z_side_m = target_range_m .* sind(elevation_deg);
scatter(x_side_m, z_side_m, 45, snr_el_cut, 'filled');
grid on;
cb_side = colorbar;
ylabel(cb_side, 'SNR (dB)');
xlabel('x (m)');
ylabel('z (m)');
title('Side View (Elevation)');
%
sgtitle(sprintf('3D Coverage (x,y,z) | Max SNR: %.2f dB at (%.2f, %.2f, %.2f) m', ...
coverage_info.max_snr, coverage_info.max_snr_x_m, coverage_info.max_snr_y_m, coverage_info.max_snr_z_m), 'FontSize', 12);
% Data Cursor Mode ( )
dcm_obj = datacursormode(fig);
dcm_obj.Enable = 'on';
set(dcm_obj, 'UpdateFcn', {@unified_datatip_callback, ...
ax_topview, ax_3d, ax_sideview, ...
x_top_m, y_top_m, snr_az_cut, azimuth_deg, ...
x_grid_m, y_grid_m, z_grid_m, snr_coverage_2d, elevation_deg, ...
x_side_m, z_side_m, snr_el_cut});
dcm_obj.SnapToDataVertex = 'on';
% Help text
fprintf('\n>>> Data Cursor Mode ENABLED <<<\n');
fprintf('Instructions:\n');
fprintf(' 1. Left-click on any plot to inspect points\n');
fprintf(' 2. Top View shows: x, y, Azimuth, SNR\n');
fprintf(' 3. 3D Surface shows: x, y, z, Azimuth, Elevation, SNR\n');
fprintf(' 4. Side View shows: x, z, Elevation, SNR\n');
fprintf(' 5. Press Escape or click "Disable Data Cursor" in Figure Tools to deactivate\n\n');
%%
fprintf('\n===== 3D COVERAGE ANALYSIS REPORT (Cartesian x,y,z) =====\n');
fprintf('Target Configuration:\n');
fprintf(' - Range: %.1f m\n', target_range_m);
fprintf(' - RCS: %.1f dBsm\n', target_rcs_dBsm);
fprintf(' - Analysis Angles: Azimuth ±90 deg, Elevation ±90 deg\n');
fprintf('\nSNR Statistics (Full Grid):\n');
fprintf(' - Mean SNR: %.2f dB\n', coverage_info.mean_snr);
fprintf(' - Max SNR: %.2f dB @ (x,y,z)=(%.2f, %.2f, %.2f) m\n', coverage_info.max_snr, coverage_info.max_snr_x_m, coverage_info.max_snr_y_m, coverage_info.max_snr_z_m);
fprintf(' - Min SNR: %.2f dB\n', coverage_info.min_snr);
fprintf(' - SNR Range: %.2f dB\n', coverage_info.max_snr - coverage_info.min_snr);
fprintf('System Parameters:\n');
fprintf(' - TX Power: %.2f dBm\n', ptx_dbm);
fprintf(' - RX Gain: %.1f dB\n', rxPathGain_dB);
fprintf(' - System NF: %.1f dB\n', system_NF_dB);
fprintf('=============================================\n\n');
end
function output_txt = unified_datatip_callback(~, event, ...
ax_topview, ax_3d, ax_sideview, ...
x_top_m, y_top_m, snr_az_cut, azimuth_deg, ...
x_grid_m, y_grid_m, z_grid_m, snr_coverage_2d, elevation_deg, ...
x_side_m, z_side_m, snr_el_cut)
% Unified datatip callback for all subplots
try
pos = event.Position;
current_ax = event.Target.Parent; % Get current axis
% Determine which subplot was clicked
if isequal(current_ax, ax_3d)
% 3D Plot
x_val = pos(1);
y_val = pos(2);
z_val = pos(3);
dist2 = (x_grid_m - x_val).^2 + (y_grid_m - y_val).^2 + (z_grid_m - z_val).^2;
[~, nearest_idx] = min(dist2(:));
[nearest_el_idx, nearest_az_idx] = ind2sub(size(dist2), nearest_idx);
snr_val = snr_coverage_2d(nearest_el_idx, nearest_az_idx);
az_val = azimuth_deg(nearest_az_idx);
el_val = elevation_deg(nearest_el_idx);
output_txt = {
['x: ' num2str(x_val, '%.3f') ' m']
['y: ' num2str(y_val, '%.3f') ' m']
['z: ' num2str(z_val, '%.3f') ' m']
['Azimuth: ' num2str(az_val, '%.2f') ' deg']
['Elevation: ' num2str(el_val, '%.2f') ' deg']
['SNR[dB]: ' num2str(snr_val, '%.3f')]
};
elseif isequal(current_ax, ax_topview)
% Top View - x, y, Azimuth, SNR
x_val = pos(1);
y_val = pos(2);
[~, nearest_az_idx] = min(abs(x_top_m - x_val).^2 + abs(y_top_m - y_val).^2);
az_val = azimuth_deg(nearest_az_idx);
snr_val = snr_az_cut(nearest_az_idx);
output_txt = {
['x: ' num2str(x_val, '%.3f') ' m']
['y: ' num2str(y_val, '%.3f') ' m']
['Azimuth: ' num2str(az_val, '%.2f') ' deg']
['SNR[dB]: ' num2str(snr_val, '%.3f')]
};
elseif isequal(current_ax, ax_sideview)
% Side View - x, z, Elevation, SNR
x_val = pos(1);
z_val = pos(2);
[~, nearest_el_idx] = min(abs(x_side_m - x_val).^2 + abs(z_side_m - z_val).^2);
el_val = elevation_deg(nearest_el_idx);
snr_val = snr_el_cut(nearest_el_idx);
output_txt = {
['x: ' num2str(x_val, '%.3f') ' m']
['z: ' num2str(z_val, '%.3f') ' m']
['Elevation: ' num2str(el_val, '%.2f') ' deg']
['SNR[dB]: ' num2str(snr_val, '%.3f')]
};
else
output_txt = 'Unknown plot';
end
catch ME
output_txt = ['Error: ' ME.message];
end
end
@@ -0,0 +1,68 @@
function PatternArray = build_antenna_patterns(NumAnt, PatternType, MaxGain, AzSquintStep, AzBeamwidth, ElBeamwidth, FileList)
% build_antenna_patterns: .
%
% []
% - FileList (): Cell Array.
% : {'tx1.mat', 'tx2.mat', ''} ( )
% ( )
az_angles_default = -90:1:90;
el_angles_default = -90:1:90;
% nargin을 FileList가
if nargin < 7
FileList = {};
end
for n = 1:NumAnt
PatternArray(n).Type = PatternType;
% -------------------------------------------------------------
% Case A: , (Load)
% -------------------------------------------------------------
if length(FileList) >= n && ~isempty(FileList{n}) && isfile(FileList{n})
filePath = FileList{n};
loadedData = load(filePath); % .mat
% (az_angles, el_angles ) .
PatternArray(n).az_angles = loadedData.az_angles;
PatternArray(n).el_angles = loadedData.el_angles;
if strcmpi(PatternType, '2D')
PatternArray(n).gain_dBi = loadedData.gain_dBi;
elseif strcmpi(PatternType, '1D')
PatternArray(n).max_gain_dBi = loadedData.max_gain_dBi;
PatternArray(n).gain_az_dBi = loadedData.gain_az_dBi;
PatternArray(n).gain_el_dBi = loadedData.gain_el_dBi;
end
%
continue; %
end
% -------------------------------------------------------------
% Case B: ( )
% -------------------------------------------------------------
PatternArray(n).az_angles = az_angles_default;
PatternArray(n).el_angles = el_angles_default;
% (Squint Angle)
% (3): (1-2)*5 = -5, (2-2)*5 = 0, (3-2)*5 = +5
squint_az = (n - (NumAnt + 1)/2) * AzSquintStep;
if strcmpi(PatternType, '2D')
[AZ, EL] = meshgrid(PatternArray(n).az_angles, PatternArray(n).el_angles);
PatternArray(n).gain_dBi = MaxGain - 3*((AZ - squint_az)/AzBeamwidth).^2 - 3*(EL/ElBeamwidth).^2;
PatternArray(n).gain_dBi(PatternArray(n).gain_dBi < -20) = -20;
elseif strcmpi(PatternType, '1D')
PatternArray(n).max_gain_dBi = MaxGain;
PatternArray(n).gain_az_dBi = MaxGain - 3*((PatternArray(n).az_angles - squint_az)/AzBeamwidth).^2;
PatternArray(n).gain_az_dBi(PatternArray(n).gain_az_dBi < -20) = -20;
PatternArray(n).gain_el_dBi = MaxGain - 3*(PatternArray(n).el_angles/ElBeamwidth).^2;
PatternArray(n).gain_el_dBi(PatternArray(n).gain_el_dBi < -20) = -20;
end
end
end
@@ -0,0 +1,81 @@
function fig = visualize_antenna_pattern(PatternArray, AntennaName)
% visualize_antenna_pattern: .
% AntennaName (: 'TX', 'RX') .
if nargin < 2
AntennaName = 'Antenna'; %
end
fig_name = sprintf('%s Radiation Pattern Cuts', AntennaName);
fig = figure('Name', fig_name, 'Position', [200, 200, 1000, 700]);
NumAnt = length(PatternArray);
colors = lines(NumAnt);
% [ Y축 ]
y_min_total = inf;
y_max_total = -inf;
for i = 1:NumAnt
pat = PatternArray(i);
if strcmpi(pat.Type, '2D')
y_min_total = min(y_min_total, min(pat.gain_dBi(:)));
y_max_total = max(y_max_total, max(pat.gain_dBi(:)));
else
y_min_total = min(y_min_total, min([pat.gain_az_dBi(:); pat.gain_el_dBi(:)]));
y_max_total = max(y_max_total, max([pat.gain_az_dBi(:); pat.gain_el_dBi(:)]));
end
end
% ---------------------------------------------------------
% [Subplot 1] : Azimuth Pattern Cut
% ---------------------------------------------------------
subplot(2, 1, 1);
hold on; grid on;
for i = 1:NumAnt
pat = PatternArray(i);
if strcmpi(pat.Type, '2D')
el0_idx = find(abs(pat.el_angles - 0) < 1e-6, 1);
az_gain_cut = pat.gain_dBi(el0_idx, :);
else
az_gain_cut = pat.gain_az_dBi;
end
plot(pat.az_angles, az_gain_cut, '-', 'Color', colors(i,:), 'LineWidth', 2.5, 'DisplayName', sprintf('%s %d', AntennaName, i));
end
title(sprintf('[%s] Azimuth Radiation Pattern Cut (at Elevation = 0^\\circ)', AntennaName), 'FontSize', 12);
xlabel('Azimuth Angle (deg)', 'FontSize', 11);
ylabel('Antenna Gain (dBi)', 'FontSize', 11);
xlim([min(PatternArray(1).az_angles), max(PatternArray(1).az_angles)]);
ylim([y_min_total - 2, y_max_total + 2]);
legend('show', 'Location', 'south');
hold off;
% ---------------------------------------------------------
% [Subplot 2] : Elevation Pattern Cut
% ---------------------------------------------------------
subplot(2, 1, 2);
hold on; grid on;
for i = 1:NumAnt
pat = PatternArray(i);
if strcmpi(pat.Type, '2D')
az0_idx = find(abs(pat.az_angles - 0) < 1e-6, 1);
el_gain_cut = pat.gain_dBi(:, az0_idx);
else
el_gain_cut = pat.gain_el_dBi;
end
plot(pat.el_angles, el_gain_cut, '-', 'Color', colors(i,:), 'LineWidth', 2.5, 'DisplayName', sprintf('%s %d', AntennaName, i));
end
title(sprintf('[%s] Elevation Radiation Pattern Cut (at Azimuth = 0^\\circ)', AntennaName), 'FontSize', 12);
xlabel('Elevation Angle (deg)', 'FontSize', 11);
ylabel('Antenna Gain (dBi)', 'FontSize', 11);
xlim([min(PatternArray(1).el_angles), max(PatternArray(1).el_angles)]);
ylim([y_min_total - 2, y_max_total + 2]);
legend('show', 'Location', 'south');
hold off;
end
+307
View File
@@ -0,0 +1,307 @@
% =========================================================================
% Tx_Main.m
% Pure MATLAB FMCW TX Simulator (Multi-Chirp Frame Generation)
% =========================================================================
clear; clc; close all;
% =============== ===============
currentFilePath = mfilename('fullpath');
currentFolder = fileparts(currentFilePath);
addpath(genpath(currentFolder));
%% =================== ===================
% figure를 on/off .
PlotToggle.range_profile = false; % Step 8 Range Profile
PlotToggle.tx_single = false; % Figure 1 Single chirp waveform
PlotToggle.tx_multi = false; % Figure 2 Multi chirp waveform
PlotToggle.tx_antenna = false; % Figure 3 TX antenna pattern
PlotToggle.rx_antenna = false; % Figure 4 RX antenna pattern
PlotToggle.rd_map = false; % Figure 5 Range-Doppler map
PlotToggle.cfar = false; % Figure 6 CFAR detections
PlotToggle.coverage = true; % Figure 7 Coverage analysis
%% =================== ===================
% 1)
RadarParams.Basic.c = 3e8; % (m/s)
RadarParams.Basic.kb = physconst('Boltzmann'); %
RadarParams.Basic.T0 = 290; % (K)
% 2)
RadarParams.Waveform.fc = 77e9; % (Hz)
RadarParams.Waveform.lambda_c = RadarParams.Basic.c / RadarParams.Waveform.fc; % (m)
RadarParams.Waveform.NumChirps = 128; %
RadarParams.Waveform.B_valid = 1e9; % (Hz)
RadarParams.Waveform.fs_adc = 10e6; % ADC (Hz)
RadarParams.Waveform.fs_waveform = 1e6; % (Hz)
% 2-1)
RadarParams.Waveform.Timing.IdleTime = 7e-6; % (s)
RadarParams.Waveform.Timing.TxStartTime = 1e-6; % Tx on (s)
RadarParams.Waveform.Timing.AdcStartTime = 6e-6; % ADC on (s)
RadarParams.Waveform.Timing.AdcSampTime = 50e-6; % ADC (s)
RadarParams.Waveform.Timing.ExcessTime = 1e-6; % (s)
RadarParams.Waveform.Timing.RampEndTime = RadarParams.Waveform.Timing.IdleTime + RadarParams.Waveform.Timing.TxStartTime + RadarParams.Waveform.Timing.AdcSampTime + RadarParams.Waveform.Timing.ExcessTime; % (s)
RadarParams.Waveform.PRI = RadarParams.Waveform.Timing.IdleTime + RadarParams.Waveform.Timing.RampEndTime; % Pulse Repetition Interval (s)
RadarParams.Waveform.PRF = 1 / RadarParams.Waveform.PRI; % Pulse Repetition Frequency (Hz)
RadarParams.Waveform.Slope = RadarParams.Waveform.B_valid / RadarParams.Waveform.Timing.AdcSampTime; % (Hz/s)
RadarParams.Waveform.f_start = RadarParams.Waveform.fc - RadarParams.Waveform.B_valid/2; % (Hz)
% 2-2) Phase noise & nonlinearity parameters
RadarParams.Waveform.nonideal.pn_level = 0.05; %
RadarParams.Waveform.nonideal.f_ripple = 300e3; %
RadarParams.Waveform.nonideal.peak_phase_error = 0.1; %
RadarParams.Waveform.nonideal.power_drop_edge = 0.8; %
RadarParams.Waveform.nonideal.enable_phase_noise = true; % phase noise on/off
RadarParams.Waveform.nonideal.enable_nonlinearity = false; % nonlinearity on/off
% Datasheet phase noise (: -89 dBc/Hz @ 1 MHz)
RadarParams.Waveform.nonideal.phaseNoise.offset_Hz = [1e6];
RadarParams.Waveform.nonideal.phaseNoise.level_dBc_Hz = [-89];
% 2-3) MIMO
RadarParams.Waveform.mimoMode = 'TDM'; % 'TDM' 'DDMA'
% 3) RF - [dBm]
RadarParams.RFOutput.PA_Profile.freqs = [76e9, 76.5e9, 77e9, 77.5e9, 78e9, 79e9];
RadarParams.RFOutput.PA_Profile.power_dBm = [10.5, 11.5, 12.0, 11.5, 10.5, 8.0];
% 4)
% 4-1) (3xN , : m)
RadarParams.Antenna.NumTx = 2;
RadarParams.Antenna.NumRx = 4;
RadarParams.Antenna.lambda = RadarParams.Waveform.lambda_c;
RadarParams.Antenna.TxPos = [ (0:RadarParams.Antenna.NumTx-1)*2*RadarParams.Antenna.lambda; zeros(1,RadarParams.Antenna.NumTx); zeros(1,RadarParams.Antenna.NumTx) ]; % 3xNumTx
RadarParams.Antenna.RxPos = [ (0:RadarParams.Antenna.NumRx-1)*0.5*RadarParams.Antenna.lambda; zeros(1,RadarParams.Antenna.NumRx); zeros(1,RadarParams.Antenna.NumRx) ]; % 3xNumRx
% 4-2) ( NumTx, NumRx )
RadarParams.Antenna.tx_files = {}; % ( )
RadarParams.Antenna.rx_files = {};
RadarParams.Antenna.TxPattern = build_antenna_patterns(RadarParams.Antenna.NumTx, '1D', 12, 5, 30, 10, RadarParams.Antenna.tx_files);
RadarParams.Antenna.RxPattern = build_antenna_patterns(RadarParams.Antenna.NumRx, '2D', 14, 2, 40, 15, RadarParams.Antenna.rx_files);
% 5)
RadarParams.Target.R = [12]; % (m)
RadarParams.Target.v = [7]; % (m/s)
RadarParams.Target.rcs = [-10000]; % RCS (dBsm)
RadarParams.Target.az = [10]; % (deg)
RadarParams.Target.el = [0]; % (deg)
RadarParams.Target.NumTargets = length(RadarParams.Target.R);
% 6)
RadarParams.Rxpath.fc_hpf = 1400; % HPF (Hz)
RadarParams.Rxpath.fc_lpf = 0.8*RadarParams.Waveform.fs_adc/2; % LPF (Hz)
RadarParams.Rxpath.adc_bits = 12; % ADC
RadarParams.Rxpath.adc_v_full_scale = 2.0; % full-scale (Vp-p)
RadarParams.Rxpath.receiver_mode = 'Real'; % 'IQ' 'Real'
RadarParams.Rxpath.load_impedance_ohm = 50; % (Ohm)
RadarParams.Rxpath.rxPathGain_dB = 50; % RX (dB)
RadarParams.Rxpath.system_NF_dB = 15; % (dB)
% 6-5) Spur parameters ( )
% Spur의 , LO leakage, , ADC , .
% SpurParams .
% spur의 , spur apply_lna_and_mixer SpurParams를 spur adc_combined에 .
RadarParams.SpurParams = struct();
RadarParams.SpurParams.lo = struct('amp',0.01,'freq',1e6); %
RadarParams.SpurParams.mixer= struct('alpha2',1e-4,'alpha3',1e-6); % 2/3
RadarParams.SpurParams.adc = struct('amp',1e-3,'freq',2e6); % ADC
RadarParams.SpurParams.switch = struct('amp',5e-4,'freq',500e3); %
RadarParams.SpurParams.pulse = struct('amp',0.0,'freq',200e3,'duty',0.1,'phase',0);%
RadarParams.SpurParams.lo_leak = struct('amp',0.05); % DC/LO leakage
RadarParams.SpurParams.clip = struct('amp',0.03,'range',20); % Clipping spur
RadarParams.SpurParams.enabled = false; % spur on/off
% 7)
RadarParams.SP.RDM.window_type_range = 'hann'; % Range FFT용
RadarParams.SP.RDM.window_type_doppler = 'chebwin'; % Doppler FFT용
num_samples_range = round(RadarParams.Waveform.fs_adc * RadarParams.Waveform.Timing.AdcSampTime);
num_chirps_doppler = RadarParams.Waveform.NumChirps;
% Range
[RadarParams.SP.RDM.window_range, RadarParams.SP.RDM.window_metrics_range] = ...
create_window_with_metrics(RadarParams.SP.RDM.window_type_range, num_samples_range, 60);
% Doppler
[RadarParams.SP.RDM.window_doppler, RadarParams.SP.RDM.window_metrics_doppler] = ...
create_window_with_metrics(RadarParams.SP.RDM.window_type_doppler, num_chirps_doppler, 60);
RadarParams.SP.CFAR.method = 'OS'; % 'CA' 'OS'
RadarParams.SP.CFAR.dimension = '2D'; % '1D' '2D'
RadarParams.SP.CFAR.axis = 'doppler'; % 1D일 : 'range' 'doppler'
RadarParams.SP.CFAR.pfa = 1e-6; % false alarm
RadarParams.SP.CFAR.train = [8, 8]; % [doppler, range] training cell (1D면 )
RadarParams.SP.CFAR.guard = [2, 2]; % [doppler, range] guard cell (1D면 )
RadarParams.SP.CFAR.rank = 0.75; % OS-CFAR rank (0~1)
RadarParams.SP.CFAR.os_scale = 15.0; % OS-CFAR
% 8) Secondary surface loss ( )
RadarParams.Antenna.SecondarySurfaceLoss_dB = 3;
% 9)
RadarParams.Coverage.R_max = 10; % (m)
RadarParams.Coverage.RCS_dBsm = 0; % RCS (dBsm)
%% 2. (TX )
% step 0. TX
Target = RadarParams.Target;
TxPattern = RadarParams.Antenna.TxPattern;
RxPattern = RadarParams.Antenna.RxPattern;
NumTx = RadarParams.Antenna.NumTx;
NumChirps = RadarParams.Waveform.NumChirps;
T_chirp = RadarParams.Waveform.Timing.IdleTime + RadarParams.Waveform.Timing.RampEndTime;
T_frame = T_chirp * NumChirps;
[t, tx_mask] = generate_waveform_timing(RadarParams);
% ================= ADC Raw Data (step 1 ~ 6) =================
% step 1. (Ptx )
TxOut = radiate_antenna(Target, TxPattern);
% step 2. (TTD )
ChannelOut = apply_channel_effects(RadarParams);
% step 3.
RxOut = receive_antenna(ChannelOut, Target, RxPattern);
% step 4-1. (TX_Gain * Space_Loss * RX_Gain)
% RF (V_tx) "시스템 전달 함수" .
SystemAmp = zeros(size(RxOut.space_loss_amp));
for k = 1:Target.NumTargets
for tx = 1:NumTx
SystemAmp(k, :, tx, :, :) = RxOut.space_loss_amp(k, :, tx, :, :) * TxOut.G_tx_amp(k, tx);
end
end
% step 4-2. RF/IF ADC
% - : /, PA , , , MIMO , non-ideal, spur를
% RX .
% - 1) adc_raw_data : [NumRx x NumTx x NumChirps x N_adc_samples]
% TX별로 ( TX )
% - 2) adc_combined : [NumRx x NumChirps x N_adc_samples]
% ADC ( RX에서 TX )
[adc_raw_data, adc_combined] = apply_lna_and_mixer(RxOut, TxOut, RadarParams);
% step 5. ADC Combined ( ADC ) HPF
adc_combined_hpf = apply_analog_hpf(adc_combined, RadarParams.Waveform.fs_adc, RadarParams.Rxpath.fc_hpf);
% step 6. ADC (Quantization & Clipping)
%
v_peak = max(abs(adc_combined_hpf(:)));
lsb_val = RadarParams.Rxpath.adc_v_full_scale / (2^RadarParams.Rxpath.adc_bits);
% ADC
adc_digital = apply_adc_quantization(adc_combined_hpf, RadarParams.Rxpath.adc_bits, RadarParams.Rxpath.adc_v_full_scale, RadarParams.Rxpath.receiver_mode);
% ================= (step 7 ~ ) =================
% step 7. Range-FFT Ideal LPF
% LPF
% (RadarParams.Rxpath.fc_lpf, RadarParams.SP.RDM.*)
% adc_digital을 [NumRx, 1, NumChirps, N_samples] (process_range_fft_lpf )
adc_digital_expanded = reshape(adc_digital, [size(adc_digital,1), 1, size(adc_digital,2), size(adc_digital,3)]);
[range_data, r_axis] = process_range_fft_lpf(adc_digital_expanded, RadarParams);
% step 8. Range Profile (1 , 1 )
if PlotToggle.range_profile
figure('Name', 'Range Profile with Ideal LPF');
plot(r_axis, 20*log10(abs(squeeze(range_data(1,1,1,:)))));
grid on; hold on;
xlabel('Range (m)');
ylabel('Magnitude (dB)');
title(['Range Profile (LPF Cut-off: ', num2str(RadarParams.Rxpath.fc_lpf/1e6), ' MHz)']);
% LPF
xline((RadarParams.Rxpath.fc_lpf * RadarParams.Basic.c)/(2*RadarParams.Waveform.Slope), '--r', 'LPF Cut-off');
end
% step 9. Doppler-FFT
[rd_cube, v_axis] = process_doppler_fft(range_data, RadarParams);
% step 10. NumRx (Noncoherent Integration) RDM
target_rd_map = integrate_nci_rdm(rd_cube, RadarParams); % [Doppler x Range]
% step 11. CFAR (OS/CA, 1D/2D, range/doppler )
[cfar_mask, cfar_threshold, cfar_detections] = detect_targets_cfar(target_rd_map, RadarParams);
%% 3. (Single Chirp & Multi-Chirp)
% [Figure 1]
T_chirp = RadarParams.Waveform.Timing.IdleTime + RadarParams.Waveform.Timing.RampEndTime;
idx_single = (t <= T_chirp);
t_single = t(idx_single);
tx_mask_single = tx_mask(idx_single);
fig_single = [];
if PlotToggle.tx_single
fig_single = visualize_tx_waveform(t_single, RadarParams.Waveform.Timing, RadarParams.Waveform.fc, RadarParams.Waveform.f_start, RadarParams.Waveform.Slope, tx_mask_single, RadarParams.Waveform.nonideal.peak_phase_error, RadarParams.Waveform.nonideal.f_ripple);
end
% [Figure 2] 퀀 MIMO
% (: mimoMode와 NumTx )
fig_multi = [];
if PlotToggle.tx_multi
fig_multi = visualize_multi_tx_waveform(t, RadarParams.Waveform.Timing, RadarParams.Waveform.fc, RadarParams.Waveform.f_start, RadarParams.Waveform.Slope, tx_mask, NumChirps, RadarParams.Waveform.mimoMode, NumTx);
end
% --- [ Figure 3, 4: ] ---
% 'TX' 'RX' .
fig_tx_ant = [];
if PlotToggle.tx_antenna
fig_tx_ant = visualize_antenna_pattern(TxPattern, 'TX');
end
fig_rx_ant = [];
if PlotToggle.rx_antenna
fig_rx_ant = visualize_antenna_pattern(RxPattern, 'RX');
end
% [Figure 5] Range-Doppler Map (NumRx NCI RDM)
fig_rd_map = [];
if PlotToggle.rd_map
fig_rd_map = visualize_rd_map_with_spurs(target_rd_map, r_axis, v_axis, Target, RadarParams.SpurParams, RadarParams, RadarParams.Waveform.Slope);
end
% [Figure 6] CFAR (2D MAP )
fig_cfar = [];
if PlotToggle.cfar
fig_cfar = figure('Name', 'CFAR Detections on NCI RDM');
imagesc(r_axis, v_axis, 20*log10(abs(target_rd_map)));
axis xy;
colormap(jet);
colorbar;
xlabel('Range (m)');
ylabel('Velocity (m/s)');
title('CFAR Detections (Circle Markers)');
hold on;
if ~isempty(cfar_detections)
det_d_idx = cfar_detections(:,1); % Doppler bin index
det_r_idx = cfar_detections(:,2); % Range bin index
det_r = r_axis(det_r_idx);
det_v = v_axis(det_d_idx);
plot(det_r, det_v, 'wo', 'MarkerSize', 7, 'LineWidth', 1.5);
end
hold off;
end
% ================= Coverage Analysis (Performance Analysis) =================
% [Figure 7] Angular Coverage : 2D SNR
% 100m, RCS 0 dBsm SNR
snr_coverage_2d = [];
azimuth_deg = [];
elevation_deg = [];
coverage_info = [];
fig_coverage = [];
if PlotToggle.coverage
[snr_coverage_2d, azimuth_deg, elevation_deg, coverage_info, fig_coverage] = analyze_coverage(RadarParams, RadarParams.Coverage.R_max, RadarParams.Coverage.RCS_dBsm);
end