1.引言
FMCW雷达信号处理建模与仿真,是诸多初学者在研究雷达的过程中一定会遇到的问题,很多初学者甚至不知道如何下手,因此调皮哥针对这个问题,给大家整理了本期文章。本文将带领大家学会如何基于FMCW进行雷达信号处理建模仿真,并结合雷达信号处理理论与MATLAB编程实践学习,从发射信号开始,到回波信号、距离维FFT、速度维FFT以及CFAR检测等内容,文章提供了对应的MATLAB代码帮助大家学习仿真,其中CFAR部分放在下期文章详细分析。
调皮哥之前写过《调皮连续波:雷达原理 | 用MATLAB信号处理是如何解算目标的距离和速度信息的?》一文,但是太过于理论,并没有从雷达的根源角度系统性地介绍雷达的基本原理,因此本期文章将结合该文更加系统地介绍,能够帮助大家形成知识体系。
2.发射信号模型
一般情况下,雷达发射信号的模型可采用线性调频连续波(LFMCW) ,发射波形的信号形式为调频连续锯齿波。线性调频的含义即调制信号频率随时间线性变化,从时域上看是一个频率随时间线性变化的波形;从频域上看, 发射信号的频率与时间成正比, 如图1所示。
图1 线性调频连续波(LFMCW)
理论上,发射信号的模型可以用公式(1)表示: 其中A是发射信号幅度,发射的chirp信号扫频带宽为 B,发射时宽为T,即调频斜率为 B/T,记为S 。所以,单周期的 FMCW 雷达发射信号模型的相位可以写成公式(2):
所谓单周期,即只有一个chirp扫频信号。但是一般我们研究FMCW雷达信号处理,需要多周期(如128个)的chirp信号,这是因为我们需要利用多个chirp来测量目标的速度。
通常,各种文献资料对于雷达发射信号的表示方法有两种,一种是实数信号表示形式,如本文所用的。另一种是复数信号表示形式,如公式(3)所示:
但是由于我们在发射信号是,只有实数信号,而并非复数信号,只是在后续的ADC芯片中采用正交采样时,才把信号变为复信号。需要注意的是,这两种表达方式都是正确的,只是在后续的处理过程中,有一些不同。本文仍旧采用传统的实数信号表达方式。
在仿真MATLAB代码中发射信号的理想模型表示为:
Tx(i) = cos(2*pi*(fc*t(i) + (slope*t(i)^2)/2)); % 发射信号 实数信号
发射信号时域图如图2所示:
plot(Tx(1:1024));xlabel('点数');ylabel('幅度');title('TX发射信号时域图');
图2 发射信号时域图
发射信号时频图如图3所示:
plot(t(1:1024),freq);xlabel('时间');ylabel('频率');title('TX发射信号时频图');
图3 发射信号时频图
3.雷达参数与目标参数设置
在MATLAB代码中设置雷达的参数和目标的参数,便于我们进行建模仿真实验。具体参数如下所示:
maxR = 200; % 雷达最大探测目标的距离rangeRes = 1; % 雷达的距离分率maxV = 70; % 雷达最大检测目标的速度fc= 77e9; % 雷达工作频率 载频c = 3e8; % 光速
r0 = 90; % 目标距离设置 (max = 200m)v0 = 10; % 目标速度设置 (min =-70m/s, max=70m/s)
B = c / (2*rangeRes); % 发射信号带宽 (y-axis) B = 150MHzTchirp = 5.5 * 2 * maxR/c; % 扫频时间 (x-axis), 5.5= sweep time should be at least 5 o 6 times the round trip timeendle_time=6.3e-6; %空闲时间slope = B / Tchirp; %调频斜率f_IFmax= (slope*2*maxR)/c ; %最高中频频率f_IF=(slope*2*r0)/c ; %当前中频频率
Nd=128; %chirp数量 Nr=1024; %ADC采样点数vres=(c/fc)/(2*Nd*(Tchirp+endle_time));%速度分辨率Fs=Nr/Tchirp; %模拟信号采样频率t=linspace(0,Nd*Tchirp,Nr*Nd); %发射信号和接收信号的采样时间
Tx=zeros(1,length(t)); %发射信号Rx=zeros(1,length(t)); %接收信号Mix = zeros(1,length(t)); %差频、差拍、拍频、中频信号
4.回波信号模型
假设静止的目标距离雷达的距离为R,电磁波在空气中传输速度为c,则接受天线接受到的信号比发射的信号延迟τ=2R/c,所以理想情况下接受天线接收到的目标回波信号模型如公式(4)所示:
由上述公式可知,回波信号具有和发射信号同样的信号形式,相对于发射信号在时间上有固定延时τ,故而回波信号的相位可以表式为公式(5):
在MATLAB中接收信号(雷达的回波信号)表示为:
Rx(i) = cos(2*pi*(fc*(t(i)-td(i)) + (slope*(t(i)-td(i))^2)/2)); %接收信号 实数信号
接收信号的时域图如图4所示:
%接收信号时域图plot(Rx(1:1024));xlabel('点数');ylabel('幅度');title('RX接收信号时域图');
图4 接收信号的时域图
接收信号与发射信号的时频图如图5所示:
%接收信号与发射信号的时频图plot(t(1:1024),freq);hold on;plot(td(1:1024)+t(1:1024),freq);xlabel('时间');ylabel('频率');title('接收信号与发射信号时频图');legend ('TX','RX');
图5 接收信号与发射信号的时频图
4.混频
将接收到回波信号Sr(t)和发射信号St(t)进行混频,并经过低通滤波器后就可以得到一个单一频率的正弦波信号,如图中黑色的“IF signal”便是中频信号的频率,如图6和图7所示。
图6 混频乘法器
图7 中频信号模型
4.1 静止目标下测距
根据公式推导,当在静止目标下时, 可以得到差频信号的相位表达式为公式(6):
此时可以很明显地看出,发射信号和单目标的回波信号的频率差为一个单频信号。根据公式(6),可以得到中频信号频率 ,如公式(7)所示:
对中频信号进行 ADC 采样, 然后做 FFT 提取信号的频率信息, 假设 FFT 得到频谱的谱峰值对应的频率为fm , 则目标的距离信息可以表示为公式(8):
4.2 运动目标情况下测距
假设,在电磁波的覆盖区域中,存在某一目标在 时刻距离发射天线为 ,以径向速度 远离雷达,以远离天线为正方向,至于为什么可看文章(《调皮连续波:雷达原理 | 讨论调频连续波雷达目标运动方向与速度正负的关系?》)
那么接收到目标的回波信号模型公式依旧与单目标一致,如公式(9)所示。
但是,τ有所改变,如公式(10)所示:
此时,通过混频后得到中频信号的相位如下所示:
很显然,根据上述公式可知,对于运动目标的信号的中频信号依然是一个近似线性调频信号。所以设中频信号的调频斜率Um,载频fm,初相ϕm分别如下:
因为光速c等于3*10^8m/s,所以忽略c的平方项,则中频信号的时宽带宽积Dm为:
其中,Bm中频信号的调频带宽,为D为发射信号时宽带宽积。上述公式可以表明,即使目标在几百米每秒的高速运动情况下,中频信号的时宽带宽积仅有原来的10^-6 倍,在毫米波雷达发射极大时宽带宽积(10^6)的信号情况下,中频信号Dm的数量级也只有10^0 。因此,可近似认为回波差拍信号是一单频信号,通过频谱分析(FFT)即能得到其中心频率。所以也可以将中频信号近似地写成公式(11):
注意:公式(11)是近似值,忽略了多普勒频移的影响,其中R为目标距离, , 。
近似中频信号公式忽略中频信号的频率与物体速度的依赖关系。在快速 FMCW 雷达中,其影响通常非常小,且在处理完成多普勒 (速度维)FFT 后,即可轻松对其进行进一步校正,这就是多普勒相位补偿。
混频过程在MATLAB中可以表示为:
Mix(i) = Tx(i).*Rx(i);%差频、差拍、拍频、中频信号
混频后的信号频图,为实信号做傅里叶变换,如图8所示:
plot(db(abs(fft(Mix(1:1024)))));%查看宽带的和频信号 将chirp的点数改为1024*256即可看到有一个门信号,但注意计算机内存。xlabel('频率');ylabel('幅度');title('中频信号频谱');
图8 中频信号频谱(实信号)
此时要注意,本文采用“数字信号的形式”来模拟发射信号与回波信号,实际上是不会通过这种方式实现的,发射信号和回波信号都是模拟信号。因此在本文中,混频之后得到的中频信号只有单一的频谱,即只保留了混频乘法器积化和差公式后的“差频”信号,“和频”信号由于MATLAB采用数字信号的形式导致采样率过低,N/Tchirp=1024/Tchirp=139MHz,对“和频信号”相当于欠采样,因此并不会保留高频信号,只保留目标的最大中频信号(27.27MHz),从而MATLAB自身形成一个低通(或带通)滤波器。
如果读者想要看到“和频”信号,则可以通过增加发射信号点数的点数,相当于增加采样率的形式来观察,但由于数据量很大,雷达载波也很高,因此不建议大家尝试。调皮哥给大家尝试了,在MATLAB代码中约为35GHz,离77GHz还很远,相当于欠采样,但还可以采集到的部分高频信号。如图9中所示,类似于门函数的频带就是和频信号。
Nr=1024*256; %中频信号频谱plot(db(abs(fft(Mix(1:1024*256)))));%查看宽带的和频信号 将chirp的点数改为1024*256即可看到有一个门信号,但注意计算机内存。xlabel('频率');ylabel('幅度');title('中频信号频谱');
图9 中频信号(差频以及和频)
5.低通(带通滤波)
本文的FMCW雷达信号处理建模仿真,由于MATLAB采样点数问题,低通滤波器相当于自带。如若按照增加采样点数的方式来实现低通滤波器的模拟,则可以实现如图10~12所示。
%% 低通滤波 截止频率30MHz 采样频率120MHzMix=lowpass(Mix(1:1024*256),30e6,120e6);plot(db(abs(fft(Mix(1:1024*256)))));xlabel('频率');ylabel('幅度');title('中频信号低通滤波器');
图10 滤波后结果
图11 低通滤波幅频响应
图12 低通滤波结果
放大后保留的频率如图13所示:
图13 低通滤波结果(放大效果)
6.正交采样(IQ复采样)
本文暂采用实数信号进行分析,因此暂不对实信号进行正交采样,后续有时间更新在进行处理。
7.距离估计
假设单个扫描周期 ADC采样点数为 1024, 采样频率为Fs=139MHz, 式中 n 为 FFT 谱峰对应的频点(Index), 根据 T、 Fs、 N 之间的关系,距离与频率关系的公式可进一步化简为:
其中,请注意,因为上述的采样点刚好是2的次幂,因此,计算之后系数K刚好为1,假如采样点数不是2的次幂,如200个点,然后进行256个点的FFT,那么K值就不能等于1。
将中频信号按照距离门点数、chirp脉冲数排列为二维矩阵,其时域谱如图14所示。
signal = reshape(Mix,Nr,Nd);
图14 距离门点数、chirp脉冲数排列为二维矩阵(时域)
第一个chirp脉冲的距离维FFT结果,实信号FFT后频谱对称,只保留一半的频谱,即512个点,如图15所示。
%% 距离维FFTsig_fft = fft(signal,Nr)./Nr;sig_fft = abs(sig_fft);sig_fft = sig_fft(1:(Nr/2),:);
figure;plot(sig_fft(:,1));xlabel('距离(频率)');ylabel('幅度')title('第一个chirp的FTF结果')
图15 第一个chirp脉冲的距离维FFT结果
距离维FFT结果谱矩阵,如图16所示。
figure;mesh(sig_fft);xlabel('距离(频率)');ylabel('chirp脉冲数')zlabel('幅度')title('距离维FTF结果')
图16 距离维FFT结果
在MATLAB中,对于距离维做FFT,目的就是找出中频频率峰值和频点,然后根据峰值所在的距离门号,解算出目标的距离。
如本文的距离门号为91,距离分辨率为1米,则真实距离为:(91-1)*1=90米,符合预先设定的目标参数值。
8.速度估计
对于两个相邻周期的信号,由于周期间隔时间Tc较短,距离分辨率有限,两个周期内距离维FFT谱中的峰值位置几乎没有发生变化,但是由于相位比距离更加敏感,即周期间微小的距离变化会引起中频信号的初相的变化。由傅里叶变换特性,信号的初相体现在峰值处的复数值对应的相位,计算相邻周期的相位差,即可得到目标的速度:
其中,两个脉冲相邻的微小距离变化为v*Tc,所以推导出两个脉冲间的速度为: 推广到多个脉冲,如128个,那么相位的变化是呈周期性的,速度维度做128点FFT后的峰值就是相位差,即获取速度对应的峰值和频点再通过解算便可获得目标的速度,这个过程就是速度维FFT。接下来进行速度解算。上述公式经过推导可以变为:
其中,相位的最大值为2π。由此,速度的解算依赖于速度分辨率和速度门号,速度门号即为速度维FFT之后拿到的峰值所对应的下标(相位变化的频点)。
其中,请注意,因为上述的采样点刚好是2的次幂,因此,计算之后系数k刚好为1,假如采样点数不是2的次幂,如200个点,然后进行256个点的速度维FFT,那么K值就不能等于1。
在MATLAB中的解算过程为:
Mix=reshape(Mix,[Nr,Nd]);sig_fft2 = fft2(Mix,Nr,Nd);
sig_fft2 = sig_fft2(1:Nr/2,1:Nd);sig_fft2 = fftshift (sig_fft2);RDM = abs(sig_fft2);RDM = 10*log10(RDM) ;doppler_axis = linspace(-100,100,Nd);range_axis = linspace(-200,200,Nr/2)*((Nr/2)/400);figure;mesh(doppler_axis,range_axis,RDM);xlabel('距离通道'); ylabel('多普勒通道'); zlabel('幅度(dB)');title('速度维FFT 距离多普勒谱');
速度维FFT的结果如图17所示:(订:图中速度维和距离维写反,注意辨别)
图17 速度维FFT结果
如本文的多普勒门号为10,速度分辨率为 vres= 1.1163m/s,则真实距离为:(10-1)* 1.1163=10.04m/s,符合预先设定的目标参数值,其中这里的速度分辨率vres是由雷达工程师自主设定的,在FMCW雷达中可以增加调频连续波的空闲时间来增加速度分辨率,但要注意与最大不模糊速度保持平衡。