线性调频信号的解线调(dechirp，去斜)处理matlab仿真

线性调频信号的回波模型

对于线性调频脉冲压缩雷达，其发射信号为：
$$\begin{array}{c}s(t)=rect(\frac{t}{T})exp(j\pi \mu t^2)\end{array}\text{\hspace{1em}?}$$
其中
$$\begin{array}{c}rect(\frac{t}{T})=\{\begin{array}{ll}1,& 0\le t\le T\\ 0,& otherwise.\end{array}\end{array}\text{\hspace{1em}?}$$
其中$t$为快时间变量，$T$为脉冲宽度，$\mu$为调频斜率，$B=T\mu$为雷达信号的带宽，搭载到载频上的发射信号为：
$$\begin{array}{c}s(t)=rect(\frac{t}{T})exp(j\pi \mu t^2)exp(j2\pi f_{c}t)\end{array}\text{\hspace{1em}?}$$
假设单目标的情况，目标位于$R_{tar}$处，接收到的回波信号表达式为
$$\begin{array}{c}{s}_r(t)=rect(\frac{t?t_0}{T})exp(j\pi \mu (t?t_0{)}^2)exp(j2\pi f_c(t?t_0))\end{array}\text{\hspace{1em}?}$$
其中$t_0=2R_{tar}/c$。

参考信号

参考信号同样是一个线性调频信号，调频斜率、载频与雷达发射的信号一样，但是脉宽大于雷达发射线性调频信号的$T$，这个脉宽其实就是一个信号处理窗，我们接收到一个脉冲重复周期的雷达回波之后，对选定的一段信号进行处理，这一段信号的长度就是窗，因为我们不知道雷达回波中的脉冲位置，所以选取的窗应该至少等于脉冲宽度才能将其“覆盖”住。参考信号可以表示为：
$$\begin{array}{c}{s}_{ref}(t)=rect(\frac{t?t_{ref}}{T_{ref}})exp(j\pi \mu (t?t_{ref}{)}^2)exp(j2\pi f_c(t?t_{ref}))\end{array}\text{\hspace{1em}?}$$
其中，$t_{ref}=2R_{ref}/c$，$c$为光速，$R_{ref}$是参考距离，$T_{ref}$是接收窗，一般而言$T_{ref}>T$。

图1 雷达回波和参考信号示意图

去斜处理

解线调就是将雷达回波和参考信号的共轭进行混频
$$\begin{array}{c}{s}_{de}(t)=s_r(t)?s_{ref}^{?}(t)=rect(\frac{t?t_0}{T})exp(j2\pi (?\mu (t_0?t_{ref}))t)exp(?f_c(t_0?t_{ref}))exp(1/2\mu (t_0^2?t_{ref}^2)))\end{array}\text{\hspace{1em}?}$$
上面式子中一共有3个指数项，且只有第一个指数项是时间$t$的函数，后面两个记为$exp(j?)$，则可以化简为：
$$\begin{array}{c}{s}_{de}(t)=s_r(t)?s_{ref}^{?}(t)=rect(\frac{t?t_0}{T})exp(j2\pi (?\mu (t_0?t_{ref}))t)exp(j?)\end{array}\text{\hspace{1em}?}$$
可以看到解线调后的信号是一个单频信号，其频率为$f=?\mu (t_0?t_{ref})$，它是调频斜率$?\mu$和相对时间$(t_0?t_{ref})$的函数，$(t_0?t_{ref})$也就是雷达回波的脉冲相对于参考信号时间窗的延时。$t_{ref}$是已知的，那么由解线调信号的频率便可以得到目标的延时。

去斜处理傅里叶变换得到脉压结果

