目录
交付物：parabolic_phase_monopulse.m
交付物：cross_polarization_analysis.m
交付物：polarization_jamming_sim.m
2.1 振幅和差单脉冲系统
2.1.1 双波束交叉与等信号方向
振幅和差单脉冲体制通过两个对称偏轴的接收波束实现高精度角度跟踪，其几何布局构成了现代雷达导引头的核心测角架构。两个部分波束相对于瞄准轴（boresight）分别偏移±θ_s角度，在远场辐射方向图上形成高斯型的功率分布重叠区域。当目标位于等信号方向（electric axis）时，两波束接收的归一化电压幅度相等；当目标偏离等信号方向时，幅度差异产生与偏离角成正比的误差电压。
部分波束偏移角的优化设计遵循信噪比损失与角度灵敏度之间的权衡准则。理论上，偏移角的选择需最大化和差通道的误差斜率k_a（error slope）同时最小化和波束的增益损失。对于3dB波束宽度为θ_BW的高斯型波束，最佳偏移角度经验确定为0.39θ_BW，该准则源于对差波束斜率与和波束增益乘积的最大化求解。设单波束电压方向图为E(θ) = E_0 exp(-2.776θ²/θ_BW²)，则和波束与差波束分别构造为：
Σ(θ)=E(θ−θs)+E(θ+θs)Δ(θ)=E(θ−θs)−E(θ+θs)
在等信号方向附近（θ≈0），利用泰勒展开可得差波束的线性近似：
Δ(θ)≈2θs⋅dθdEθ=θs⋅θ
和波束在瞄准轴处取得最大值：
Σ(0)=2E(θs)=2E0exp(−θBW22.776θs2)
波导桥（waveguide hybrid junction）或魔T（magic-T）网络实现了射频端的和差形成。四端口环形波导桥（rat-race coupler）的 circumference 设计为1.5λ，当两路信号从对称端口输入时，和端口（sum port）输出同相叠加信号，差端口（difference port）输出180°反相信号。该无源网络的隔离度通常优于20dB，确保Σ与Δ通道的幅相一致性。和波束用于发射与距离跟踪，差波束用于角度误差提取，二者的正交性通过波导结构的严格对称性保证。
2.1.2 误差信号提取与AGC归一化
角度误差信号的提取依赖于相位检波器（phase detector）对和差通道信号的比较处理。接收机中频输出的和信号Σ(t)与差信号Δ(t)在相位检波器中混频，其输出电压与两信号的复相关性的实部成正比。设和差通道的相位差为ψ_s - ψ_d（理想情况下应为0°或180°），则相位检波器输出为：
S(θ)=Kpd⋅∣Σ∣2Re{Σ∗⋅Δ}≈Kpd⋅ka⋅θ
其中K_pd为相位检波器增益常数，分母|Σ|²实现自动增益控制（AGC）的归一化功能。该归一化消除了目标距离、雷达截面积（RCS）及大气衰减变化导致的回波功率波动，确保角度误差电压仅与偏离角相关。
自动增益控制环路通过检测和通道信号强度调节中频放大器增益，使和通道输出维持恒定电平。AGC环路的带宽设计需满足稳定性与响应速度的双重约束：带宽过窄导致对快速起伏目标的跟踪滞后，带宽过宽则引入热噪声增益调制。对于导弹导引头应用，AGC时间常数通常设定为毫秒量级，介于距离门往返时间与目标角度闪烁带宽之间。AGC归一化引入的额外噪声系数约为1-2dB，该代价换取了角度跟踪通道的大动态范围（通常>60dB）。
2.1.3 和差通道幅相不平衡校正
实际射频通道的幅相不平衡（amplitude and phase mismatch）是单脉冲测角误差的主要来源。设和通道传递函数为H_Σ(ω) = A_Σ exp(jφ_Σ)，差通道为H_Δ(ω) = A_Δ exp(jφ_Δ)，则归一化误差信号修正为：
Sactual(θ)=AΣAΔ⋅∣Σ∣2Re{Δ⋅Σ∗exp(j(ϕΔ−ϕΣ))}
幅相失配导致误差曲线零点偏移（bias）和斜率变化（scale factor error）。正交混频架构中的I/Q失衡（imbalance）进一步恶化通道匹配度，需通过数字校准算法补偿。
基于最大似然估计的通道校准算法利用注入的校准信号或外部点源目标，估计和差通道的复增益比。设接收信号矢量为x = [x_Σ, x_Δ]^T，理想信号模型为s = [s_Σ, s_Δ]^T，则失配矩阵M的最大似然估计通过最小化代价函数实现：
M^=argMmini=1∑N∥xi−Msi∥2
其中M为2×2复对角矩阵，对角元素分别为和差通道的增益与相位误差。校准后的信号通过逆矩阵M^{-1}重构，恢复理想的和差特性。该算法要求校准源覆盖系统带宽内的多个频点，以捕获色散性失配。
2.2 相位和差单脉冲系统
2.2.1 双天线基线干涉测角
相位和差单脉冲体制通过空间分离的两个天线单元接收信号的相位差实现角度测量，适用于长基线高精度测角场景。两天线间距（基线长度）为L，目标入射波与基线法线夹角为θ，则两天线接收信号的相位差由路径差ΔR = L sin(θ)决定：
Δψ=λ2πLsinθ
该相位差通过共模本振（common LO）下的混频器提取，经相位检波器转换为与角度相关的电压。相位干涉仪的测角精度与基线长度成正比，然而基线增长引入相位模糊（phase ambiguity）问题：当|Δψ| > π时，无法区分真实角度与镜像角度。无模糊测角范围受限于：
L≤2sinθmaxλ
对于宽覆盖应用（如θ_max = 60°），基线长度受限为0.577λ，严重制约测角精度。多基线干涉仪（multi-baseline interferometer）通过长短基线组合解决模糊问题：长基线提供高精度，短基线（通常L_short = λ/2）提供无模糊粗测。中国余数定理（Chinese Remainder Theorem）或最优解卷积算法用于从多组相位测量中重构无模糊角度估计。
2.2.2 抛物面反射器中的相位和差实现
相位和差单脉冲在反射面天线中的实现依赖于孔径面的相位分割。S-200导弹导引头采用的金属隔板（metal septum）将抛物面反射器的孔径沿俯仰或方位维度分为两个区域，两区域的反射场相位相反（180°相位差），形成差通道；全孔径同相合成形成和通道。隔板位于抛物面焦点与反射面之间，其几何位置决定和差波束的交叉电平与隔离度。
基于物理光学（Physical Optics）的孔径场分析采用镜像法（method of images）计算隔板存在时的电流分布。设抛物面焦点位于F，隔板平面将孔径分为区域A与区域B，则区域A的镜面电流J_A与区域B的电流J_B在隔板边缘满足边界条件n×(E_A - E_B) = 0。和波束对应J_A + J_B的辐射积分，差波束对应J_A - J_B。隔板深度d与抛物面焦距f的比值d/f影响差波束的零点深度，优化设计使差波束零点低于-30dB以确保跟踪精度。
2.3 交叉极化辐射与极化干扰
2.3.1 反射器天线交叉极化机理
抛物面反射器天线的交叉极化（cross-polarization）辐射源于曲面几何对入射场极化矢量的旋转作用。当馈源辐射的线极化波经抛物面反射后，除主极化分量（co-polar）外，在45°对角平面产生正交的交叉极化分量。对于对称抛物面，交叉极化场分布呈现奇对称特性，在瞄准轴处为零，在约0.7倍波束宽度处达到峰值。
主极化与交叉极化方向图通过 Ludwig-3 定义区分。设馈源辐射电场为E_f = E_θ θ̂ + E_φ φ̂，经抛物面反射后的远场通过表面电流积分计算：
Erad=4πrjωμe−jkr∫S[J−(J⋅r^)r^]ejk⋅r′dS′
交叉极化电平（XPD）定义为峰值交叉极化增益与峰值主极化增益之比，高性能单脉冲雷达要求XPD优于-40dB。焦距与口径比（f/D）是设计关键：f/D增大降低交叉极化但增加结构尺寸，通常选择f/D在0.3-0.4之间平衡电气与机械约束。
2.3.2 极化干扰对单脉冲的影响
极化失配（polarization mismatch）对单脉冲跟踪的影响体现在和差通道的能量耦合与相位畸变。当干扰信号（jamming）的极化状态与雷达接收极化正交时，理想情况下和通道接收为零；然而反射器天线的交叉极化响应导致干扰能量泄漏至差通道，产生虚假角度误差。
线极化干扰作用于圆极化雷达时，极化失配损失3dB，但更严重的是交叉极化方向图的畸变（squint）现象。极化旋转使差波束的零值点偏离瞄准轴，误差曲线斜率在某些极化态下变号，导致跟踪环路收敛至错误角度甚至发散（即S(φ)→∞的失稳判据）。极化自适应（polarization adaptation）或极化分集接收是抑制此类干扰的有效对策。
第二部分 代码实现
交付物：sum_difference_patterns.m
脚本功能：生成高斯型振幅和差单脉冲波束方向图，实现0.39θ_BW最优偏移角计算，输出和波束增益、差波束斜率与和差损失曲线，支持波束宽度归一化显示。
使用方式：设定波束宽度参数与偏移角扫描范围，直接运行脚本生成方向图与性能曲线。修改theta_s_frac可验证不同偏移角（0.2-0.5倍波束宽度）对斜率与损失的影响。
matlab
复制
function sum_difference_patterns()
% sum_difference_patterns: Amplitude Monopulse Sum/Difference Beam Pattern Analysis
% Optimizes squint angle for maximum error slope with minimum sum loss
% Implements 0.39*theta_BW design criterion
%
% Usage:
% Run sum_difference_patterns() to execute complete analysis
% Adjust beamwidth_theta and squint_range for parametric studies
%
% Outputs:
% - Figure 1: Normalized Sum and Difference beam patterns
% - Figure 2: Error slope vs Squint angle (optimal at ~0.39)
% - Figure 3: Sum pattern loss vs Squint angle
% System Parameters
beamwidth_theta = 3.0; % 3dB beamwidth (degrees)
squint_optimal = 0.39 * beamwidth_theta; % Optimal squint angle (degrees)
% Angular grid (normalized to beamwidth)
theta_norm = linspace(-1.5, 1.5, 1000); % +/- 1.5 beamwidths
theta_deg = theta_norm * beamwidth_theta; % Absolute angle in degrees
% Gaussian beam pattern (voltage)
% E(theta) = exp(-2.776 * (theta/theta_BW)^2) for power pattern -3dB at theta=theta_BW/2
% For voltage pattern: -3dB corresponds to 1/sqrt(2), so exponent is -1.386*(theta/(0.5*theta_BW))^2
% Standard form: E(theta) = exp(-4*ln(2)*(theta/theta_BW)^2) = exp(-2.773*(theta/theta_BW)^2)
E_single = @(theta, theta_s) exp(-2.773 * ((theta - theta_s)/beamwidth_theta).^2);
% Generate patterns for optimal squint
E_left = E_single(theta_deg, -squint_optimal);
E_right = E_single(theta_deg, squint_optimal);
% Sum and Difference beams
Sigma = E_left + E_right;
Delta = E_left - E_right;
% Normalize
Sigma_norm = Sigma / max(Sigma);
Delta_norm = Delta / max(abs(Delta));
% Figure 1: Beam Patterns
figure('Name', 'Amplitude Monopulse Beam Patterns', 'Position', [100 100 1400 600]);
subplot(1, 2, 1);
plot(theta_norm, abs(E_left), 'b--', 'LineWidth', 1.5); hold on;
plot(theta_norm, abs(E_right), 'b--', 'LineWidth', 1.5);
plot(theta_norm, abs(Sigma_norm), 'r-', 'LineWidth', 2);
plot(theta_norm, abs(Delta_norm), 'g-', 'LineWidth', 2);
xline(0, 'k:', 'LineWidth', 1);
title(sprintf('Normalized Beam Patterns (Squint = %.2f° = 0.39\\theta_{BW})', squint_optimal));
xlabel('Angle (\theta/\theta_{BW})');
ylabel('Normalized Voltage');
legend('Left Beam', 'Right Beam', 'Sum (\Sigma)', 'Difference (\Delta)', 'Boresight');
grid on; xlim([-1.5 1.5]); ylim([0 1.1]);
% Detail view near boresight
subplot(1, 2, 2);
idx = abs(theta_norm) < 0.3; % Zoom to +/- 0.3 beamwidths
plot(theta_norm(idx), abs(Sigma_norm(idx)), 'r-', 'LineWidth', 2); hold on;
plot(theta_norm(idx), abs(Delta_norm(idx)), 'g-', 'LineWidth', 2);
plot(theta_norm(idx), angle(Delta_norm(idx))*(180/pi)/90, 'm--', 'LineWidth', 1.5);
title('Detail Near Boresight');
xlabel('Angle (\theta/\theta_{BW})');
ylabel('Normalized Amplitude / Phase (deg/90)');
legend('|\Sigma|', '|\Delta|', '\angle\Delta (scaled)');
grid on;
% Figure 2: Parametric Analysis of Squint Angle
figure('Name', 'Squint Angle Optimization', 'Position', [150 150 1200 900]);
squint_range = linspace(0.1, 0.6, 50) * beamwidth_theta; % 0.1 to 0.6 beamwidths
error_slopes = zeros(size(squint_range));
sum_gains = zeros(size(squint_range));
sum_losses_dB = zeros(size(squint_range));
for i = 1:length(squint_range)
theta_s = squint_range(i);
El = E_single(theta_deg, -theta_s);
Er = E_single(theta_deg, theta_s);
Sig = El + Er;
Dif = El - Er;
% Calculate error slope at boresight (theta=0)
% Slope k = d(Delta)/d(theta) at theta=0 normalized by Sum(0)
dtheta = theta_deg(2) - theta_deg(1);
% Central difference approximation for derivative at boresight
idx0 = find(abs(theta_deg) < 0.01, 1, 'first');
if isempty(idx0), idx0 = floor(length(theta_deg)/2); end
% Linear fit near boresight for robust slope estimation
idx_fit = abs(theta_deg) < 0.1 * beamwidth_theta;
p_Delta = polyfit(theta_deg(idx_fit), abs(Dif(idx_fit)), 1);
error_slopes(i) = p_Delta(1) / max(Sig); % Normalized slope
% Sum pattern gain (relative to single beam peak = 1)
sum_gains(i) = max(Sig);
sum_losses_dB(i) = -20*log10(sum_gains(i)/2); % Loss relative to ideal sum of 2
end
% Normalize slope to maximum
slopes_norm = error_slopes / max(error_slopes);
% Subplot 1: Error Slope vs Squint
subplot(2, 2, 1);
plot(squint_range/beamwidth_theta, slopes_norm, 'b-', 'LineWidth', 2);
hold on;
[max_slope, idx_max] = max(slopes_norm);
plot(squint_range(idx_max)/beamwidth_theta, max_slope, 'ro', 'MarkerSize', 10, 'MarkerFaceColor', 'r');
xline(0.39, 'g--', 'LineWidth', 2, 'Label', '0.39 criterion');
title('Normalized Error Slope vs Squint Angle');
xlabel('Squint Angle (\theta_s/\theta_{BW})');
ylabel('Normalized Slope k_a');
grid on;
text(squint_range(idx_max)/beamwidth_theta*0.8, 0.9, sprintf('Optimum: %.3f', squint_range(idx_max)/beamwidth_theta));
% Subplot 2: Sum Pattern Loss
subplot(2, 2, 2);
plot(squint_range/beamwidth_theta, sum_losses_dB, 'r-', 'LineWidth', 2);
hold on;
xline(0.39, 'g--', 'LineWidth', 2);
title('Sum Pattern Loss vs Squint Angle');
xlabel('Squint Angle (\theta_s/\theta_{BW})');
ylabel('Loss (dB)');
grid on;
text(0.5, max(sum_losses_dB)*0.8, sprintf('Loss at 0.39: %.2f dB', interp1(squint_range/beamwidth_theta, sum_losses_dB, 0.39)));
% Subplot 3: Product of Slope and Sum Gain (Figure of Merit)
subplot(2, 2, 3);
FOM = slopes_norm .* (sum_gains/2); % Trade-off metric
plot(squint_range/beamwidth_theta, FOM, 'm-', 'LineWidth', 2);
hold on;
[max_FOM, idx_FOM] = max(FOM);
plot(squint_range(idx_FOM)/beamwidth_theta, max_FOM, 'ko', 'MarkerSize', 10, 'MarkerFaceColor', 'k');
xline(0.39, 'g--', 'LineWidth', 2);
title('Figure of Merit: Slope \times Sum Gain');
xlabel('Squint Angle (\theta_s/\theta_{BW})');
ylabel('Normalized FOM');
grid on;
text(squint_range(idx_FOM)/beamwidth_theta, max_FOM*0.9, sprintf('Max FOM at %.3f', squint_range(idx_FOM)/beamwidth_theta));
% Subplot 4: Difference Pattern Null Depth
subplot(2, 2, 4);
null_depths = zeros(size(squint_range));
for i = 1:length(squint_range)
theta_s = squint_range(i);
El = E_single(theta_deg, -theta_s);
Er = E_single(theta_deg, theta_s);
Dif = El - Er;
% Find minimum near boresight
idx_null = abs(theta_deg) < 0.05;
null_depths(i) = 20*log10(min(abs(Dif(idx_null))) / max(abs(Dif)));
end
plot(squint_range/beamwidth_theta, null_depths, 'k-', 'LineWidth', 2);
title('Difference Pattern Null Depth');
xlabel('Squint Angle (\theta_s/\theta_{BW})');
ylabel('Null Depth (dB)');
grid on;
ylim([-100 0]);
% Summary output
fprintf('\n=== Amplitude Monopulse Analysis Summary ===\n');
fprintf('Beamwidth: %.2f degrees\n', beamwidth_theta);
fprintf('Optimal squint angle: %.3f * beamwidth = %.3f degrees\n', squint_range(idx_max)/beamwidth_theta, squint_range(idx_max));
fprintf('Normalized error slope at optimum: %.4f\n', max_slope);
fprintf('Sum pattern loss at 0.39 criteria: %.2f dB\n', interp1(squint_range/beamwidth_theta, sum_losses_dB, 0.39));
fprintf('Difference null depth: %.1f dB\n', interp1(squint_range/beamwidth_theta, null_depths, 0.39));
fprintf('===========================================\n');
end交付物：agc_normalization.m
脚本功能：实现自动增益控制（AGC）环路的离散时间仿真，模拟和差单脉冲接收机中AGC对角度误差信号归一化的动态过程，分析不同输入功率下的收敛时间与稳态误差。
使用方式：设定目标距离剖面（range profile）与AGC环路带宽，运行脚本观察AGC增益调整、归一化误差输出及瞬态响应。支持阶跃输入与动态起伏目标模拟。
matlab
复制
function agc_normalization()
% agc_normalization: AGC Loop Simulation for Monopulse Angle Tracking
% Models automatic gain control dynamics and normalization effects
%
% Usage:
% Run agc_normalization() for standard scenario
% Modify agc_tau for different loop bandwidths (0.1ms to 10ms)
%
% Outputs:
% - Time-domain AGC response to power variations
% - Normalized angle error output vs true error
% - Convergence time analysis for different input levels
% Simulation Parameters
fs = 1e6; % Sampling rate (1 MHz, typical IF bandwidth)
T_sim = 0.05; % Simulation time (50 ms)
t = 0:1/fs:T_sim;
N = length(t);
% AGC Loop Parameters
agc_tau = 0.002; % AGC time constant (seconds) - 2ms typical
agc_alpha = 1 - exp(-1/(fs*agc_tau)); % Digital loop gain
agc_setpoint = 1.0; % Desired output level (normalized)
% Scenario: Target with varying range (power changes) and constant angle error
% Range profile: starts at 10km, closes to 5km, then steps change
range_km = 10 - 4*(t/T_sim) + 0.5*randn(1,N); % Closing target with noise
range_km(range_km < 2) = 2; % Minimum range limit
% Received power (4th power law + fluctuations)
P_received = (10./range_km).^4; % Normalized power (20km reference)
P_received = P_received * (1 + 0.2*sin(2*pi*50*t)); % 50Hz scintillation
% True angle error (constant 0.1 beamwidths)
theta_error_true = 0.1; % Constant 10% of beamwidth
beamwidth = 3.0; % degrees
% Monopulse error slope (linear region)
k_error = 1.5; % Error voltage per beamwidth
% Unnormalized error signal (before AGC)
error_voltage_raw = k_error * theta_error_true * sqrt(P_received);
% AGC Loop Simulation (First-order feedback loop)
agc_gain = ones(1, N); % Initial gain
vga_out = zeros(1, N); % Variable gain amplifier output
normalized_error = zeros(1, N);
for n = 2:N
% VGA output (input * current gain)
vga_out(n) = sqrt(P_received(n)) * agc_gain(n-1);
% AGC Detector (envelope detection, simplified)
detected = abs(vga_out(n));
% Loop filter (integrator with leakage)
agc_gain(n) = agc_gain(n-1) + agc_alpha * (agc_setpoint - detected) / detected;
% Prevent negative gain
if agc_gain(n) < 0.01, agc_gain(n) = 0.01; end
if agc_gain(n) > 100, agc_gain(n) = 100; end
% Normalized angle error (error_voltage / sum_amplitude)
% Sum amplitude is approximated by VGA output (assuming sum channel controlled)
sum_amplitude = vga_out(n);
if sum_amplitude > 0.01
normalized_error(n) = error_voltage_raw(n) / sum_amplitude;
else
normalized_error(n) = 0;
end
end
% Figure 1: Time Domain Response
figure('Name', 'AGC Loop Dynamics', 'Position', [100 100 1400 900]);
subplot(4, 1, 1);
plot(t*1000, range_km, 'b-', 'LineWidth', 1.5);
title('Target Range Profile');
xlabel('Time (ms)');
ylabel('Range (km)');
grid on;
subplot(4, 1, 2);
plot(t*1000, 10*log10(P_received), 'r-', 'LineWidth', 1.5);
title('Received Power (dB)');
xlabel('Time (ms)');
ylabel('Power (dB)');
grid on;
subplot(4, 1, 3);
plot(t*1000, 20*log10(agc_gain), 'g-', 'LineWidth', 1.5);
title('AGC Gain (dB)');
xlabel('Time (ms)');
ylabel('Gain (dB)');
grid on;
subplot(4, 1, 4);
plot(t*1000, normalized_error, 'm-', 'LineWidth', 1.5);
hold on;
yline(k_error * theta_error_true, 'r--', 'LineWidth', 2, 'Label', 'Ideal Error');
title('Normalized Angle Error Output');
xlabel('Time (ms)');
ylabel('Normalized Error');
grid on;
legend('AGC Output', 'True Error');
% Figure 2: Convergence Analysis for Different Power Levels
figure('Name', 'AGC Convergence Analysis', 'Position', [150 150 1200 600]);
power_steps = [0.01, 0.1, 1.0, 10]; % Different input power levels
colors = lines(length(power_steps));
subplot(1, 2, 1);
for p_idx = 1:length(power_steps)
P_in = power_steps(p_idx);
agc_gain_step = ones(1, N);
for n = 2:N
detected = P_in * agc_gain_step(n-1);
agc_gain_step(n) = agc_gain_step(n-1) + agc_alpha * (agc_setpoint - detected) / detected;
if agc_gain_step(n) < 0.01, agc_gain_step(n) = 0.01; end
end
plot(t*1000, 20*log10(agc_gain_step), 'Color', colors(p_idx,:), 'LineWidth', 2);
hold on;
end
title('AGC Gain Convergence for Different Input Powers');
xlabel('Time (ms)');
ylabel('AGC Gain (dB)');
legend(arrayfun(@(x) sprintf('P_{in} = %.2f', x), power_steps, 'UniformOutput', false), 'Location', 'east');
grid on;
% Settling time calculation
subplot(1, 2, 2);
settling_times = zeros(size(power_steps));
for p_idx = 1:length(power_steps)
P_in = power_steps(p_idx);
agc_gain_step = ones(1, N);
settled = false;
for n = 2:N
detected = P_in * agc_gain_step(n-1);
agc_gain_step(n) = agc_gain_step(n-1) + agc_alpha * (agc_setpoint - detected) / detected;
% Check settling (within 5% of final)
final_gain = agc_setpoint / P_in;
if ~settled && abs(agc_gain_step(n) - final_gain) < 0.05*final_gain && n > 10
settling_times(p_idx) = t(n) * 1000; % ms
settled = true;
end
end
if ~settled, settling_times(p_idx) = T_sim * 1000; end
end
bar(power_steps, settling_times);
title('AGC Settling Time vs Input Power');
xlabel('Input Power Level (normalized)');
ylabel('Settling Time (ms)');
set(gca, 'XScale', 'log');
grid on;
% Figure 3: Angle Error vs Power Level (AGC vs No AGC comparison)
figure('Name', 'AGC Normalization Performance', 'Position', [200 200 1000 500]);
power_range = logspace(-2, 2, 100); % 0.01 to 100
theta_true = 0.15; % Beamwidths
% Without AGC (raw error)
error_raw = k_error * theta_true * sqrt(power_range);
% With AGC (normalized)
% After settling, VGA output is constant, so normalized error is constant
error_agc = k_error * theta_true * ones(size(power_range));
subplot(1, 2, 1);
semilogx(power_range, error_raw, 'b-', 'LineWidth', 2); hold on;
semilogx(power_range, error_agc, 'r--', 'LineWidth', 2);
title('Angle Error: Raw vs AGC Normalized');
xlabel('Received Power (normalized)');
ylabel('Error Voltage');
legend('Without AGC', 'With AGC');
grid on;
% AGC induced noise (gain modulation noise)
subplot(1, 2, 2);
% Thermal noise through AGC
k_boltz = 1.38e-23;
T0 = 290;
B = 1e6; % Bandwidth
noise_figure = 3; % dB
F_lin = 10^(noise_figure/10);
noise_power = k_boltz * T0 * B * F_lin;
% Signal-to-noise ratio
SNR = power_range / noise_power;
% AGC noise figure degradation (simplified model)
agc_noise_factor = 1 + 0.5./SNR; % AGC introduces ~0.5/SNR additional noise
agc_noise_dB = 10*log10(agc_noise_factor);
semilogx(SNR, agc_noise_dB, 'k-', 'LineWidth', 2);
title('AGC Noise Figure Degradation');
xlabel('Input SNR (linear)');
ylabel('Additional Noise (dB)');
grid on;
fprintf('\n=== AGC Normalization Analysis ===\n');
fprintf('AGC Time Constant: %.2f ms\n', agc_tau*1000);
fprintf('Loop Bandwidth (approx): %.2f Hz\n', 1/(2*pi*agc_tau));
fprintf('Settling time for 10x power step: %.2f ms\n', interp1(power_steps, settling_times, 10));
fprintf('=================================\n');
end交付物：channel_calibration.m
脚本功能：实现和差通道幅相不平衡的离线校准与实时补偿，基于最大似然估计与正交I/Q数据，估计通道间复增益误差并重构理想和差信号。
使用方式：生成或加载包含校准源（点目标或注入信号）的原始和差通道IQ数据，调用校准函数估计失配矩阵，应用补偿系数输出校正后的和差波束。
matlab
复制
function channel_calibration()
% channel_calibration: Amplitude and Phase Mismatch Correction for Sum/Difference Channels
% Implements maximum likelihood estimation of channel imbalances
% Supports I/Q calibration and quadrature error compensation
%
% Usage:
% Run channel_calibration() for demonstration with synthetic data
% Input: raw_sum_IQ, raw_delta_IQ matrices (N x num_pulses)
% Output: calibrated sum/delta with imbalance metrics
%
% Features:
% - Complex gain ratio estimation
% - I/Q imbalance correction (gain and orthogonality)
% - Calibration signal injection simulation
% Generate synthetic calibration data
% Scenario: Point target at boresight, varying SNR, with known imbalances
N_samples = 1000; % Range bins
N_pulses = 64; % Coherent processing interval
% True channel characteristics (unknown to algorithm, ground truth)
true_sum_gain = 1.0;
true_delta_gain = 0.95; % 0.5dB amplitude imbalance
true_delta_phase = 8 * pi/180; % 8 degrees phase offset
true_I_phase = 2 * pi/180; % I/Q orthogonality error in sum channel
true_Q_gain = 1.02; % I/Q amplitude imbalance
% Generate ideal signals (calibration source at boresight)
% Sum channel: strong signal
SNR_sum_dB = 20;
signal_sum = ones(N_samples, 1) * sqrt(10^(SNR_sum_dB/10));
noise_sum = (randn(N_samples, N_pulses) + 1j*randn(N_samples, N_pulses))/sqrt(2);
% Delta channel: near null (small leakage), then add imbalance
null_depth = -30; % dB
signal_delta_leakage = signal_sum * 10^(null_depth/20);
% Add synthetic noise and imbalances
raw_sum = (true_sum_gain * signal_sum * (1+randn(1,N_pulses)*0.01)) + noise_sum;
raw_delta = (true_delta_gain * exp(1j*true_delta_phase) * signal_sum * 10^(null_depth/20) * (1+randn(1,N_pulses)*0.01)) + ...
noise_sum * 0.8; % Slightly different noise realization
% Apply I/Q imbalances to sum channel (crosstalk between I and Q)
raw_sum_I = real(raw_sum) * cos(true_I_phase) - imag(raw_sum) * sin(true_I_phase);
raw_sum_Q = real(raw_sum) * sin(true_I_phase) + imag(raw_sum) * cos(true_I_phase) * true_Q_gain;
raw_sum = raw_sum_I + 1j*raw_sum_Q;
% Calibration Algorithm 1: Maximum Likelihood Estimation of Complex Gain Ratio
% Model: Delta = rho * Sum + noise, where rho is complex imbalance factor
% ML Estimate: rho = E[Delta * conj(Sum)] / E[|Sum|^2]
% Calculate correlation matrices
R_sum_sum = mean(mean(raw_sum .* conj(raw_sum)));
R_delta_sum = mean(mean(raw_delta .* conj(raw_sum)));
% Estimate imbalance
rho_est = R_delta_sum / R_sum_sum;
gain_imbalance_dB = 20*log10(abs(rho_est));
phase_imbalance_deg = angle(rho_est) * 180/pi;
% Compensate Delta channel
delta_compensated = raw_delta / rho_est;
% Calibration Algorithm 2: I/Q Orthogonality and Gain Balance
% Separate I/Q components
sum_I = real(raw_sum);
sum_Q = imag(raw_sum);
% Estimate I/Q correlation (should be zero for perfect orthogonality)
IQ_corr = mean(sum_I(:) .* sum_Q(:));
I_var = var(sum_I(:));
Q_var = var(sum_Q(:));
IQ_gain_ratio = sqrt(Q_var / I_var);
IQ_phase_err = asin(IQ_corr / sqrt(I_var * Q_var)) * 180/pi;
% Correct I/Q (Gram-Schmidt orthogonalization)
sum_I_corr = sum_I;
sum_Q_corr = (sum_Q - mean(sum_I(:).*sum_Q(:))/mean(sum_I(:).^2) * sum_I) / IQ_gain_ratio;
sum_compensated = sum_I_corr + 1j*sum_Q_corr;
% Figure 1: I/Q Constellation Diagrams
figure('Name', 'Channel Calibration Analysis', 'Position', [100 100 1400 900]);
% Raw Sum I/Q
subplot(3, 2, 1);
scatter(real(raw_sum(:)), imag(raw_sum(:)), 10, 'b', 'filled', 'MarkerFaceAlpha', 0.3);
axis equal; grid on;
title(sprintf('Raw Sum I/Q\\nGain Err: %.2f dB, Phase Err: %.2f°', 20*log10(IQ_gain_ratio), IQ_phase_err));
xlabel('I (In-phase)'); ylabel('Q (Quadrature)');
% Corrected Sum I/Q
subplot(3, 2, 2);
scatter(real(sum_compensated(:)), imag(sum_compensated(:)), 10, 'r', 'filled', 'MarkerFaceAlpha', 0.3);
axis equal; grid on;
title('Corrected Sum I/Q (Orthogonalized)');
xlabel('I'); ylabel('Q');
% Raw Delta vs Sum (phase relationship)
subplot(3, 2, 3);
scatter(real(raw_sum(:)), real(raw_delta(:)), 10, 'b', 'filled', 'MarkerFaceAlpha', 0.3);
hold on;
scatter(imag(raw_sum(:)), imag(raw_delta(:)), 10, 'g', 'filled', 'MarkerFaceAlpha', 0.3);
grid on;
title(sprintf('Raw Sum vs Delta Correlation\\nEstimated Gain Imabalance: %.2f dB', gain_imbalance_dB));
xlabel('Sum Component'); ylabel('Delta Component');
legend('In-phase', 'Quadrature');
% Corrected Delta vs Sum
subplot(3, 2, 4);
scatter(real(sum_compensated(:)), real(delta_compensated(:)), 10, 'r', 'filled', 'MarkerFaceAlpha', 0.3);
hold on;
scatter(imag(sum_compensated(:)), imag(delta_compensated(:)), 10, 'm', 'filled', 'MarkerFaceAlpha', 0.3);
grid on;
title('Compensated Sum vs Delta (Orthogonal)');
xlabel('Sum'); ylabel('Delta');
% Angle error comparison (simulated off-boresight target)
subplot(3, 2, 5);
theta_scan = linspace(-0.5, 0.5, 100); % Beamwidths
% Simulate response with and without calibration
error_raw = zeros(size(theta_scan));
error_cal = zeros(size(theta_scan));
for i = 1:length(theta_scan)
% Synthetic off-boresight signal
sig = exp(-2.773*theta_scan(i)^2); % Gaussian beam
% Raw (with imbalance)
delta_raw_sim = true_delta_gain * exp(1j*true_delta_phase) * sig * theta_scan(i) * 1.5; % error slope
sum_raw_sim = sig;
error_raw(i) = real(delta_raw_sim / sum_raw_sim);
% Calibrated
delta_cal_sim = delta_raw_sim / rho_est;
sum_cal_sim = sum_raw_sim;
error_cal(i) = real(delta_cal_sim / sum_cal_sim);
end
plot(theta_scan, error_raw, 'b-', 'LineWidth', 2); hold on;
plot(theta_scan, error_cal, 'r--', 'LineWidth', 2);
plot(theta_scan, 1.5*theta_scan, 'k:', 'LineWidth', 1.5);
title('Monopulse Error Curve: Raw vs Calibrated');
xlabel('Angle (beamwidths)');
ylabel('Error Output');
legend('Raw (Imbalanced)', 'Calibrated', 'Ideal Linear');
grid on;
% Error statistics
subplot(3, 2, 6);
bias_raw = mean(error_raw - 1.5*theta_scan');
bias_cal = mean(error_cal - 1.5*theta_scan');
slope_raw = polyfit(theta_scan, error_raw, 1);
slope_cal = polyfit(theta_scan, error_cal, 1);
bar_data = [slope_raw(1), slope_cal(1); bias_raw, bias_cal];
bar(bar_data);
set(gca, 'XTickLabel', {'Slope', 'Bias'});
legend('Raw', 'Calibrated');
title('Error Curve Parameters');
ylabel('Value');
grid on;
% Figure 2: Calibration Metrics Table
figure('Name', 'Calibration Metrics', 'Position', [200 200 800 400]);
axis off;
metrics = {
'Parameter', 'True Value', 'Estimated Value', 'Error';
'Amp Imbalance (dB)', sprintf('%.2f', 20*log10(true_delta_gain)), sprintf('%.2f', gain_imbalance_dB), sprintf('%.2f', abs(20*log10(true_delta_gain) - gain_imbalance_dB));
'Phase Imbalance (deg)', sprintf('%.2f', true_delta_phase*180/pi), sprintf('%.2f', phase_imbalance_deg), sprintf('%.2f', abs(true_delta_phase*180/pi - phase_imbalance_deg));
'I/Q Gain Ratio', sprintf('%.3f', true_Q_gain), sprintf('%.3f', IQ_gain_ratio), sprintf('%.3f', abs(true_Q_gain - IQ_gain_ratio));
'I/Q Phase Error (deg)', sprintf('%.2f', true_I_phase*180/pi), sprintf('%.2f', IQ_phase_err), sprintf('%.2f', abs(true_I_phase*180/pi - IQ_phase_err));
'Error Slope', '1.50', sprintf('%.2f', slope_cal(1)), sprintf('%.2f', abs(1.5 - slope_cal(1)));
'Zero Bias', '0.00', sprintf('%.4f', bias_cal), sprintf('%.4f', abs(bias_cal))
};
y_pos = 0.9;
for i = 1:size(metrics, 1)
if i == 1
font_weight = 'bold';
font_size = 12;
else
font_weight = 'normal';
font_size = 11;
end
line_str = sprintf('%-20s %-15s %-15s %-10s', metrics{i,1}, metrics{i,2}, metrics{i,3}, metrics{i,4});
text(0.05, y_pos, line_str, 'FontName', 'Courier', 'FontSize', font_size, 'FontWeight', font_weight);
y_pos = y_pos - 0.15;
end
fprintf('\n=== Channel Calibration Summary ===\n');
fprintf('Estimated Gain Imbalance: %.2f dB (True: %.2f dB)\n', gain_imbalance_dB, 20*log10(true_delta_gain));
fprintf('Estimated Phase Imbalance: %.2f degrees (True: %.2f degrees)\n', phase_imbalance_deg, true_delta_phase*180/pi);
fprintf('I/Q Gain Ratio: %.3f (True: %.3f)\n', IQ_gain_ratio, true_Q_gain);
fprintf('Error Slope Correction: %.3f -> %.3f\n', slope_raw(1), slope_cal(1));
fprintf('====================================\n');
end交付物：phase_interferometer.m
脚本功能：实现双基线相位干涉仪测角系统，包含长/短基线配置下的相位差计算、相位模糊解算算法（基于中国余数定理），以及基线长度优化分析。
使用方式：设定波长、长短基线长度与目标角度范围，运行脚本生成相位差曲线、解模糊成功概率与测角精度随信噪比变化曲线。支持蒙特卡洛仿真验证解算鲁棒性。
matlab
复制
function phase_interferometer()
% phase_interferometer: Dual-Baseline Phase Interferometer Angle Estimation
% Implements phase ambiguity resolution using Chinese Remainder Theorem
% Analyzes baseline length constraints and angle estimation precision
%
% Usage:
% Run phase_interferometer() for demonstration
% Adjust lambda, L_short, L_long for different frequency/baseline configurations
%
% Outputs:
% - Phase difference vs angle for multiple baselines
% - Ambiguity resolution success rate
% - Angle estimation error vs SNR
% System Parameters
lambda = 0.03; % Wavelength (10 GHz, X-band)
c = 3e8; % Speed of light
% Baseline Configuration (dual-baseline)
L_short = lambda / 2; % Short baseline (unambiguous but coarse)
L_long = 5 * lambda; % Long baseline (ambiguous but fine)
% Angular range
theta_max = 60; % Maximum scan angle (degrees)
theta = linspace(-theta_max, theta_max, 1000);
theta_rad = deg2rad(theta);
% Phase calculation: Delta_phi = 2*pi*L*sin(theta)/lambda
phi_short = (2*pi*L_short*sin(theta_rad)/lambda); % Unambiguous range +/- pi
phi_long = (2*pi*L_long*sin(theta_rad)/lambda); % Highly ambiguous
% Wrap phase to [-pi, pi]
phi_short_wrapped = wrapToPi(phi_short);
phi_long_wrapped = wrapToPi(phi_long);
% Figure 1: Phase Patterns and Ambiguity
figure('Name', 'Phase Interferometer Characteristics', 'Position', [100 100 1400 900]);
subplot(3, 2, 1);
plot(theta, phi_short_wrapped * 180/pi, 'b-', 'LineWidth', 2);
title(sprintf('Short Baseline Phase (L = %.2f\\lambda = %.3fm)', L_short/lambda, L_short));
xlabel('Angle (degrees)');
ylabel('Phase (degrees)');
grid on; ylim([-180 180]);
subplot(3, 2, 2);
plot(theta, phi_long_wrapped * 180/pi, 'r-', 'LineWidth', 1);
title(sprintf('Long Baseline Phase (L = %.2f\\lambda = %.3fm) - Wrapped', L_long/lambda, L_long));
xlabel('Angle (degrees)');
ylabel('Phase (degrees)');
grid on; ylim([-180 180]);
% Unwrapped phase comparison
subplot(3, 2, 3);
plot(theta, phi_short * 180/pi, 'b-', 'LineWidth', 2); hold on;
plot(theta, phi_long * 180/pi, 'r--', 'LineWidth', 1.5);
title('Unwrapped Phase Comparison');
xlabel('Angle (degrees)');
ylabel('Phase (degrees)');
legend(sprintf('Short (L=%.1f\\lambda)', L_short/lambda), sprintf('Long (L=%.1f\\lambda)', L_long/lambda));
grid on;
% Phase slope (sensitivity)
subplot(3, 2, 4);
dtheta = theta(2) - theta(1);
slope_short = gradient(phi_short_wrapped * 180/pi, dtheta);
slope_long = gradient(phi_long_wrapped * 180/pi, dtheta);
plot(theta, abs(slope_short), 'b-', 'LineWidth', 2); hold on;
plot(theta, abs(slope_long), 'r-', 'LineWidth', 2);
title('Phase Sensitivity (deg/degree)');
xlabel('Angle (degrees)');
ylabel('|d\\phi/d\\theta|');
legend('Short Baseline', 'Long Baseline');
grid on;
% Ambiguity resolution demonstration
subplot(3, 2, 5);
% CRT-based resolution: N_long*lambda/L_long - N_short*lambda/L_short = integer ambiguity
% Simplified: use short baseline to resolve long baseline integer ambiguity
N_ambiguities_long = floor(2*L_long/lambda); % Number of ambiguity zones
% Simulate measurement at specific angle
theta_true = 25; % degrees
phi_s_meas = (2*pi*L_short*sind(theta_true)/lambda) + 0.01*randn(1); % Add noise
phi_l_meas = (2*pi*L_long*sind(theta_true)/lambda) + 0.01*randn(1);
% Wrap measured phases
phi_s_w = wrapToPi(phi_s_meas);
phi_l_w = wrapToPi(phi_l_meas);
% Resolve ambiguity
% Estimate coarse angle from short baseline
theta_coarse = asin(phi_s_w * lambda / (2*pi*L_short)) * 180/pi;
% Determine integer ambiguity for long baseline
% k = round((2*L_long/lambda)*(sin(theta_coarse*pi/180) - sin(theta_true*pi/180)))
k_est = round((phi_l_w - (2*pi*L_long/lambda)*sind(theta_coarse)) / (2*pi));
% Corrected long baseline phase
phi_l_corrected = phi_l_w + 2*pi*k_est;
theta_fine = asin(phi_l_corrected * lambda / (2*pi*L_long)) * 180/pi;
% Plot ambiguity zones
theta_range = linspace(-theta_max, theta_max, 1000);
for k = -N_ambiguities_long:N_ambiguities_long
theta_amb = asind((phi_l_w + 2*pi*k) * lambda / (2*pi*L_long));
if ~isreal(theta_amb), continue; end
if abs(theta_amb) <= theta_max
xline(theta_amb, 'g:', 'LineWidth', 1, 'Alpha', 0.5);
hold on;
end
end
xline(theta_true, 'r-', 'LineWidth', 2, 'Label', 'True Angle');
xline(theta_coarse, 'b--', 'LineWidth', 2, 'Label', 'Coarse Estimate');
xline(theta_fine, 'm-', 'LineWidth', 2, 'Label', 'Resolved Angle');
title(sprintf('Ambiguity Resolution at %.1f° (k=%d)', theta_true, k_est));
xlabel('Angle (degrees)');
xlim([-theta_max theta_max]);
grid on;
% Subplot 6: CRT Success Rate vs Phase Noise
subplot(3, 2, 6);
phase_noise_std = linspace(0, 0.5, 50); % Standard deviation in radians
success_rate = zeros(size(phase_noise_std));
N_trials = 1000;
for n_idx = 1:length(phase_noise_std)
successes = 0;
for trial = 1:N_trials
% Random true angle
theta_t = (rand() * 2 * theta_max - theta_max);
% Noisy measurements
phi_s = wrapToPi((2*pi*L_short*sind(theta_t)/lambda) + phase_noise_std(n_idx)*randn());
phi_l = wrapToPi((2*pi*L_long*sind(theta_t)/lambda) + phase_noise_std(n_idx)*randn());
% Resolution attempt
theta_c = asin(phi_s * lambda / (2*pi*L_short)) * 180/pi;
k_guess = round((phi_l - (2*pi*L_long/lambda)*sind(theta_c)) / (2*pi));
phi_l_corr = phi_l + 2*pi*k_guess;
theta_f = asin(phi_l_corr * lambda / (2*pi*L_long)) * 180/pi;
% Check success (within 0.1 degree)
if abs(theta_f - theta_t) < 0.1
successes = successes + 1;
end
end
success_rate(n_idx) = successes / N_trials;
end
plot(phase_noise_std * 180/pi, success_rate * 100, 'b-', 'LineWidth', 2);
title('Ambiguity Resolution Success Rate');
xlabel('Phase Noise Std Dev (degrees)');
ylabel('Success Rate (%)');
grid on; ylim([0 105]);
% Figure 2: Angle Estimation Precision (CRLB and Monte Carlo)
figure('Name', 'Interferometer Precision Analysis', 'Position', [150 150 1200 600]);
% Cramer-Rao Lower Bound for phase interferometer
% var(theta) >= (lambda/(2*pi*L*cos(theta)))^2 / (2*SNR)
SNR_dB = 0:2:40;
SNR_lin = 10.^(SNR_dB/10);
theta_test = 30; % degrees
% Short baseline CRLB
crlb_short = ((lambda/(2*pi*L_short*cosd(theta_test)))^2) ./ (2*SNR_lin);
% Long baseline CRLB (unambiguous phase assumption)
crlb_long = ((lambda/(2*pi*L_long*cosd(theta_test)))^2) ./ (2*SNR_lin);
subplot(1, 2, 1);
semilogy(SNR_dB, sqrt(crlb_short) * 180/pi, 'b-', 'LineWidth', 2); hold on;
semilogy(SNR_dB, sqrt(crlb_long) * 180/pi, 'r-', 'LineWidth', 2);
title(sprintf('Angle Estimation CRLB (\\theta=%.1f°)', theta_test));
xlabel('SNR (dB)');
ylabel('Standard Deviation (degrees)');
legend(sprintf('Short L=%.1f\\lambda', L_short/lambda), sprintf('Long L=%.1f\\lambda', L_long/lambda));
grid on;
% Monte Carlo verification
subplot(1, 2, 2);
snr_test = [10, 20, 30]; % dB
errors_short = zeros(length(snr_test), 1000);
errors_long = zeros(length(snr_test), 1000);
for s_idx = 1:length(snr_test)
snr_lin = 10^(snr_test(s_idx)/10);
sigma_phi = 1/sqrt(2*snr_lin); % Phase std dev
for trial = 1:1000
% True phase
phi_s_true = (2*pi*L_short*sind(theta_test)/lambda);
phi_l_true = (2*pi*L_long*sind(theta_test)/lambda);
% Noisy phase
phi_s_noisy = phi_s_true + sigma_phi*randn();
phi_l_noisy = phi_l_true + sigma_phi*randn();
% Estimate angle
theta_s_est = asin(phi_s_noisy * lambda / (2*pi*L_short)) * 180/pi;
theta_l_est = asin(phi_l_noisy * lambda / (2*pi*L_long)) * 180/pi;
errors_short(s_idx, trial) = theta_s_est - theta_test;
errors_long(s_idx, trial) = theta_l_est - theta_test;
end
end
% Box plot of errors
boxplot_data = [errors_short(1,:)', errors_long(1,:)', ...
errors_short(2,:)', errors_long(2,:)', ...
errors_short(3,:)', errors_long(3,:)'];
boxplot(boxplot_data, 'Labels', {'L_s 10dB', 'L_l 10dB', 'L_s 20dB', 'L_l 20dB', 'L_s 30dB', 'L_l 30dB'});
title('Monte Carlo Error Distribution');
ylabel('Angle Error (degrees)');
grid on;
fprintf('\n=== Phase Interferometer Analysis ===\n');
fprintf('Wavelength: %.3f m (%.1f GHz)\n', lambda, c/lambda/1e9);
fprintf('Short Baseline: %.3f m (%.2f\\lambda) - Unambiguous\n', L_short, L_short/lambda);
fprintf('Long Baseline: %.3f m (%.2f\\lambda) - %d ambiguity zones\n', L_long, L_long/lambda, N_ambiguities_long);
fprintf('Coarse resolution (short): %.2f degrees RMS @ 20dB SNR\n', sqrt(crlb_short(11))*180/pi);
fprintf('Fine resolution (long): %.4f degrees RMS @ 20dB SNR\n', sqrt(crlb_long(11))*180/pi);
fprintf('=====================================\n');
end交付物：parabolic_phase_monopulse.m
脚本功能：基于物理光学近似实现金属隔板分区抛物面天线的和差波束计算，采用镜像法分析孔径场分布，模拟S-200型导引头天线的相位和差特性。
使用方式：设定抛物面焦距、口径直径、隔板位置与偏移量，运行脚本计算和差方向图、差波束零点深度与交叉极化电平。支持隔板位置参数化扫描优化。
matlab
复制
function parabolic_phase_monopulse()
% parabolic_phase_monopulse: Phase Monopulse Feed in Parabolic Reflector
% Models metal septum (partition) configuration for phase comparison monopulse
% Uses Physical Optics (PO) approximation with method of images
% Reference: S-200 type seeker antenna design
%
% Usage:
% Run parabolic_phase_monopulse() for standard S-200 geometry
% Adjust D, f, septum_offset for parametric design studies
%
% Outputs:
% - Sum and Difference patterns from partitioned aperture
% - Aperture field distributions (magnitude and phase)
% - Septum position optimization curves
% Antenna Geometry (S-200 like parameters)
D = 0.6; % Diameter (meters)
f = 0.24; % Focal length (meters) - f/D = 0.4
lambda = 0.03; % X-band wavelength
% Septum (metal partition) parameters
septum_height = 0.05; % Height from focal plane (m)
septum_offset = 0.0; % Lateral offset from axis (m)
septum_tilt = 0; % Tilt angle (degrees)
% Aperture sampling
N_points = 201;
x_aperture = linspace(-D/2, D/2, N_points);
y_aperture = linspace(-D/2, D/2, N_points);
[X, Y] = meshgrid(x_aperture, y_aperture);
% Paraboloid surface z = (x^2 + y^2)/(4f)
R_aperture = sqrt(X.^2 + Y.^2);
valid_aperture = R_aperture <= D/2;
% Feed illumination (simplified Gaussian)
% In phase monopulse, feed is offset or partitioned
theta_feed = atan(R_aperture / f);
E_feed = exp(-2.773 * (theta_feed / (pi/4)).^2); % -3dB at 45 degrees
% Phase center for sum and difference
% Sum: uniform phase across aperture
% Difference: 180 degree phase split by septum
% Septum boundary (vertical partition along x=0 plane with offset)
septum_boundary = X > septum_offset;
% Aperture field for Sum (uniform phase)
E_sum = E_feed .* valid_aperture;
Phase_sum = zeros(size(X));
% Aperture field for Difference (180 deg phase shift across septum)
E_diff = E_feed .* valid_aperture;
Phase_diff = zeros(size(X));
Phase_diff(septum_boundary) = pi; % 180 degree shift
% Complex aperture fields
U_sum = E_sum .* exp(1j * Phase_sum);
U_diff = E_diff .* exp(1j * Phase_diff);
% Far-field calculation (2D FFT for speed, scaled appropriately)
% Physical Optics: E_far ~ FFT(aperture_field)
N_fft = 1024;
U_sum_padded = zeros(N_fft, N_fft);
U_diff_padded = zeros(N_fft, N_fft);
start_idx = floor((N_fft - N_points)/2) + 1;
U_sum_padded(start_idx:start_idx+N_points-1, start_idx:start_idx+N_points-1) = U_sum;
U_diff_padded(start_idx:start_idx+N_points-1, start_idx:start_idx+N_points-1) = U_diff;
% Far field pattern (simplified 2D)
E_far_sum = fftshift(fft2(U_sum_padded));
E_far_diff = fftshift(fft2(U_diff_padded));
% Angular coordinates
theta_max = 30; % degrees
theta_far = linspace(-theta_max, theta_max, N_fft);
% Extract cuts (azimuth plane, y=0)
center_idx = N_fft/2 + 1;
pattern_sum = abs(E_far_sum(center_idx, :)) / max(abs(E_far_sum(:)));
pattern_diff = abs(E_far_diff(center_idx, :)) / max(abs(E_far_sum(:))); % Normalized to sum peak
% Convert to dB
pattern_sum_dB = 20*log10(pattern_sum + eps);
pattern_diff_dB = 20*log10(abs(pattern_diff) + eps);
% Figure 1: Aperture Field Distributions
figure('Name', 'Parabolic Phase Monopulse Analysis', 'Position', [100 100 1400 1000]);
% Aperture amplitude
subplot(3, 2, 1);
imagesc(x_aperture*100, y_aperture*100, E_sum);
axis equal tight;
title('Sum Channel Aperture Amplitude');
xlabel('x (cm)'); ylabel('y (cm)');
colorbar;
colormap('jet');
% Difference amplitude (same)
subplot(3, 2, 2);
imagesc(x_aperture*100, y_aperture*100, E_diff);
axis equal tight;
title('Difference Channel Aperture Amplitude');
xlabel('x (cm)'); ylabel('y (cm)');
colorbar;
% Sum phase
subplot(3, 2, 3);
imagesc(x_aperture*100, y_aperture*100, Phase_sum * 180/pi);
axis equal tight;
title('Sum Channel Aperture Phase');
xlabel('x (cm)'); ylabel('y (cm)');
colorbar; caxis([-180 180]);
% Difference phase (septum partition visible)
subplot(3, 2, 4);
imagesc(x_aperture*100, y_aperture*100, Phase_diff * 180/pi);
axis equal tight;
title('Difference Channel Aperture Phase (Septum Effect)');
xlabel('x (cm)'); ylabel('y (cm)');
colorbar; caxis([-180 180]);
hold on;
contour(X*100, Y*100, double(septum_boundary), [0.5 0.5], 'w--', 'LineWidth', 2);
% Far field patterns
subplot(3, 2, 5);
plot(theta_far, pattern_sum_dB, 'b-', 'LineWidth', 2); hold on;
plot(theta_far, pattern_diff_dB, 'r-', 'LineWidth', 2);
xline(0, 'k:', 'LineWidth', 1);
yline(-3, 'g--', 'LineWidth', 1);
title('Far-Field Radiation Patterns (Azimuth Cut)');
xlabel('Angle (degrees)');
ylabel('Gain (dB)');
legend('Sum (\Sigma)', 'Difference (\Delta)', 'Location', 'southwest');
grid on; ylim([-40 5]);
% Difference pattern detail near null
subplot(3, 2, 6);
idx_null = abs(theta_far) < 5;
plot(theta_far(idx_null), pattern_diff_dB(idx_null), 'r-', 'LineWidth', 2);
title('Difference Pattern Null Detail');
xlabel('Angle (degrees)');
ylabel('Gain (dB)');
grid on;
% Calculate null depth
null_depth = min(pattern_diff_dB(idx_null));
text(0, null_depth+5, sprintf('Null Depth: %.1f dB', null_depth), 'HorizontalAlignment', 'center');
% Figure 2: Septum Position Optimization
figure('Name', 'Septum Position Optimization', 'Position', [150 150 1200 800]);
septum_offsets = linspace(-0.1, 0.1, 21); % +/- 10cm offset
null_depths = zeros(size(septum_offsets));
sum_gains = zeros(size(septum_offsets));
diff_slopes = zeros(size(septum_offsets));
for idx = 1:length(septum_offsets)
offset = septum_offsets(idx);
septum_b = X > offset;
% Recalculate difference field
U_d = E_feed .* valid_aperture .* exp(1j * pi * double(septum_b));
U_d_padded = zeros(N_fft, N_fft);
U_d_padded(start_idx:start_idx+N_points-1, start_idx:start_idx+N_points-1) = U_d;
E_f_d = fftshift(fft2(U_d_padded));
pat_d = abs(E_f_d(center_idx, :)) / max(abs(E_f_d(:)));
pat_d_dB = 20*log10(pat_d + eps);
% Metrics
null_depths(idx) = min(pat_d_dB(abs(theta_far) < 2));
sum_gains(idx) = 20*log10(max(abs(E_far_sum(:))));
% Slope at boresight (approximate from first few points)
slope_idx = find(abs(theta_far) < 1 & theta_far > 0, 5);
if ~isempty(slope_idx)
diff_slopes(idx) = mean(gradient(pat_d_dB(slope_idx), theta_far(2)-theta_far(1)));
end
end
subplot(2, 2, 1);
plot(septum_offsets*100, null_depths, 'b-o', 'LineWidth', 2);
title('Difference Null Depth vs Septum Offset');
xlabel('Septum Offset (cm)');
ylabel('Null Depth (dB)');
grid on;
yline(-30, 'r--', 'Target');
subplot(2, 2, 2);
plot(septum_offsets*100, sum_gains, 'r-s', 'LineWidth', 2);
title('Sum Pattern Gain vs Septum Offset');
xlabel('Septum Offset (cm)');
ylabel('Gain (dBi)');
grid on;
subplot(2, 2, 3);
plot(septum_offsets*100, abs(diff_slopes), 'g-d', 'LineWidth', 2);
title('Difference Slope vs Septum Offset');
xlabel('Septum Offset (cm)');
ylabel('Slope Magnitude (dB/deg)');
grid on;
% Optimal position identification
subplot(2, 2, 4);
% Combined metric: deep null + high slope
metric = -null_depths + abs(diff_slopes)/10; % Weighted combination
plot(septum_offsets*100, metric, 'm-', 'LineWidth', 2);
hold on;
[opt_val, opt_idx] = max(metric);
plot(septum_offsets(opt_idx)*100, opt_val, 'ro', 'MarkerSize', 10, 'MarkerFaceColor', 'r');
title('Optimization Metric (Null Depth + Slope)');
xlabel('Septum Offset (cm)');
ylabel('Combined Metric');
grid on;
text(septum_offsets(opt_idx)*100, opt_val*0.9, sprintf('Opt: %.1f cm', septum_offsets(opt_idx)*100));
% Figure 3: 3D Pattern Visualization
figure('Name', '3D Radiation Patterns', 'Position', [200 200 1000 800]);
% Create 2D angular grid
theta_2d = linspace(-20, 20, 100);
[THETA, PHI] = meshgrid(theta_2d, theta_2d);
% Interpolate patterns for 3D display (simplified)
% Note: Full 3D requires full FFT2 output, here we approximate with cut interpolation
pattern_sum_2d = pattern_sum_dB(center_idx-50:center_idx+50, center_idx-50:center_idx+50);
surf(THETA, PHI, pattern_sum_2d, 'EdgeColor', 'none');
title('Sum Pattern 3D Visualization');
xlabel('Azimuth (deg)'); ylabel('Elevation (deg)'); zlabel('Gain (dB)');
colorbar; colormap('jet'); caxis([-40 0]);
view(45, 30);
fprintf('\n=== Parabolic Phase Monopulse Analysis ===\n');
fprintf('Reflector Diameter: %.2f m\n', D);
fprintf('Focal Length: %.2f m (f/D = %.2f)\n', f, f/D);
fprintf('Septum Offset: %.1f cm\n', septum_offset*100);
fprintf('Difference Null Depth: %.1f dB\n', null_depth);
fprintf('Optimal Septum Offset: %.1f cm\n', septum_offsets(opt_idx)*100);
fprintf('========================================\n');
end交付物：cross_polarization_analysis.m
脚本功能：计算抛物面反射器天线的交叉极化辐射特性，基于几何光学与表面电流积分，分析主极化与交叉极化方向图，评估不同焦距/口径比（f/D）下的交叉极化电平。
使用方式：设定抛物面几何参数、馈源极化类型（线极化或圆极化）与f/D比值范围，运行脚本生成主极化/交叉极化方向图、交叉极化鉴别率（XPD）曲线及45度平面的极化纯度分析。
matlab
复制
function cross_polarization_analysis()
% cross_polarization_analysis: Cross-Polarization Radiation in Offset Parabolic Reflectors
% Calculates co-polar and cross-polar patterns using Physical Optics
% Analyzes XPD (Cross-Polarization Discrimination) vs f/D ratio
%
% Usage:
% Run cross_polarization_analysis() for standard configuration
% Modify f_over_D_range for parametric study of focal length effects
%
% Outputs:
% - Co-polar and cross-polar radiation patterns
% - XPD levels vs angle and vs f/D ratio
% - 45-degree plane polarization purity analysis
% Antenna Parameters
D = 0.5; % Aperture diameter (meters)
lambda = 0.03; % Wavelength (X-band)
f_over_D = 0.4; % Focal length to diameter ratio
f = f_over_D * D; % Focal length
% Feed characteristics (Potter horn approximation)
% E-plane and H-plane patterns slightly different
theta_feed_max = atan(D/(2*f)); % Half-angle subtended by reflector
% Analysis grid
theta_scan = linspace(0, 30, 301); % 0 to 30 degrees
phi_cuts = [0, 45, 90]; % Principal and diagonal planes
% Physical Optics Calculation
% Surface current integration for paraboloid
% Cross-polarization arises from curvature-induced polarization rotation
N_surface = 101;
r_prime = linspace(0, D/2, N_surface);
phi_prime = linspace(0, 2*pi, N_surface);
[R, PHI] = meshgrid(r_prime, phi_prime);
% Paraboloid coordinates
rho = R;
z = rho.^2 / (4*f);
x = rho .* cos(PHI);
y = rho .* sin(PHI);
% Surface normal
dz_drho = rho / (2*f);
n_mag = sqrt(1 + dz_drho.^2);
n_x = -dz_drho ./ n_mag .* cos(PHI);
n_y = -dz_drho ./ n_mag .* sin(PHI);
n_z = 1 ./ n_mag;
% Incident field from feed (assuming x-polarized feed)
% Geometric optics approximation for feed illumination
r_feed = sqrt(x.^2 + y.^2 + (f-z).^2);
theta_feed = atan(rho ./ (f-z));
% Feed pattern (axisymmetric for simplicity, but with spillover)
E_feed = cos(theta_feed).^2; % cos^2(q*theta) model, q chosen for -10dB at edge
E_feed(theta_feed > theta_feed_max) = 0; % Reflector edge truncation
% Incident E-field (x-polarized at feed)
E_inc_x = E_feed;
E_inc_y = zeros(size(E_feed));
E_inc_z = zeros(size(E_feed));
% Surface current calculation (J = 2*n x H_inc)
% For x-polarized incidence, cross-polar component appears in y-direction
% after reflection from curved surface
% Far-field pattern calculation via surface integration
% Simplified: Aperture projection method for cross-polar estimation
% Co-polar (x) and cross-polar (y) aperture fields
% For centered paraboloid with linear polarization:
% Co-polar follows illumination taper
% Cross-polar is zero on-axis, peaks at ~0.7*theta_BW in 45-deg plane
% Aperture projection
x_aperture = x;
y_aperture = y;
U_co = E_feed; % Approximate co-polar illumination
% Cross-polar generation (empirical model based on f/D)
% XPD approx 20*log10(1.14 * f/D) for -10dB edge taper (Ruze formula approximation)
XPD_nominal = 20*log10(1.14 * f_over_D);
% Cross-polar field distribution (peaks at 45-deg plane)
U_cross = E_feed * 10^(-XPD_nominal/20) .* sin(2*PHI); % sin(2*phi) variation
% Far-field via FFT-based integration
N_far = 1024;
dx = D / N_surface;
dy = D / N_surface;
% Pad and transform
U_co_padded = zeros(N_far, N_far);
U_cross_padded = zeros(N_far, N_far);
start = floor((N_far - N_surface)/2) + 1;
U_co_padded(start:start+N_surface-1, start:start+N_surface-1) = U_co;
U_cross_padded(start:start+N_surface-1, start:start+N_surface-1) = U_cross;
E_co_far = fftshift(fft2(U_co_padded));
E_cross_far = fftshift(fft2(U_cross_padded));
% Angular mapping
theta_far = linspace(-30, 30, N_far);
center = N_far/2 + 1;
% Extract cuts
idx_0 = center; % phi=0 plane (E-plane for x-pol)
idx_45 = center + floor(N_far/8); % Approximate phi=45
idx_90 = center + floor(N_far/4); % phi=90 (H-plane)
pattern_co_dB = 20*log10(abs(E_co_far(center, :)) / max(abs(E_co_far(:))) + eps);
pattern_cross_dB = 20*log10(abs(E_cross_far(center, :)) / max(abs(E_co_far(:))) + eps);
% Figure 1: Radiation Patterns
figure('Name', 'Cross-Polarization Analysis', 'Position', [100 100 1400 900]);
% Co-polar pattern
subplot(3, 2, 1);
plot(theta_far, pattern_co_dB, 'b-', 'LineWidth', 2);
title(sprintf('Co-Polar Pattern (f/D = %.2f)', f_over_D));
xlabel('Angle (degrees)');
ylabel('Gain (dBi)');
grid on; ylim([-40 0]);
% Cross-polar pattern
subplot(3, 2, 2);
plot(theta_far, pattern_cross_dB, 'r-', 'LineWidth', 2);
title(sprintf('Cross-Polar Pattern (XPD_{nom} = %.1f dB)', XPD_nominal));
xlabel('Angle (degrees)');
ylabel('Gain (dBi)');
grid on; ylim([-60 -20]);
% Aperture distribution - Co-polar
subplot(3, 2, 3);
imagesc(x_aperture(1,:)*100, y_aperture(:,1)*100, U_co);
axis equal tight;
title('Aperture Field - Co-Polar');
xlabel('x (cm)'); ylabel('y (cm)');
colorbar; colormap('jet');
% Aperture distribution - Cross-polar
subplot(3, 2, 4);
imagesc(x_aperture(1,:)*100, y_aperture(:,1)*100, U_cross);
axis equal tight;
title('Aperture Field - Cross-Polar (sin(2\\phi) variation)');
xlabel('x (cm)'); ylabel('y (cm)');
colorbar;
% XPD vs Angle
subplot(3, 2, 5);
XPD_curve = pattern_co_dB - pattern_cross_dB;
plot(theta_far, XPD_curve, 'g-', 'LineWidth', 2);
title('Cross-Polar Discrimination (XPD) vs Angle');
xlabel('Angle (degrees)');
ylabel('XPD (dB)');
grid on; ylim([0 50]);
yline(40, 'r--', 'Requirement');
% Diagonal plane pattern (45 degrees)
subplot(3, 2, 6);
% Extract 45-degree cut (approximate)
theta_45 = linspace(-30, 30, N_far);
pattern_45_co = 20*log10(abs(diag(E_co_far, floor(N_far/8))) / max(abs(E_co_far(:))) + eps);
pattern_45_cross = 20*log10(abs(diag(E_cross_far, floor(N_far/8))) / max(abs(E_co_far(:))) + eps);
plot(theta_45, pattern_45_co, 'b-', 'LineWidth', 2); hold on;
plot(theta_45, pattern_45_cross, 'r-', 'LineWidth', 2);
title('45-Degree Plane Patterns');
xlabel('Angle (degrees)');
ylabel('Gain (dBi)');
legend('Co-polar', 'Cross-polar');
grid on; ylim([-60 0]);
% Figure 2: f/D Parametric Study
figure('Name', 'f/D Ratio Optimization', 'Position', [150 150 1200 600]);
f_D_range = 0.25:0.05:0.6;
XPD_peak = zeros(size(f_D_range));
XPD_boresight = zeros(size(f_D_range));
for f_idx = 1:length(f_D_range)
f_D = f_D_range(f_idx);
% Empirical formula for cross-polar peak (based on Ludwig and others)
% XPD degrades as f/D decreases
XPD_peak(f_idx) = 20*log10(1.14 * f_D) + 25; % Approximate model
XPD_boresight(f_idx) = 50; % Ideal on-axis (infinity in theory, limited by symmetry)
end
subplot(1, 2, 1);
plot(f_D_range, XPD_peak, 'b-o', 'LineWidth', 2);
title('Peak Cross-Polar Level vs f/D Ratio');
xlabel('f/D');
ylabel('XPD (dB)');
grid on;
yline(40, 'r--', 'Typical Requirement');
text(0.4, 35, 'Target: -40 dB');
% Beamwidth vs f/D
subplot(1, 2, 2);
% HPBW approx 70*lambda/D for large reflectors, slight variation with f/D due to illumination
HPBW = 70 * lambda/D * (1 + 0.01*(0.4 - f_D_range)*10); % Minor variation
plot(f_D_range, HPBW, 'r-s', 'LineWidth', 2);
title('Beamwidth vs f/D Ratio');
xlabel('f/D');
ylabel('HPBW (degrees)');
grid on;
% Figure 3: Polarization Ellipse Visualization
figure('Name', 'Polarization Purity', 'Position', [200 200 800 800]);
% Sample polarization ellipses at different angles
angles_sample = [0, 5, 10, 15];
colors = lines(length(angles_sample));
for i = 1:length(angles_sample)
ang = angles_sample(i);
% Axial ratio (AR) degradation due to cross-pol
% AR = (1 + |E_cross/E_co|) / (1 - |E_cross/E_co|)
E_c = 10^(interp1(theta_far, pattern_co_dB, ang)/20);
E_x = 10^(interp1(theta_far, pattern_cross_dB, ang)/20);
AR = (1 + E_x/E_c) / (1 - E_x/E_c);
% Draw polarization ellipse
subplot(2, 2, i);
t_ell = linspace(0, 2*pi, 100);
a = 1; % Major axis
b = 1/AR; % Minor axis
x_ell = a * cos(t_ell);
y_ell = b * sin(t_ell);
plot(x_ell, y_ell, 'LineWidth', 2, 'Color', colors(i,:));
axis equal; grid on;
title(sprintf('%.0f°: AR = %.2f (XPD=%.1f dB)', ang, AR, 20*log10(E_c/E_x)));
xlabel('E_{co}'); ylabel('E_{cross}');
end
fprintf('\n=== Cross-Polarization Analysis ===\n');
fprintf('Reflector Diameter: %.2f m\n', D);
fprintf('f/D Ratio: %.2f\n', f_over_D);
fprintf('Focal Length: %.2f m\n', f);
fprintf('Estimated Peak XPD: %.1f dB\n', XPD_nominal);
fprintf('Cross-Polar Peak Location: ~%.1f degrees\n', 0.7 * 70*lambda/D);
fprintf('==================================\n');
end交付物：polarization_jamming_sim.m
脚本功能：仿真极化干扰对单脉冲雷达角度跟踪的影响，模拟线极化干扰源与圆极化雷达接收机之间的极化失配效应，计算极化损失因子与角度误差牵引曲线。
使用方式：设定雷达极化类型（圆极化或线极化）、干扰极化角度与极化比（极化纯度），运行脚本生成极化失配损失曲线、和差通道能量互换现象及跟踪失稳判据分析。
matlab
复制
function polarization_jamming_sim()
% polarization_jamming_sim: Polarization Jamming Effects on Monopulse Tracking
% Simulates angle pull-off due to cross-polarization coupling
% Models linear jamming vs circular radar and vice versa
%
% Usage:
% Run polarization_jamming_sim() for standard cross-polar jamming scenario
% Adjust jam_polarization_angle and radar_pol_type for different scenarios
%
% Outputs:
% - Polarization mismatch loss vs jammer aspect angle
% - Sum/Delta channel energy interchange visualization
% - Angle error bias due to polarization squint
% Radar Parameters
radar_pol_type = 'circular'; % 'circular' or 'linear'
radar_pol_tilt = 0; % Tilt angle for linear (degrees)
radar_XPD = 30; % Radar cross-polar isolation (dB)
% Jammer Parameters
jam_power = 1; % Normalized jammer power
jam_polarization_angle = linspace(0, 180, 181); % Jammer E-field angle (degrees)
jam_XPD = 20; % Jammer polarization purity (dB)
% Monopulse Parameters
squint_angle = 0.39; % Beamwidths
beamwidth = 3; % Degrees
% Calculate Stokes Parameters and polarization efficiency
% For circular radar receiving linear jamming
% Radar polarization vectors
if strcmp(radar_pol_type, 'circular')
% Right-hand circular (RHC)
r_radar = [1; -1j] / sqrt(2); % Jones vector for RHC
l_radar = [1; 1j] / sqrt(2); % Left-hand circular (cross-pol)
else
% Linear horizontal
r_radar = [cosd(radar_pol_tilt); sind(radar_pol_tilt)];
l_radar = [-sind(radar_pol_tilt); cosd(radar_pol_tilt)]; % Orthogonal
end
% Jammer polarization (linear, rotating)
jam_pol = [cosd(jam_polarization_angle); sind(jam_polarization_angle)];
% Polarization mismatch loss (channel transmission coefficients)
% Co-pol channel: |r_radar' * jam_pol|^2
% Cross-pol channel: |l_radar' * jam_pol|^2
co_pol_response = abs(r_radar' * jam_pol).^2;
cross_pol_response = abs(l_radar' * jam_pol).^2;
% Include XPD limitations (radar doesn't have perfect isolation)
% Actual received power in co-pol and cross-pol channels
co_pol_rx = co_pol_response + (10^(-radar_XPD/10)) * cross_pol_response;
cross_pol_rx = cross_pol_response + (10^(-radar_XPD/10)) * co_pol_response;
% Normalized mismatch loss
mismatch_loss_dB = 10*log10(co_pol_rx + eps);
% Monopulse processing with polarization contamination
% Normal sum and delta patterns
theta = linspace(-beamwidth, beamwidth, 100); % Scan angle
% Ideal patterns (Gaussian)
E_left = exp(-2.773 * ((theta/beamwidth + squint_angle)).^2);
E_right = exp(-2.773 * ((theta/beamwidth - squint_angle)).^2);
Sigma_ideal = E_left + E_right;
Delta_ideal = E_left - E_right;
% With polarization contamination (cross-pol pattern has different phase/squint)
% Cross-pol pattern typically has null on-axis and peaks off-axis
Sigma_cross = 0.1 * (1 - exp(-2.773 * (theta/(0.8*beamwidth)).^2)); % Dished shape
% Angle error calculation for different jammer polarizations
angle_error_bias = zeros(size(jam_polarization_angle));
for j = 1:length(jam_polarization_angle)
% Weighted sum of co-pol and cross-pol responses
alpha = co_pol_rx(j); % Co-pol weight
beta = cross_pol_rx(j); % Cross-pol weight (contamination)
% Composite patterns
Sigma_total = alpha * Sigma_ideal + beta * Sigma_cross;
% Delta pattern distorted by cross-pol (different spatial distribution)
Delta_distorted = alpha * Delta_ideal + 0.5 * beta * Sigma_ideal .* sign(theta);
% Monopulse error voltage (normalized)
error_voltage = real(Delta_distorted ./ Sigma_total);
% Find zero crossing (tracking point)
% Interpolate to find where error = 0
zero_crossings = find(diff(sign(error_voltage)) ~= 0);
if ~isempty(zero_crossings)
% Linear interpolation for precise zero
idx = zero_crossings(1);
theta_z = theta(idx) - error_voltage(idx) * (theta(idx+1) - theta(idx)) / (error_voltage(idx+1) - error_voltage(idx));
angle_error_bias(j) = theta_z; % Bias from true boresight
else
angle_error_bias(j) = NaN; % No tracking solution (break lock)
end
end
% Figure 1: Polarization Mismatch Loss
figure('Name', 'Polarization Jamming Analysis', 'Position', [100 100 1400 900]);
subplot(3, 2, 1);
plot(jam_polarization_angle, mismatch_loss_dB, 'b-', 'LineWidth', 2);
title(sprintf('Polarization Mismatch Loss (%s radar)', radar_pol_type));
xlabel('Jammer Polarization Angle (degrees)');
ylabel('Received Power (dB)');
grid on;
yline(-3, 'r--', '3dB Loss');
% Co-pol and Cross-pol responses
subplot(3, 2, 2);
plot(jam_polarization_angle, 10*log10(co_pol_rx+eps), 'g-', 'LineWidth', 2); hold on;
plot(jam_polarization_angle, 10*log10(cross_pol_rx+eps), 'r-', 'LineWidth', 2);
title('Channel Responses');
xlabel('Jammer Polarization Angle (degrees)');
ylabel('Response (dB)');
legend('Co-Pol Channel', 'Cross-Pol Channel');
grid on;
% Angle error bias
subplot(3, 2, 3);
plot(jam_polarization_angle, angle_error_bias, 'm-', 'LineWidth', 2);
title('Tracking Angle Bias due to Polarization');
xlabel('Jammer Polarization Angle (degrees)');
ylabel('Angle Error (degrees)');
grid on;
yline(0, 'k:', 'LineWidth', 1);
yline(beamwidth/2, 'r--', 'Half Beamwidth');
% Identify critical regions (error > 0.3 beamwidths = break lock)
break_lock = abs(angle_error_bias) > 0.3 * beamwidth;
hold on;
plot(jam_polarization_angle(break_lock), angle_error_bias(break_lock), 'ro', 'MarkerSize', 4);
title('Tracking Angle Bias (Red dots indicate break-lock)');
% Sum/Delta interchange visualization (specific angles)
subplot(3, 2, 4);
critical_angles = [0, 45, 90, 135];
colors = lines(length(critical_angles));
for c_idx = 1:length(critical_angles)
ang = critical_angles(c_idx);
alpha = co_pol_rx(jam_polarization_angle == ang);
beta = cross_pol_rx(jam_polarization_angle == ang);
Sigma_comp = alpha * Sigma_ideal + beta * Sigma_cross;
Delta_comp = alpha * Delta_ideal + 0.5 * beta * Sigma_ideal .* sign(theta);
plot(theta, Sigma_comp/max(Sigma_comp), 'LineWidth', 1.5, 'Color', colors(c_idx,:), 'LineStyle', '-'); hold on;
plot(theta, Delta_comp/max(Delta_comp), 'LineWidth', 1.5, 'Color', colors(c_idx,:), 'LineStyle', '--');
end
title('Sum/Delta Patterns at Critical Polarizations');
xlabel('Angle (degrees)');
ylabel('Normalized Amplitude');
legend(arrayfun(@(x) sprintf('%d° Sum', x), critical_angles, 'UniformOutput', false), 'Location', 'best');
grid on;
% 3D Polarization sphere visualization (Poincare Sphere projection)
subplot(3, 2, 5);
% Stokes parameters for visualization
S1 = cosd(2*jam_polarization_angle);
S2 = sind(2*jam_polarization_angle);
S3 = zeros(size(jam_polarization_angle)); % Linear only
plot(S1, S2, 'b-', 'LineWidth', 2);
axis equal; grid on;
title('Poincare Sphere Equatorial Projection');
xlabel('S1 (Horizontal/Vertical)');
ylabel('S2 (45/-45 deg)');
hold on;
plot(cosd(2*radar_pol_tilt), sind(2*radar_pol_tilt), 'ro', 'MarkerSize', 10, 'MarkerFaceColor', 'r');
text(cosd(2*radar_pol_tilt)+0.1, sind(2*radar_pol_tilt), 'Radar Pol');
% Power exchange diagram
subplot(3, 2, 6);
% Sankey-style diagram showing power flow
bar_data = [max(co_pol_rx), max(cross_pol_rx), max(co_pol_rx)*10^(-radar_XPD/10), max(cross_pol_rx)*10^(-radar_XPD/10)];
bar(bar_data);
set(gca, 'XTickLabel', {'Jam->Co', 'Jam->Cross', 'Leak->Co', 'Leak->Cross'});
title('Power Flow in Polarization Channels');
ylabel('Relative Power');
% Figure 2: Jamming Effectiveness Map
figure('Name', 'Jamming Effectiveness', 'Position', [150 150 1000 600]);
% 2D map: Jammer angle vs Jammer power (showing break-lock regions)
jam_angles = 0:5:180;
jam_powers = -20:2:20; % dB relative
[J_ANG, J_PWR] = meshgrid(jam_angles, jam_powers);
effectiveness = zeros(size(J_ANG));
for i = 1:size(J_ANG, 1)
for j = 1:size(J_ANG, 2)
% Calculate if break-lock occurs
idx = find(abs(jam_polarization_angle - J_ANG(i,j)) < 2.5, 1);
if ~isempty(idx) && ~isnan(angle_error_bias(idx))
if abs(angle_error_bias(idx)) > 0.3 * beamwidth
effectiveness(i,j) = 1; % Break lock
elseif abs(angle_error_bias(idx)) > 0.1 * beamwidth
effectiveness(i,j) = 0.5; % Significant bias
end
end
end
end
imagesc(jam_angles, jam_powers, effectiveness);
set(gca, 'YDir', 'normal');
colormap([0.8 0.8 0.8; 0.9 0.6 0.2; 0.8 0.2 0.2]);
colorbar('Ticks', [0.25, 0.75], 'TickLabels', {'Minor', 'Break Lock'});
title('Jamming Effectiveness Map (Angle Pull-Off)');
xlabel('Jammer Polarization Angle (degrees)');
ylabel('Jammer Power (dB)');
hold on;
contour(J_ANG, J_PWR, effectiveness, [0.5 0.5], 'w--', 'LineWidth', 2);
fprintf('\n=== Polarization Jamming Analysis ===\n');
fprintf('Radar Type: %s\n', radar_pol_type);
fprintf('Radar XPD: %d dB\n', radar_XPD);
fprintf('Maximum Angle Bias: %.2f degrees\n', max(abs(angle_error_bias), [], 'omitnan'));
fprintf('Break-lock threshold: %.1f beamwidths\n', 0.3);
fprintf('Critical polarization angles: ');
crit_idx = find(abs(angle_error_bias) > 0.3*beamwidth);
if ~isempty(crit_idx)
fprintf('%.0f°, ', jam_polarization_angle(crit_idx(1:min(5,end))));
end
fprintf('\n=====================================\n');
end
255
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
