From 775668afdb1bf71eb3f6349b2c0c0c09cdfcf15f Mon Sep 17 00:00:00 2001 From: YKG Date: Mon, 2 Mar 2026 14:48:21 +0900 Subject: [PATCH] Initial ARSS. Need to improve and add functions. --- 01. Tx Path/apply_mimo_mode.m | 31 +++ 01. Tx Path/generate_waveform_timing.m | 26 ++ 01. Tx Path/radiate_antenna.m | 25 ++ 01. Tx Path/visualize_multi_tx_waveform.m | 145 +++++++++++ 01. Tx Path/visualize_tx_waveform.m | 78 ++++++ 02. Channel/apply_channel_effects.m | 59 +++++ 03. Rx Path/apply_adc_quantization.m | 39 +++ 03. Rx Path/apply_analog_hpf.m | 39 +++ 03. Rx Path/apply_lna_and_mixer.m | 203 +++++++++++++++ .../generate_phase_noise_from_datasheet.m | 58 +++++ 03. Rx Path/generate_spur.m | 83 ++++++ 03. Rx Path/receive_antenna.m | 21 ++ 04. Signal Processing/detect_targets_cfar.m | 140 ++++++++++ 04. Signal Processing/integrate_nci_rdm.m | 18 ++ 04. Signal Processing/process_doppler_fft.m | 45 ++++ 04. Signal Processing/process_range_fft_lpf.m | 48 ++++ .../visualize_rd_map_with_spurs.m | 159 ++++++++++++ 99. Antenna Pattern/build_antenna_patterns.m | 68 +++++ .../visualize_antenna_pattern.m | 81 ++++++ Main.m | 244 ++++++++++++++++++ 20 files changed, 1610 insertions(+) create mode 100644 01. Tx Path/apply_mimo_mode.m create mode 100644 01. Tx Path/generate_waveform_timing.m create mode 100644 01. Tx Path/radiate_antenna.m create mode 100644 01. Tx Path/visualize_multi_tx_waveform.m create mode 100644 01. Tx Path/visualize_tx_waveform.m create mode 100644 02. Channel/apply_channel_effects.m create mode 100644 03. Rx Path/apply_adc_quantization.m create mode 100644 03. Rx Path/apply_analog_hpf.m create mode 100644 03. Rx Path/apply_lna_and_mixer.m create mode 100644 03. Rx Path/generate_phase_noise_from_datasheet.m create mode 100644 03. Rx Path/generate_spur.m create mode 100644 03. Rx Path/receive_antenna.m create mode 100644 04. Signal Processing/detect_targets_cfar.m create mode 100644 04. Signal Processing/integrate_nci_rdm.m create mode 100644 04. Signal Processing/process_doppler_fft.m create mode 100644 04. Signal Processing/process_range_fft_lpf.m create mode 100644 04. Signal Processing/visualize_rd_map_with_spurs.m create mode 100644 99. Antenna Pattern/build_antenna_patterns.m create mode 100644 99. Antenna Pattern/visualize_antenna_pattern.m create mode 100644 Main.m diff --git a/01. Tx Path/apply_mimo_mode.m b/01. Tx Path/apply_mimo_mode.m new file mode 100644 index 0000000..88fd350 --- /dev/null +++ b/01. Tx Path/apply_mimo_mode.m @@ -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 \ No newline at end of file diff --git a/01. Tx Path/generate_waveform_timing.m b/01. Tx Path/generate_waveform_timing.m new file mode 100644 index 0000000..7255b8f --- /dev/null +++ b/01. Tx Path/generate_waveform_timing.m @@ -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 \ No newline at end of file diff --git a/01. Tx Path/radiate_antenna.m b/01. Tx Path/radiate_antenna.m new file mode 100644 index 0000000..4609b77 --- /dev/null +++ b/01. Tx Path/radiate_antenna.m @@ -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 \ No newline at end of file diff --git a/01. Tx Path/visualize_multi_tx_waveform.m b/01. Tx Path/visualize_multi_tx_waveform.m new file mode 100644 index 0000000..7be1e86 --- /dev/null +++ b/01. Tx Path/visualize_multi_tx_waveform.m @@ -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 \ No newline at end of file diff --git a/01. Tx Path/visualize_tx_waveform.m b/01. Tx Path/visualize_tx_waveform.m new file mode 100644 index 0000000..290beb8 --- /dev/null +++ b/01. Tx Path/visualize_tx_waveform.m @@ -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 \ No newline at end of file diff --git a/02. Channel/apply_channel_effects.m b/02. Channel/apply_channel_effects.m new file mode 100644 index 0000000..f4947d6 --- /dev/null +++ b/02. Channel/apply_channel_effects.m @@ -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 \ No newline at end of file diff --git a/03. Rx Path/apply_adc_quantization.m b/03. Rx Path/apply_adc_quantization.m new file mode 100644 index 0000000..7c7ef11 --- /dev/null +++ b/03. Rx Path/apply_adc_quantization.m @@ -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 \ No newline at end of file diff --git a/03. Rx Path/apply_analog_hpf.m b/03. Rx Path/apply_analog_hpf.m new file mode 100644 index 0000000..bb595d4 --- /dev/null +++ b/03. Rx Path/apply_analog_hpf.m @@ -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 \ No newline at end of file diff --git a/03. Rx Path/apply_lna_and_mixer.m b/03. Rx Path/apply_lna_and_mixer.m new file mode 100644 index 0000000..1ee8d31 --- /dev/null +++ b/03. Rx Path/apply_lna_and_mixer.m @@ -0,0 +1,203 @@ +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; + if isfield(phase_noise_cfg,'enabled') && phase_noise_cfg.enabled && ... + isfield(phase_noise_cfg,'offset_Hz') && isfield(phase_noise_cfg,'level_dBc_Hz') + 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; + 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) 변환 (50옴 시스템 기준) + Vtx_inst = sqrt(10.^((P_inst_dBm - 30) / 10) * 50); + + % --- 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 * 50); + + % --- 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=1→TX1, m=2→TX2, m=3→TX1... + 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 \ No newline at end of file diff --git a/03. Rx Path/generate_phase_noise_from_datasheet.m b/03. Rx Path/generate_phase_noise_from_datasheet.m new file mode 100644 index 0000000..4dce8ce --- /dev/null +++ b/03. Rx Path/generate_phase_noise_from_datasheet.m @@ -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 diff --git a/03. Rx Path/generate_spur.m b/03. Rx Path/generate_spur.m new file mode 100644 index 0000000..3f549e4 --- /dev/null +++ b/03. Rx Path/generate_spur.m @@ -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 (0–1), 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 element‑wise 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 diff --git a/03. Rx Path/receive_antenna.m b/03. Rx Path/receive_antenna.m new file mode 100644 index 0000000..e9c9aca --- /dev/null +++ b/03. Rx Path/receive_antenna.m @@ -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 \ No newline at end of file diff --git a/04. Signal Processing/detect_targets_cfar.m b/04. Signal Processing/detect_targets_cfar.m new file mode 100644 index 0000000..9f0456d --- /dev/null +++ b/04. Signal Processing/detect_targets_cfar.m @@ -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 diff --git a/04. Signal Processing/integrate_nci_rdm.m b/04. Signal Processing/integrate_nci_rdm.m new file mode 100644 index 0000000..d5ea49a --- /dev/null +++ b/04. Signal Processing/integrate_nci_rdm.m @@ -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 diff --git a/04. Signal Processing/process_doppler_fft.m b/04. Signal Processing/process_doppler_fft.m new file mode 100644 index 0000000..da810ca --- /dev/null +++ b/04. Signal Processing/process_doppler_fft.m @@ -0,0 +1,45 @@ +function [rd_map, doppler_axis] = process_doppler_fft(range_profile, RadarParams) + % 입력: + % - range_profile: [NumRx, NumTx, NumChirps, NumRangeBins] + % - RadarParams: 메인 파라미터 구조체 + NumChirps = RadarParams.Waveform.NumChirps; + window_type = RadarParams.SP.RDM.window_type_doppler; + + [~, ~, ~, ~] = size(range_profile); + % 중심 주파수에서의 파장 + lambda = RadarParams.Waveform.lambda_c; % = c/fc + + % 1. 처프 간 반복 주기 (PRI, Pulse Repetition Interval) + T_pri = RadarParams.Waveform.PRI; + + % 2. Doppler-FFT용 윈도우 함수 (사용자 선택 가능) + % 도플러 방향(3번째 차원)으로 사이드로브를 억제합니다. + if strcmpi(window_type, 'none') + win_doppler = ones(1, NumChirps); + elseif strcmpi(window_type, 'hamming') + win_doppler = hamming(NumChirps)'; + elseif strcmpi(window_type, 'blackman') + win_doppler = blackman(NumChirps)'; + elseif strcmpi(window_type, 'hann') + win_doppler = hann(NumChirps)'; + else % default: 'chebwin' + win_doppler = chebwin(NumChirps, 60)'; % 60dB 사이드로브 억제 + 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 \ No newline at end of file diff --git a/04. Signal Processing/process_range_fft_lpf.m b/04. Signal Processing/process_range_fft_lpf.m new file mode 100644 index 0000000..8176a86 --- /dev/null +++ b/04. Signal Processing/process_range_fft_lpf.m @@ -0,0 +1,48 @@ +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; + window_type = RadarParams.SP.RDM.window_type_range; + + [~, ~, ~, N_samples] = size(adc_raw_data); + c = RadarParams.Basic.c; + + % 1. 윈도우 함수 적용 + if strcmpi(window_type, 'none') + win = ones(1, N_samples); + elseif strcmpi(window_type, 'hamming') + win = hamming(N_samples)'; + elseif strcmpi(window_type, 'blackman') + win = blackman(N_samples)'; + elseif strcmpi(window_type, 'chebwin') + win = chebwin(N_samples, 60)'; + else % default: 'hann' + win = hann(N_samples)'; + 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 \ No newline at end of file diff --git a/04. Signal Processing/visualize_rd_map_with_spurs.m b/04. Signal Processing/visualize_rd_map_with_spurs.m new file mode 100644 index 0000000..ca38e09 --- /dev/null +++ b/04. Signal Processing/visualize_rd_map_with_spurs.m @@ -0,0 +1,159 @@ +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 핸들 + + fig = figure('Name', '2D Range-Doppler Map'); + imagesc(r_axis, v_axis, 20*log10(abs(target_rd_map))); + axis xy; % y축 방향(속도)을 위로 정렬 + colorbar; + xlabel('Range (m)'); + ylabel('Velocity (m/s)'); + title('Range-Doppler Map (Single Channel)'); + 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; + enable_datasheet_phase_noise = 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 + if isfield(ni,'phaseNoise') && isstruct(ni.phaseNoise) && isfield(ni.phaseNoise,'enabled') + enable_datasheet_phase_noise = logical(ni.phaseNoise.enabled); + end + + show_nonideal_overlay = show_nonideal_overlay && ... + (enable_phase_noise || enable_nonlinearity || enable_datasheet_phase_noise); + 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 diff --git a/99. Antenna Pattern/build_antenna_patterns.m b/99. Antenna Pattern/build_antenna_patterns.m new file mode 100644 index 0000000..bff9327 --- /dev/null +++ b/99. Antenna Pattern/build_antenna_patterns.m @@ -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 \ No newline at end of file diff --git a/99. Antenna Pattern/visualize_antenna_pattern.m b/99. Antenna Pattern/visualize_antenna_pattern.m new file mode 100644 index 0000000..61bed4e --- /dev/null +++ b/99. Antenna Pattern/visualize_antenna_pattern.m @@ -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 \ No newline at end of file diff --git a/Main.m b/Main.m new file mode 100644 index 0000000..1800d75 --- /dev/null +++ b/Main.m @@ -0,0 +1,244 @@ +% ========================================================================= +% 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)); + +%% =================== 파라미터 입력 =================== + +% 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 = false; % 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]; +RadarParams.Waveform.nonideal.phaseNoise.enabled = false; + +% 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, 80, 120]; % 타겟 거리 (m) +RadarParams.Target.v = [7, -5, 0]; % 타겟 속도 (m/s) +RadarParams.Target.rcs = [10, 5, 20]; % 타겟 RCS (dBsm) +RadarParams.Target.az = [10, -10, 0]; % 방위각 (deg) +RadarParams.Target.el = [0, 0, 5]; % 고각 (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.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용 윈도우 +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 임계 스케일 + + +%% 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번 처프) +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'); + +% 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 = 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); + +% [Figure 2] 다중 처프 프레임 시퀀스 및 MIMO 변조 확인 +% (수정됨: mimoMode와 NumTx 변수를 추가로 전달) +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); + +% --- [ Figure 3, 4: 안테나 방사 패턴 단면 도시 ] --- +% 'TX'와 'RX'라는 이름을 넘겨주어 그래프 타이틀과 범례를 구분합니다. +fig_tx_ant = visualize_antenna_pattern(TxPattern, 'TX'); +fig_rx_ant = visualize_antenna_pattern(RxPattern, 'RX'); + +% [Figure 5] Range-Doppler Map 시각화 (NumRx NCI RDM) +fig_rd_map = visualize_rd_map_with_spurs(target_rd_map, r_axis, v_axis, Target, RadarParams.SpurParams, RadarParams, RadarParams.Waveform.Slope); + +% [Figure 6] CFAR 탐지 결과 시각화 (2D MAP 위 원 마커) +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; + + + +