我么对解线调信号进行傅里叶变换
S d e ( f ) = F { s d e ( t ) } = ∫ s d e ( t ) e x p ( ? j 2 π f t ) d t = T r e f e x p ( j φ ( f ) ) s i n c ( T r e f ( f ? f t a r ) ) \begin{equation} S_{de}(f)=\mathcal{F}\{s_{de}(t) \}=\int s_{de}(t)exp(-j2\pi ft)dt=T_{ref}exp(j\varphi (f))sinc(T_{ref}(f-f_{tar})) \end{equation} Sde?(f)=F{sde?(t)}=∫sde?(t)exp(?j2πft)dt=Tref?exp(jφ(f))sinc(Tref?(f?ftar?))??
其中$f_{tar}$和$R_{tar}$具有以下关系：
$$\begin{array}{c}{R}_{tar}=R_{ref}+(c?f_{tar})/(2\mu )\end{array}\text{\hspace{1em}?}$$

整个流程如下图所示，解线调之后由频域图可以得到和匹配滤波相同的结果，它们都是sinc函数。

图2 解线调流程图及时频图

解线调仿真

下面仿真的例子：
线性调频信号的带宽为2MHz，采样频率为8MHz，雷达发射线性调频信号脉宽50us，解线调脉宽100us。

%%Dechirp仿真，解线性调频
% 雷达信号采用线性调频信号，
% 主要内容：线性调频信号的生成、雷达回波的模拟、解线调参考信号生成、解线调操作、结果绘制
% Author: baidan @beijing 2023.1.11
%=========================================================================%
%                        雷达参数设置                                     %
%=========================================================================%
clear all;close all;clc;
C = 3.0e8;              %光速（m/s）
BandWidth = 2.0e6;      %雷达发射信号带宽，带宽=B=1/tau，tau是脉冲宽度
TimeWidth = 50.0e-6;    %雷达发射信号的脉冲时宽
k = BandWidth/TimeWidth;%调频斜率
PRT = 100e-6;           %雷达发射脉冲重复周期（s），100us对应1/2*100*300=15000米最大无模糊距离
Fs = BandWidth*4;       %采样频率 采样点数为PRT*FS=480
NoisePower = 1;         %噪声功率
SNR = 0;                %信噪比，以dB为单位
SigPower = 10^(SNR/10); %目标功率，无量纲
TargetDistance = 2500;  %目标距离，单位：
TargetDelay = 2*TargetDistance/C;           %延时时间，由目标的距离换算成时间，即线性调频信号的时间
mDelayNumber = fix(Fs*2*TargetDistance/C);  %把目标距离换算成采样点（距离门）
plot_enable_H = 1;      %绘图控制符
plot_enable_L = 0;      %绘图控制符
%=========================================================================%
%                            调用函数产生线性调频信号                     %
%=========================================================================%
[LFMPulse,targetEchoPRT,matchedFilterCoeff,pulseNumber,PRTNumber] = GenerateLFMSignal(BandWidth,TimeWidth,PRT,Fs,SigPower,TargetDistance,plot_enable_L);%调用函数
fprintf('Target appear in the %.1f us, distence is %.1f m, at %dth sampling point\n',TargetDelay*1e6,TargetDistance,mDelayNumber);
%=========================================================================%
% %                        产生功率为1高斯白噪声并添加噪声                  %
% %=========================================================================%
noise = GaussianNoise(PRT,Fs);
%mixedSignal = targetEchoPRT + noise; %添加噪声
mixedSignal = targetEchoPRT; %不添加噪声
%=========================================================================%
%                                   产生参考信号                          %
%=========================================================================%
[referSig,sigSamNumber] = ReferenceSignal(k,PRT,Fs,plot_enable_H); %四个参数分别是调频斜率，脉宽，采样率，绘图控制符
sref = referSig; %进行共轭
sde = targetEchoPRT.*(sref); %雷达回波和参考信号相乘
Sde = fftshift(fft(sde))*2/1600; % 对dechirp之后的信号进行fft
f=Fs*(0:PRTNumber-1)/PRTNumber-Fs/2;% -fs/2到 fs/2
figure;
subplot(4,1,1);plot(real(targetEchoPRT),'linewidth',1);axis tight;title('Target echo '); %绘制雷达回波
subplot(4,1,2);plot(real(sref),'linewidth',1);axis tight;title('Reference signal'); %绘制参考信号
subplot(4,1,3);plot(real(sde),'linewidth',1);axis tight;title('Dechirp signal'); %绘制参考信号
subplot(4,1,4);plot(f,(abs(Sde)),'linewidth',1);axis tight;title('fft of Dechirp signal'); %绘制参考信号

