2024.6.4
理论部分
算法用途
Keystone变换是一种雷达信号处理技术,主要用于校正目标距离走动现象。
原理
这项技术的主要原理是将接收到的雷达信号从频域转换到Keystone频谱,然后进行时间徙动校正,以消除距离徙动的影响。
在Keystone变换中,首先对原始的多普勒频谱进行转换,得到Keystone频谱,接着对这个频谱进行时间徙动校正,最终得到校正后的多普勒频谱。这种变换在距离-多普勒处理中非常重要,它能够提高雷达目标检测和定位的精度。
Keystone变换的实现方法有多种,包括DFT+IFFT算法、Chirp-Z变换法、sinc函数内插法和三阶Lagrange内插法等。这些方法各有特点,适用于不同的应用场景。例如,DFT+IFFT算法和Chirp-Z变换法在物理意义上虽然不同,但具体实现过程是相同的。
此外,Keystone变换不仅适用于脉冲雷达(PD雷达),也适用于宽带阵列数据的预处理,通过将各频点的相位对齐到一个基准频率上,从而提升目标成像和高速机动目标检测的精度。
具体可见:https://blog.csdn.net/m0_54016576/article/details/138974645?spm=1001.2014.3001.5502
MATLAB代码
%{
1) 测试KeyStone变换插值结果是否正确
2) 使用DFT+IFFT,sinc插值和Chirp-Z变换三种方法
%}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all; close all; clc;
dbstop if error;
%---------------------系统参数设置-----------------------------------------
settings.fc = 1e9; % 载波频率 --- 1GHz
settings.c = 3e8; % 光速
settings.lambda = settings.c/settings.fc; % 波长
% 信号参数
settings.PRF = 2e3; % PRF --- 2KHz
settings.PRI = 1/settings.PRF; % PRI
settings.B = 15e6; % 带宽 --- 15MHz
settings.tau = 10e-6; % 脉宽 --- 10us
settings.mu = settings.B/settings.tau; % 调频率
settings.BT = settings.B*settings.tau; % 时宽带宽积
settings.Pm = 512; % 慢时间采样点数
settings.CPI = settings.Pm*settings.PRI; % CPI
settings.d = settings.tau/settings.PRI; % 占空比
settings.SegNum = 32; % 分段数
settings.SegLen = settings.Pm/settings.SegNum; % 段内脉冲数
settings.Tp = settings.SegLen*settings.PRI; % 段内时间
%Pm_mat = reshape(1:settings.Pm, [settings.SegLen, settings.SegNum]);
%settings.SegInd = Pm_mat(1,:); % 段落起始脉冲索引
% 天线数目
settings.TraNum = 1; % 发射天线数目
settings.RecNum = 4; % 接收天线数目
% 快时间采样频率
settings.fs = 40e6; % 采样频率
settings.ts = 1/settings.fs; % 采样间隔
settings.samples = floor(settings.tau*settings.fs);
% 脉内采样点数
% 1) 考虑没有速度模糊的情况
%----------------------构建S(t,f)矩阵--------------------------------------
Va = settings.PRF*settings.lambda/4; % 150m/s
Num_f = 1024; % range sampling cells
Num_s = settings.Pm;
Rmax = (Num_f - settings.samples)*settings.c/settings.fs/2;
% 对回波斜距进行建模
R0 = 1e3; % 初始位置
Vt = 100; % 径向速度
R_t = R0 + Vt.*(-settings.Pm/2:settings.Pm/2-1)*settings.PRI;
if max(R_t) > Rmax || min(R_t) < 0
disp('-_-');
return
end
%-------------------------产生回波-----------------------------------------
Time_f = linspace(0,(Num_f-1)/settings.fs, Num_f) - settings.tau/2;
% 减掉tau/2是因为现在信号快时间的范围是-t/2到t/2
% 计算延时
tau = 2*R_t./settings.c;
% 回波延时矩阵
Ti_mat = repmat(Time_f, Num_s,1) - tau';
% 计算回波相位(不考虑脉内多普勒频移)
Phasepoint = pi*settings.mu*Ti_mat.^2 - 2*pi*settings.fc*tau';
% 生成回波采样矩阵
Sr_mat(:,:) = exp(1i.*Phasepoint) ...
.*(-settings.tau/2<Ti_mat&Ti_mat<settings.tau/2);
%----------------------对回波信号进行脉压----------------------------------
% 参考信号建模
t = linspace(-settings.tau/2,settings.tau/2,settings.samples);
s_ref = exp(1i*pi*settings.mu.*t.^2);
s_ref = [s_ref, zeros(1,Num_f - settings.samples)];
% 参考信号频域模型
H_ref = fftshift(ones(Num_s,1)*fft(s_ref),2);
% 接收回波脉冲压缩
S_f = fftshift(fft(Sr_mat,[],2),2);
St_f = S_f.*conj(H_ref);
% IFFT转换到(tm,t)上
St_r = ifft(St_f,[],2);
% 画图
figure(1)
imagesc(abs(St_r));
xlabel('range cell');
ylabel('pulse index');
axis([240 290 -inf inf]);
%-------------------DFT+IFFT方法-------------------------------------------
% 快时间频率
fr = linspace(-settings.fs/2,settings.fs/2,Num_f);
% step1: DFT的过程
%
alpha = (settings.fc + fr)./settings.fc;
S_lk = zeros(Num_s,Num_f);
for k = -Num_s/2:Num_s/2 - 1
tem = 0;
for n = -Num_s/2:Num_s/2 - 1
tem = tem + St_f(k+Num_s/2+1,:).*exp(-1i*2*pi.*alpha*k*n/Num_s);
end % for m = -Num_s/2:Num_s/2 - 1
S_lk(k+Num_s/2+1,:) = tem;
end % for index = -settings.Pm/2:settings.Pm/2-1
% step2: IFFT的过程
% 直接用IFFT函数计算
S_lm = ifft(S_lk);
s_lm = ifft(S_lm,[],2);
figure(2)
imagesc(abs(s_lm));
xlabel('range cell');
ylabel('pulse index');
axis([240 290 -inf inf]);
title({'DFT+IFFT方法','对齐到轨迹中点'});
%
% %
% % %------------使用sinc插值的方法完成KeyStone变换--------------------------
% % 开始sinc插值
Stm_f = zeros(Num_s,Num_f);
fr = linspace(-settings.fs/2,settings.fs/2,Num_f);
%
row_vec = 1:Num_s;
% 脉压频域结果沿快时间FFTSHIFT
% St_f = fftshift(St_f,2);
for n = 1:Num_f
theta = settings.fc/(fr(n) + settings.fc);
for row = 1:Num_s
Stm_f(row,n) = sinc(theta*(row-1) - (row_vec-1))*St_f(:,n);
end % for row = 1:Num_s
end % for n = 1:Num_f
% 对f做IFFT
Stm_r = ifft(Stm_f,[],2);
figure(3)
imagesc(abs(Stm_r));
xlabel('range cell');
ylabel('pulse index');
axis([240 290 -inf inf]);
title('sinc插值方法');
%-------------------Chirp-Z变换方法----------------------------------------
%--- 分块调试
fr = linspace(-settings.fs/2,settings.fs/2,Num_f);
% 构造L
L = 2^ceil(log2(2*Num_s -1));
% 构造W
phi0 = ((settings.fc + fr)./settings.fc)*(2*pi/Num_s);
W = exp(1i*phi0);% exp(-1i*phi0);
% 开辟G和H矩阵的存储空间
G = zeros(L,1);
H = zeros(L,1);
% 变换后结果的存储空间
X_f = zeros(Num_s,Num_f);
n = 1:Num_s;
n = n';
n2 = Num_s+1:L;
n2 = n2';
for fk = 1:Num_f
% 构造G矩阵
G(1:Num_s) = St_f(:,fk).*(W(fk).^(n.^2/2));
% 构造H矩阵
H(1:Num_s) = W(fk).^(-n.^2/2);
H(Num_s+1:end) = W(fk).^(-(L-n2).^2/2);
% 傅里叶变换
Gf = fft(G);
Hf = fft(H);
% 计算序列v
v = ifft(Gf.*Hf);
% 用v的前Num_s个值作为权值,代入公式
X_f(:,fk) = v(1:Num_s).*(W(fk).^(n.^2/2));
% 直接对X_f做逆变换得到的结果是错的
% z_k是z平面内的采样点,我们需要的是x(n)序列
% 所以应该对X_f在慢时间域上做IFFT
% 然后再从快时间频率f上做IFFT
end % for fk = 1:Num_f
X_tf = ifft(X_f,[],1);
X_k = ifft(X_tf,[],2);
figure(4)
imagesc(abs(X_k));
xlabel('range cell');
ylabel('pulse index');
axis([240 290 -inf inf]);
title({'Chirp-Z变换方法',['Vt = ', num2str(Vt)]});
% 2) 考虑存在速度模糊的情况
% --- testKeyStone2---------