%% 主脚本:相位编码信号分析
clear; clc; close all;
%% 参数设置(保持不变)
fc=10e9;
% lightspeed=299792458;
lightspeed=3e8;
lambda=lightspeed/fc;
pd = 0.9;
pfa = 1e-6;
trcs = 1.5;
%% 创建阵列(保持不变)
Array = phased.URA('Size',[38 38],...
'Lattice','Rectangular','ArrayNormal','x');
Array.ElementSpacing = [0.5 0.5]*lambda;
rwind = hamming(38);
cwind = hamming(38);
taper = rwind*cwind.';
Array.Taper = taper.';
arraygain = phased.ArrayGain('SensorArray',Array,'PropagationSpeed',physconst('LightSpeed'));
ag = arraygain(fc,[0;0]);% 天线法向增益
Elem = phased.IsotropicAntennaElement;
Elem.BackBaffled = true;
Elem.FrequencyRange = [9.99e9 10.01e9];
Array.Element = Elem;
figure;
pattern(Array.Element,fc,'PropagationSpeed',physconst('LightSpeed'),...
'Type','powerdb');% 单个天线基本阵子方向图
figure;
pattern(Array,fc,'PropagationSpeed',physconst('LightSpeed'),...
'Type','powerdb');% 阵列天线方向图
%% 获取方向图数据
az = -60:0.01:60; % 方位角范围
el = 0; % 固定俯仰角
P = pattern(Array, fc, az, el, 'Type', 'powerdb', ...
'PropagationSpeed', physconst('LightSpeed'));
P = P(:); % 确保为列向量
%% 核心:寻找所有峰值(主瓣+旁瓣)
% 1. 找到主瓣峰值
[main_peak_val, main_peak_idx] = max(P);
main_peak_az = az(main_peak_idx);
% 2. 确定主瓣区域(根据3dB波束宽度)
half_power = main_peak_val - 3;
left_idx = find(P(1:main_peak_idx) <= half_power, 1, 'last');
right_idx = main_peak_idx + find(P(main_peak_idx:end) <= half_power, 1, 'first') - 1;
beamwidth = az(right_idx) - az(left_idx);
mainlobe_start_idx2 = max(1, main_peak_idx - round(beamwidth/0.01/2));
mainlobe_end_idx2 = min(length(az), main_peak_idx + round(beamwidth/0.01/2));
% 3. 排除主瓣区域
mainlobe_start_idx = max(1, main_peak_idx - round(3.4*beamwidth/0.01/2));
mainlobe_end_idx = min(length(az), main_peak_idx + round(3.4*beamwidth/0.01/2));
% 4. 标记非主瓣区域
sidelobe_mask = true(size(P));
sidelobe_mask(mainlobe_start_idx:mainlobe_end_idx) = false;
% 5. 在非主瓣区域寻找所有旁瓣峰值
sidelobe_data = P(sidelobe_mask);
sidelobe_az = az(sidelobe_mask);
% 使用findpeaks函数寻找旁瓣峰值
[peak_values, peak_indices] = findpeaks(sidelobe_data, ...
'MinPeakDistance', 50,... % 至少50个点(0.5度)的间隔
'MinPeakHeight', main_peak_val - 60); % 只找高于主瓣-60dB的峰值
% 转换回原始索引
sidelobe_peak_indices = find(sidelobe_mask);
sidelobe_peak_indices = sidelobe_peak_indices(peak_indices);
sidelobe_peak_az = az(sidelobe_peak_indices);
% 6. 确定关键旁瓣参数
if ~isempty(peak_values)
% 最高旁瓣(SLL)
[max_sll_val, max_sll_idx] = max(peak_values);
max_sll_az = sidelobe_peak_az(max_sll_idx);
max_sll_db = max_sll_val - main_peak_val; % 相对于主瓣的dB值
% 第一旁瓣(最靠近主瓣的旁瓣)
[~, first_sll_idx] = min(abs(sidelobe_peak_az - main_peak_az));
first_sll_val = peak_values(first_sll_idx);
first_sll_az = sidelobe_peak_az(first_sll_idx);
first_sll_db = first_sll_val - main_peak_val;
else
warning('未找到有效旁瓣峰值,请调整角度范围或参数');
max_sll_db = -inf;
first_sll_db = -inf;
end
%% 绘图显示阵列方向图结果
figure;
plot(az, P, 'b', 'LineWidth', 1.2);
hold on;
% 标记主瓣峰值
plot(main_peak_az, main_peak_val, 'ro', 'MarkerSize', 8, 'MarkerFaceColor', 'r');
text(main_peak_az+1, main_peak_val, '主瓣峰值', 'Color', 'r');
% 标记主瓣区域
fill([az(mainlobe_start_idx2) az(mainlobe_end_idx2) az(mainlobe_end_idx2) az(mainlobe_start_idx2)], ...
[min(P) min(P) max(P) max(P)], 'y', 'FaceAlpha', 0.2);
text(az(main_peak_idx), min(P)+5, '主瓣区域', 'HorizontalAlignment', 'center');
% 标记所有旁瓣峰值
if ~isempty(peak_values)
plot(sidelobe_peak_az, peak_values, 'go', 'MarkerSize', 6, 'MarkerFaceColor', 'g');
% 标记最高旁瓣
plot(max_sll_az, max_sll_val, 'mo', 'MarkerSize', 8, 'MarkerFaceColor', 'm');
text(max_sll_az+1, max_sll_val, sprintf('最高旁瓣: %.1fdB', max_sll_db), 'Color', 'm');
end
grid on;
xlabel('方位角(度)');
ylabel('功率(dB)');
title('阵列方向图及旁瓣峰值标记');
xlim([min(az) max(az)]);
ylim([main_peak_val-70 main_peak_val+5]);
%% 输出阵列参数结果
fprintf('主瓣峰值: %.2f dB\n', main_peak_val);
fprintf('波束宽度: %.2f 度\n', beamwidth);
fprintf('最高旁瓣电平(SLL): %.2f dBc\n', max_sll_db); % dBc表示相对于载波(主瓣)
if ~isempty(peak_values)
fprintf('第一旁瓣电平: %.2f dBc\n', first_sll_db);
end
%% 设计线性调频信号
Waveform = phased.LinearFMWaveform('SampleRate',10e+06,...
'PulseWidth',5e-6,'PRF',5000,'SweepBandwidth',10e6,...
'SweepDirection','Up','SweepInterval','Positive'...
,'Envelope','Rectangular','NumPulses',1,'FrequencyOffset',10e9);
Fs = Waveform.SampleRate;
prf = Waveform.PRF;
NumOfPulses = 6; % 用于多普勒分析的脉冲数
tb_product = Waveform.SweepBandwidth * Waveform.PulseWidth;% 时宽带宽积
delta_R = lightspeed / (2*Waveform.SweepBandwidth);% 距离分辨率
delta_v = (lambda * prf) / (2*NumOfPulses);% 速度分辨率
compression_gain = 10*log10(tb_product); % 脉冲压缩增益(dB)
snr_min_raw = albersheim(pd, pfa, NumOfPulses); % 未考虑压缩增益的SNR
snr_min_actual = snr_min_raw - compression_gain; % 实际所需SNR
% 生成基带信号
x = Waveform();
% 匹配滤波
coeff = getMatchedFilter(Waveform);
Compression = phased.MatchedFilter('Coefficients',coeff(:,1), ...
'SpectrumWindow','None');
y = Compression(x);
%% 设计发射机和接收机
transmitter=phased.Transmitter('PeakPower',1,...
'Gain',20,'InUseOutputPort',true,...
'CoherentOnTransmit',true);
receiver=phased.ReceiverPreamp(...
"Gain",20,'NoiseFigure',10.0,...
'SampleRate',Fs,'EnableInputPort',true);
% 噪声参数计算
F_dB = 7; % 噪声系数
T0 = 290; % 参考温度 (K)
L_dB = 10; % 系统损失
k = 1.38e-23; % 玻尔兹曼常数
max_range = 20e3; % 最大探测距离
P_noise = k * T0 * Waveform.SweepBandwidth * db2pow(F_dB); % 噪声功率
% 峰值功率和最大作用距离计算
peak_power = ((4*pi)^3*k*T0*Waveform.SweepBandwidth*db2pow(F_dB)*max_range^4*...
db2pow(snr_min_raw))/(db2pow(2*(ag+transmitter.Gain))*trcs*lambda^2);
Rmax=((transmitter.PeakPower * db2pow(ag+transmitter.Gain) * db2pow(ag+receiver.Gain)...
* lambda^2 * trcs) ./ (P_noise ...
* (4 * pi)^3 * db2pow(snr_min_raw))).^(1/4);
transmitter.PeakPower = peak_power/1444;
%% 目标和传感器运动设置
initialAz = 60;
endAz = -60;
scanstep = -4; % 扫描步进角度
scangrid = initialAz+scanstep/2:scanstep:endAz;
numscans = length(scangrid);
pulsenum = NumOfPulses*numscans;
% 目标初始位置和速度
tgtpos = [[19900*cosd(0); 19900*sind(0); 0],[19950*cosd(2); 19950*sind(2); 0],[20000*cosd(4); 20000*sind(4); 0]];
tgtvel = [[10; 0; 0],[20; 0; 0],[30; 0; 0]];% 目标真实径向速度
% tgtpos = [19900*cosd(45); 19900*sind(45); 0];
% tgtvel = [-15; 0; 0];
tgtmotion = phased.Platform('InitialPosition',tgtpos,'Velocity',tgtvel);
tgtrcs = [1.5 1.5 1.5];% 目标RCS
% tgtrcs = 1.5;
target = phased.RadarTarget('MeanRCS',tgtrcs,'OperatingFrequency',fc);
% 传感器位置(静止)
sensorpos = [0; 0; 0];
sensorvel = [0; 0; 0];
sensormotion = phased.Platform('InitialPosition',sensorpos,'Velocity',sensorvel);
% 初始目标距离和角度
[tgtrng,tgtang] = rangeangle(tgtmotion.InitialPosition,...
sensormotion.InitialPosition);
numtargets = length(target.MeanRCS);
% 导向矢量和波束形成器
steeringvec = phased.SteeringVector('SensorArray',Array,...
'PropagationSpeed',lightspeed);
beamformer = phased.PhaseShiftBeamformer('SensorArray',Array,...
'OperatingFrequency',fc,'PropagationSpeed',lightspeed,...
'DirectionSource','Input port');
% 传播信道
channel = phased.FreeSpace(...
'SampleRate',Fs,...
'TwoWayPropagation',true,...
'OperatingFrequency',fc);
% 快时间网格和距离轴
fast_time_grid = unigrid(0, 1/Fs, 1/prf, '[)');
range_axis = lightspeed * fast_time_grid / 2; % 距离轴
% 接收脉冲存储
waveform_samples = length(Waveform());
rxpulses = zeros(waveform_samples, pulsenum, 'double');
% 辐射器和收集器
radiator = phased.Radiator(...
'Sensor',Array,...
'PropagationSpeed', lightspeed,...
'OperatingFrequency', fc,...
'CombineRadiatedSignals', true,...
'WeightsInputPort',true);
collector = phased.Collector(...
'Sensor',Array,...
'PropagationSpeed', lightspeed,...
'OperatingFrequency', fc,...
'Wavefront', 'Plane',...
'WeightsInputPort',false);
%% 回波模拟
for m = 1:pulsenum
fprintf('脉冲数: %d \n', m);
% 更新位置
[sensorpos,sensorvel] = sensormotion(1/prf);
[tgtpos,tgtvel] = tgtmotion(1/prf);
% 目标角度计算
[tgtrng,tgtang] = rangeangle(tgtpos,sensorpos);
% 扫描角度和导向矢量
scanid = floor((m-1)/NumOfPulses) + 1;
sv = steeringvec(fc,scangrid(scanid));
% w = conj(sv);
w = conj(sv).*taper(:);
% 发射信号和传播
pulse = Waveform();
% fprintf('波形生成功率: %.2f \n', sum(abs(pulse).^2)/400);
[txsig,txstatus] = transmitter(pulse);
% fprintf('发射波形功率: %.2f \n', sum(abs(txsig).^2)/400);
txsig = radiator(txsig,tgtang,w);
% fprintf('发射天线功率: %.2f \n', sum(abs(txsig).^2)/400);
% fprintf('传输路径功率理论计算值: %.30f \n',sum(abs(txsig).^2)/400/((4*pi*19900/lambda).^4));
% fprintf('传输路径幅度理论计算值: %.30f \n',sum(abs(txsig))/400/((4*pi*19900/lambda).^2));
% a0=sum(abs(txsig).^2)/400/((4*pi*19900/lambda).^4);
% a1=sum(abs(txsig))/400/((4*pi*19900/lambda).^2);
txsig = channel(txsig,sensorpos,tgtpos,sensorvel,tgtvel);
% fprintf('传输路径功率实验仿真值: %.30f \n', sum(abs(txsig).^2)/400);
% fprintf('传输路径幅度实验仿真值: %.30f \n',sum(abs(txsig))/400);
% b0=sum(abs(txsig).^2)/400;
% b1=sum(abs(txsig))/400;
% fprintf('功率衰减系数:%.30f \n',a0/b0);
% fprintf('幅度衰减系数:%.30f \n',a1/b1);
% 目标反射和接收处理
tgtsig = target(txsig);
% fprintf('目标反射功率: %.30f \n', sum(abs(tgtsig).^2)/400);
rxsig = collector(tgtsig,tgtang);
% fprintf('接收天线功率: %.30f \n', sum(sum(abs(rxsig).^2))/400);
rxsig = receiver(rxsig,~(txstatus>0));
% fprintf('接收波形功率: %.30f \n', sum(sum(abs(rxsig).^2))/400);
rxpulses(:,m) = beamformer(rxsig,[scangrid(scanid);0]);
% fprintf('接收平均功率: %.30f \n', sum(sum(abs(rxpulses).^2))/180);
end
matchingcoeff = getMatchedFilter(Waveform);
matchedfilter = phased.MatchedFilter(...
'Coefficients',matchingcoeff,...
'GainOutputPort',true);
compressed= matchedfilter(rxpulses);
% compressed= reshape(compressed,[],NumOfPulses,numscans);
% % 脉冲压缩(提升距离分辨率)
% compressed = zeros(size(rxpulses));
% for m = 1:pulsenum
% compressed(:,m) = Compression(rxpulses(:,m));
% end
% 匹配滤波后,保存STC前的数据
compressed_before_STC = compressed;
% 应用STC
range_gates = range_axis; % 获取距离门向量
tvg = phased.TimeVaryingGain('RangeLoss', 2*fspl(range_gates, lambda));
compressed = tvg(compressed); % 应用STC增益补偿
% 绘制STC前后的距离-幅度曲线
figure;
subplot(2,1,1);
plot(range_axis, 20*log10(max(abs(compressed_before_STC), [], 2) + 1e-10));
title('STC前:距离-幅度曲线');
xlabel('距离(m)'); ylabel('幅度(dB)');
grid on;
subplot(2,1,2);
plot(range_axis, 20*log10(max(abs(compressed), [], 2) + 1e-10));
title('STC后:距离-幅度曲线');
xlabel('距离(m)'); ylabel('幅度(dB)');
grid on;
% 转换为分贝幅度
magnitude_dB = 20*log10(abs(compressed) + 1e-10); % +1e-10避免log(0)
% 绘制热力图
figure;
imagesc(1:pulsenum, range_axis, magnitude_dB);
axis xy; % 翻转y轴,使距离从下到上递增
xlabel('脉冲编号(慢时间)');
ylabel('距离(米)');
title('rxpulses的距离-慢时间热力图');
colorbar; % 颜色对应幅度(dB)
colormap('jet'); % 红→黄→蓝表示强度从高到低
% 1. 确保已完成匹配滤波(使用前序步骤的compressed_pulses)
% 2. 准备坐标轴数据
[X, Y] = meshgrid(1:pulsenum, range_axis); % X:脉冲编号(矩阵),Y:距离(矩阵)
Z = 20*log10(abs(compressed) + 1e-12); % Z:幅度(dB),与X、Y维度匹配
% 3. 绘制三维网格图(兼顾细节与可读性)
figure('Position', [100, 100, 1000, 800]);
mesh(X, Y, Z, 'EdgeColor', 'interp', 'FaceColor', 'none'); % 线框插值着色,无填充
colormap('jet'); % 红→黄表示高幅度,蓝→青表示低幅度
colorbar;
caxis([max(Z(:))-60, max(Z(:))]); % 限制动态范围,突出目标
% 4. 坐标轴与标题设置
xlabel('脉冲编号(慢时间)', 'FontSize', 10);
ylabel('距离(米)', 'FontSize', 10);
zlabel('幅度(dB)', 'FontSize', 10);
title('匹配滤波后rxpulses的三维网格图(距离-脉冲编号-幅度)', 'FontSize', 12);
view(30, 45); % 视角:方位角30°,仰角45°(可调整)
grid on;
axis tight;
%% 多目标距离检测(基于距离轴峰值)
% 1. 计算距离轴上的幅度(跨所有脉冲的最大值)
max_amp_per_range = max(abs(compressed), [], 2); % 按距离维度取最大值(维度:[快时间,1])
% 2. 设置检测阈值(基于噪声电平,此处用最大幅度的30%作为阈值)
noise_floor = prctile(max_amp_per_range, 10); % 估计噪声底(10%分位值)
detection_threshold = max(3*noise_floor, 0.3*max(max_amp_per_range)); % 检测阈值
% 3. 找到所有距离峰值(多目标)
[range_peak_amps, range_peak_idxs] = findpeaks(max_amp_per_range, ...
'MinPeakHeight', detection_threshold, ... % 仅保留超过阈值的峰值
'MinPeakDistance', 5); % 最小距离间隔(大于距离分辨率15m,此处取5个距离点)
% 4. 按距离从小到大排序(便于对应目标编号)
[~, sorted_idx] = sort(range_peak_idxs);
range_peak_idxs = range_peak_idxs(sorted_idx); % 排序后的距离索引
measured_ranges = range_axis(range_peak_idxs); % 测量的目标距离(m)
num_targets = length(measured_ranges); % 检测到的目标数量
%% 1. 构建角度-距离-多普勒三维谱
% 1.1 按扫描角度拆分数据(每个角度对应NumOfPulses个脉冲)
angle_data = cell(numscans, 1); % 存储每个角度的回波数据(快时间×慢时间)
for scan_idx = 1:numscans
pulse_start = (scan_idx - 1)*NumOfPulses + 1;
pulse_end = scan_idx*NumOfPulses;
angle_data{scan_idx} = compressed(:, pulse_start:pulse_end); % 维度:[快时间, NumOfPulses]
end
% 1.2 计算每个角度的多普勒谱(快时间×多普勒频率)
doppler_freq = (-NumOfPulses/2 : NumOfPulses/2 - 1) * (prf / NumOfPulses); % 多普勒频率轴(Hz)
angle_range_doppler = zeros(size(compressed, 1), length(doppler_freq), numscans); % 三维谱:[距离, 多普勒, 角度]
for scan_idx = 1:numscans
% 对当前角度的慢时间数据做FFT(获取多普勒信息)
doppler_fft = fftshift(fft(angle_data{scan_idx}, [], 2), 2); % 慢时间FFT并移位
angle_range_doppler(:, :, scan_idx) = abs(doppler_fft); % 存储幅度
end
% 转换为dB(便于可视化和检测)
angle_range_doppler_dB = 20*log10(angle_range_doppler + 1e-10);
%% 2. 三维峰值检测(区分多目标)
% 2.1 设置检测阈值(基于噪声估计)
noise_floor = prctile(angle_range_doppler_dB(:), 5); % 5%分位值估计噪声底
detection_threshold = noise_floor + 15; % 高于噪声15dB作为阈值
% 2.2 三维峰值搜索(使用imregionalmax找局部最大值)
peak_mask = imregionalmax(angle_range_doppler_dB, 8); % 8邻域局部最大值
peak_coords = find(peak_mask); % 峰值线性索引
[range_idxs, doppler_idxs, angle_idxs] = ind2sub(size(peak_mask), peak_coords); % 转换为三维坐标
% 2.3 筛选超过阈值的峰值,并按能量排序
peak_values = angle_range_doppler_dB(peak_mask);
valid_idx = peak_values >= detection_threshold;
range_idxs = range_idxs(valid_idx);
doppler_idxs = doppler_idxs(valid_idx);
angle_idxs = angle_idxs(valid_idx);
peak_values = peak_values(valid_idx);
% 按能量从高到低排序(确保目标被优先检测)
[~, sorted_idx] = sort(peak_values, 'descend');
range_idxs = range_idxs(sorted_idx);
doppler_idxs = doppler_idxs(sorted_idx);
angle_idxs = angle_idxs(sorted_idx);
% 2.4 保留有效目标(假设已知目标数量为3)
num_targets = 3;
if length(range_idxs) < num_targets
warning('检测到的目标数量不足%d个', num_targets);
num_targets = length(range_idxs);
end
range_idxs = range_idxs(1:num_targets);
doppler_idxs = doppler_idxs(1:num_targets);
angle_idxs = angle_idxs(1:num_targets);
%% 3. 计算每个目标的参数
% 存储测量值与真实值
measured = struct('range', [], 'az', [], 'velocity', []);
true_params = struct('range', [19900, 19950, 20000], ...
'az', [0, 2, 4], ...
'velocity', [10, 20, 30]);
for t = 1:num_targets
% 3.1 距离计算(对应range_axis)
measured.range(t) = range_axis(range_idxs(t));
% 3.2 角度计算(对应扫描角度scangrid)
measured.az(t) = scangrid(angle_idxs(t));
% 3.3 速度计算(多普勒频率→速度)
fd = doppler_freq(doppler_idxs(t)); % 多普勒频率
measured.velocity(t) = (fd * lambda) / 2; % 径向速度公式
end
% 按距离排序(匹配目标编号1-3)
[~, sort_idx] = sort(measured.range);
measured.range = measured.range(sort_idx);
measured.az = measured.az(sort_idx);
measured.velocity = measured.velocity(sort_idx);
%% 4. 结果输出
fprintf('\n===== 三目标(角度不同)测量结果 =====\n');
fprintf('目标编号 | 真实距离(m) | 测量距离(m) | 距离误差(m) | 真实角度(°) | 测量角度(°) | 角度误差(°) | 真实速度(m/s) | 测量速度(m/s) | 速度误差(m/s)\n');
fprintf('----------------------------------------------------------------------------------------------------------------------------------------\n');
for t = 1:num_targets
range_err = abs(measured.range(t) - true_params.range(t));
az_err = abs(measured.az(t) - true_params.az(t));
vel_err = abs(measured.velocity(t) - true_params.velocity(t));
fprintf('%-8d | %-12.1f | %-12.1f | %-11.2f | %-12.1f | %-12.1f | %-11.2f | %-14.1f | %-14.1f | %-14.2f\n', ...
t, true_params.range(t), measured.range(t), range_err, ...
true_params.az(t), measured.az(t), az_err, ...
true_params.velocity(t), measured.velocity(t), vel_err);
end
%% 5. 可视化:角度-距离谱(直观展示目标分离)
figure;
xx=max(angle_range_doppler_dB);
imagesc(scangrid, range_axis, max(angle_range_doppler_dB, [], 2)); % 对多普勒维度取最大值
axis xy; colorbar; colormap('jet');
xlabel('方位角(°)'); ylabel('距离(m)');
title('角度-距离谱(多目标分离)');
hold on;
% 标记检测到的目标
plot(measured.az, measured.range, 'ro', 'MarkerSize', 8, 'MarkerFaceColor', 'r', 'DisplayName', '测量目标');
plot(true_params.az, true_params.range, 'go', 'MarkerSize', 8, 'MarkerFaceColor', 'g', 'DisplayName', '真实目标');
legend();
%% 6. 可视化:距离-多普勒谱(验证速度区分)
% 提取目标角度对应的距离-多普勒谱
figure;
for t = 1:num_targets
subplot(num_targets, 1, t);
scan_idx = find(scangrid == measured.az(t), 1); % 目标角度对应的扫描索引
plot(doppler_freq, angle_range_doppler_dB(range_idxs(t), :, scan_idx), 'b-');
hold on;
plot(doppler_freq(doppler_idxs(t)), angle_range_doppler_dB(range_idxs(t), doppler_idxs(t), scan_idx), ...
'ro', 'MarkerSize', 6, 'MarkerFaceColor', 'r');
xlabel('多普勒频率(Hz)'); ylabel('幅度(dB)');
title(sprintf('目标%d:距离=%.0fm,角度=%.1f°', t, measured.range(t), measured.az(t)));
grid on;
end
把以下cfar阈值计算加入进去并应用
npower=noisepow(noise_bw);
threshold=npower*db2pow(npwgnthesh(pfa,int_pulsenum,'noncoherent'));
threshold=threshold*db2pow(mfgain);
最新发布