其中的子函数，产生线性调频信号：

function [LFMPulse,targetEchoPRT,matchedFilterCoeff,pulseNumber,PRTNumber] = GenerateLFMSignal(bandWidth,pulseDuration,PRTDuration,samplingFrequency,signalPower,targetDistece,plotEnableHigh)
% 该函数用于产生线性调频信号，以及雷达的目标反射回波，仅产生单个回波
%  Author:壮神 2023.9.21 @Beijing
% Input : 
%   * bandWidth: 信号带宽 ，参考值：2.0e6 表示2MHz
%   * pulseDuration：脉冲持续时间，参考值：40.0e-6 表示40ms
%   * PRTDuration：脉冲重复周期，参考值：240ms
%   * samplingFrequency：采样频率，参考值：2倍的信号带宽
%   * signalPower：信号能量，参考值：1
%   * targetDistece：目标距离，最大无模糊距离由脉冲重复周期决定。计算公式：1/2*PRTDuration*光速
%   * plotEnableHigh： 绘图控制符，1：打开绘图，0：关闭绘图
% Output : 
%    * LFMPulse：线性调频信号
%    * targetEchoPRT： 目标反射回波
%    * matchedFilterCoeff： 匹配滤波器系数
%    * pulseNumber：当前采样率下线性调频信号的采样点数
%    * PRTNumber：1个PRT对应的采样点数
C = 3.0e8;      %光速（m/s）
BandWidth = bandWidth;  %雷达发射信号带宽，带宽=B=1/tau，tau是脉冲宽度
TimeWidth = pulseDuration; %雷达发射信号的脉冲时宽

PRT = PRTDuration;       %雷达发射脉冲重复周期（s），240us对应1/2*240*300=360000米最大无模糊距离
Fs = samplingFrequency;         %采样频率
SampleNumber = fix(Fs*PRT);
%=========================================================================%
%                        目标参数设置                                     %
%=========================================================================%
SigPower = signalPower;           %目标功率，无量纲
TargetDistance = targetDistece; %目标距离，单位：m
DelayNumber = fix(Fs*2*TargetDistance/C); %把目标距离换算成采样点（距离门）
fprintf('目标出现在距离门%d\n',DelayNumber);
%=========================================================================%
%                        产生线性调频信号、匹配滤波器                     %
%=========================================================================%
number = fix(Fs*TimeWidth); %回波采样点数=脉压系数长度=暂态点数目+1
if rem(number,2)~=0
    nember = nember + 1;
end
Chirp = zeros(1,number);
for i = -fix(number/2):fix(number/2)-1
    Chirp(i+fix(number/2)+1)=exp(1j*(pi*(BandWidth/TimeWidth)*(i/Fs)^2));%产生复ChIrp信号
end
coeff = conj(fliplr(Chirp)); %把Chirp信号翻转并把复数共轭，产生脉压系数
%=========================================================================%
%                      绘制线性调频信号                                   %
%=========================================================================%
if plotEnableHigh == 1
    figure;
    plot(real(Chirp)); %绘制线性调频信号
    xlabel('Sampling points'); ylabel('Amplitude');title('线性调频信号实部');
end
SignalTemp = zeros(1,SampleNumber); %1个PRT
SignalTemp(DelayNumber+1:DelayNumber+number) = sqrt(SigPower)*Chirp;%将线性调频信号按照距离进行延时
if plotEnableHigh == 1
    figure;
    plot(real(SignalTemp)); %绘制1个完整的PRT的雷达回波信号
    xlabel('Range bin'); ylabel('Amplitude');title('雷达回波的实部');
end
%=========================================================================%
%                          进行脉冲压缩                                   %
%=========================================================================%
Echo = SignalTemp; % 目标回波
pc_time0 = conv(Echo,coeff); % 回波和滤波器卷积的结果
pc_time1 = pc_time0(number:number+SampleNumber-1); %去掉暂态点
realTargetRange = find(abs(pc_time1)==max(abs(pc_time1)))-1; %由脉压结果目标距离
fprintf('The target range bin is  %d',realTargetRange);
if plotEnableHigh == 1
    figure; %时域脉压结果
    subplot(2,1,1);plot(abs(pc_time0),'r-');
    xlabel('Range bin'); ylabel('Amplitude');title('时域脉压结果');
    subplot(2,1,2);plot(abs(pc_time1),'r-');
    xlabel('Range bin'); ylabel('Amplitude');title('去掉暂态点的时域脉压结果');
end
%=========================================================================%
%                              返回参数                                   %
%=========================================================================%
LFMPulse = Chirp; %线性调频信号
targetEchoPRT = SignalTemp; %目标反射回波
matchedFilterCoeff = coeff; %匹配滤波器系数
pulseNumber = number; %线性调频信号的采样点数
PRTNumber = SampleNumber; %目标反射回波的采样点数
end

其中的子函数，产生解线调参考信号：

function [referSig,sigSamNumber] = ReferenceSignal(chirpRate,windowDuration,samplingFrequency,plotEnableHigh)
% 该函数用于产生dechirp雷达的参考信号
%  Author:王壮壮 2024.1.14 @shenzhen
% Input : 
%   * chirpRate: 调频斜率
%   * windowDuration：接收信号的时间窗，一般要大于雷达的线性调频信号的脉宽
%   * samplingFrequency：采样频率，参考值：2倍的信号带宽
%   * plotEnableHigh： 绘图控制符，1：打开绘图，0：关闭绘图
% Output : 
%    * referSig：参考信号，也是一个线性调频信号，但是脉宽要大于雷达发射的LFM信号
%    * pulseNumber：当前采样率下线性调频信号的采样点数
C = 3.0e8;      %光速（m/s）
TimeWidth = windowDuration; %接收窗的时宽
Fs = samplingFrequency;         %采样频率
%=========================================================================%
%                              产生线性调频信号                           %
%=========================================================================%
number = fix(Fs*TimeWidth); %采样点数
if rem(number,2)~=0
    nember = nember + 1;
end
% Chirp = zeros(1,number);
for i = -fix(number/2):fix(number/2)-1
    Chirp(i+fix(number/2)+1)=exp(1j*(pi*chirpRate*(i/Fs)^2));%产生复ChIrp信号
end
%=========================================================================%
%                              返回参数                                   %
%=========================================================================%
referSig = conj(Chirp); %对复chirp信号进行共轭得到参考信号
sigSamNumber = number; %参考信号的采样点数
%=========================================================================%
%                        分析参考信号的频域                               %
%=========================================================================%
CHIRP = fftshift(fft(referSig))*2/sigSamNumber;  %对参考信号进行fft得到它的频域
frequency_index = Fs*(0:sigSamNumber-1)/sigSamNumber-Fs/2; %频率下标
%=========================================================================%
%                        绘制参考信号                                     %
%=========================================================================%
if plotEnableHigh == 1
    figure;
    subplot(2,1,1);plot((1:sigSamNumber),real(Chirp)); %绘制参考信号的时域波形
    xlabel('Sampling points'); ylabel('Amplitude');title('参考信号实部');
    subplot(2,1,2);plot(frequency_index,abs(CHIRP)); %绘制参考信号的频谱图
    xlabel('Frequency/Hz'); ylabel('Amplitude');title('参考信号的频谱');
end
end

绘制图形如下：

图3 参考信号时频图

图4 解线调信号时频图

总结

去斜处理的主要原理为用发射副本信号与回波混频得到基带信号，再用傅里叶变换提取基带信号频率，通过时域到频域的转换将时域脉冲信号压缩为频域上的一个峰。匹配滤波和解线调其实都是一个sinc函数乘以一个复指数项的形式，两者的sinc函数一个位于快时间域一个位于距离频率域。去斜处理在一些特殊场合，不仅运算简单，而且可以简化设备，已经广泛应用于SAR和ISAR中作脉冲压缩。