【原】基于Python的FMCW雷达工作原理仿真（附代码）

调皮连续波 2022-10-04 发布于贵州

展开全文

本文编辑：@调皮连续波，保持关注调皮哥，获得更多雷达学习资料和建议！

大家好，我是调皮哥，今天继续给大家分享干货，助力大家轻松、快乐、有方向地学习雷达。

通常，我们都是采用Matlab进行雷达系统设计，或者进行相关的算法仿真。Matlab的确是一个十分强大的工具，但随着越来越多的雷达er开始使用Python进行雷达设计，我们不得不对基于Python的FMCW雷达信号处理仿真做一次详细的介绍，以便于后来的雷达初学者能够了解到更全面的信息，开拓视野。

本文的主要内容如下：

1.距离检测原理和仿真
2.速度检测原理和仿真

一、距离检测

(一)距离检测原理

调制连续波雷达也称为调频连续波(FMCW) 雷达，是一种使用频率调制技术测量目标距离的系统。

在频率调制中，电磁波的频率随时间线性增加。换言之，发射频率将以恒定的速率变化，这种频率随时间线性增加的信号称为“啁啾（Chirp）”。FMCW系统测量发射和反射信号频率间的瞬时差 $Δf$ ，也叫做中频信号、差频信号、差拍信号，它与时间差 $Δt$ 成正比，可用于估计目标的距离。

下图(左)显示了啁啾信号的频率与时间的关系，右图显示了频率随时间线性增加的啁啾信号的幅度与时间的关系。啁啾信号的特征在于起始频率( $f _c$ )、带宽( $B$ )和持续时间( $T _c$ )，啁啾信号的斜率定义啁啾信号上升的速率，斜率为 $S = B/T_c$ 。

雷达前方的单个目标会产生一个恒定频率的中频信号（IF signal），频率由下式给出： $\mathrm{IF}_{\text {frequency }}=\frac{2 S d}{c}$ 其中， $d$ 是目标与雷达的距离，单位为 $m$ ， $c$ 是光速，单位为 m/s。因此，我们可以通过测量IF信号经过（DFT）FFT后的频率来估计目标的距离。

本文关于原理部分的介绍十分简单，关于更细节的内容可以阅读之前分享的文章：雷达原理 | 用MATLAB信号处理是如何解算目标的距离和速度信息的？

（二）距离检测Python仿真

这里我用的Python开发平台是Spyder最新版本，属于开源软件。虽然没有Pycharm那么强大，库没有那么多，但对于简单的程序仿真完全足够了，软件安装也简单，没有Python基础的人也就一个星期就学会这个工具了，有了Matlab基础的人直接就可以上手了。（Spyder下载链接：https://www./）

教程以及代码和之前分享的Matlab差不多，详情可见：干货：FMCW雷达系统信号处理建模与仿真（含MATLAB代码）

步骤1：雷达参数设置

设置雷达系统的基本参数。

#Radar parameters setting
maxR = 200 #Maximum range
rangeRes = 1 #Range resolution 
maxV = 70 #Maximum speed
fc = 77e9 #Carrier frequency
c = 3e8 #Speed of light

r0 = 100 #Target distance
v0 = 70 #Target speed

B = c/(2*rangeRes) #Bandwidth
Tchirp = 5.5*2*maxR/c #Chirp time
endle_time = 6.3e-6 
slope = B/Tchirp #Chirp slope
f_IFmax = (slope*2*maxR)/c #Maximum IF frequency
f_IF = (slope*2*r0)/c #Current IF frequency

Nd = 128 #Number of chirp
Nr = 1024 #Numnber ADC sampling points
vres = (c/fc)/(2*Nd*(Tchirp+endle_time)) #Speed resolution 
Fs = Nr/Tchirp #Sampling rate

步骤2：Tx 信号模型

假设Tx信号是频率随时间线性变化的余弦信号。

t = np.linspace(0,Nd*Tchirp,Nr*Nd) #Time of Tx and Rx
angle_freq = fc*t+(slope*t*t)/2 #Tx signal angle speed
freq = fc + slope*t #Tx frequency
Tx = np.cos(2*np.pi*angle_freq) #Waveform of Tx

步骤3：Rx 信号模型

Rx波形可以从Tx波形和延迟时间计算出来。

td = 2*r0/c
freqRx = fc + slope*(t)
Rx = np.cos(2*np.pi*(fc*(t-td) + (slope*(t-td)*(t-td))/2))

步骤4：中频信号模型

根据处理，假设中频信号可以用cos((2*pi*wt*t-2*pi*wr*t))表示。

IF_angle_freq = fc*t+(slope*t*t)/2 - ((fc*(t-td) + (slope*(t-td)*(t-td))/2))
freqIF = slope*td
IFx = np.cos(-(2*np.pi*(fc*(t-td) + (slope*(t-td)*(t-td))/2))+(2*np.pi*angle_freq))

步骤5：中频信号的FFT

在这一步中，我们通过中频信号的 FFT 计算中频信号的频率。

IF_angle_freq = fc*t+(slope*t*t)/2 - ((fc*(t-td) + (slope*(t-td)*(t-td))/2))
freqIF = slope*td
IFx = np.cos(-(2*np.pi*(fc*(t-td) + (slope*(t-td)*(t-td))/2))+(2*np.pi*angle_freq))

步骤6：随时间变化的频谱图

在这一步中，将计算随时间变化的频谱图。

我们可以看到，单帧周期内由于目标位移引起的IF信号频率变化在频谱图中很难区分，因此我们需要通过相位变化来检测微小的位移和速度。

完整的python源代码如下：

import numpy as np
import matplotlib.pyplot as plt
from numpy import fft
from mpl_toolkits.mplot3d import Axes3D

#Radar parameters setting
maxR = 200 
rangeRes = 1
maxV = 70 
fc = 77e9
c = 3e8 
r0 = 100 
v0 = 70

B = c/(2*rangeRes)
Tchirp = 5.5*2*maxR/c
endle_time = 6.3e-6
slope = B/Tchirp 
f_IFmax = (slope*2*maxR)/c
f_IF = (slope*2*r0)/c 

Nd = 128 
Nr = 1024
vres = (c/fc)/(2*Nd*(Tchirp+endle_time)) 
Fs = Nr/Tchirp 
#Tx = np.zeros(1,len(t))
#Rx = np.zeros(1,len(t)) 
#Mix = np.zeros(1,len(t)) 

#Tx波函数参数
t = np.linspace(0,Nd*Tchirp,Nr*Nd) #发射信号和接收信号的采样时间
angle_freq = fc*t+(slope*t*t)/2 #角频率
freq = fc + slope*t #频率
Tx = np.cos(2*np.pi*angle_freq) #发射波形函数
plt.subplot(4,2,1)
plt.plot(t[0:1024],Tx[0:1024])
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.title('Tx Signal')
plt.subplot(4,2,3)
plt.plot(t[0:1024],freq[0:1024])
plt.xlabel('Time')
plt.ylabel('Frequency')
plt.title('Tx F-T')

r0 = r0+v0*t

#Rx波函数参数
td = 2*r0/c
tx = t
freqRx = fc + slope*(t)
Rx = np.cos(2*np.pi*(fc*(t-td) + (slope*(t-td)*(t-td))/2)) #接受波形函数
plt.subplot(4,2,2)
plt.plot(t[0:1024],Rx[0:1024])
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.title('Rx Signal')
plt.subplot(4,2,3)
plt.plot(t[0:1024]+td[0:1024],freqRx[0:1024])
plt.xlabel('Time')
plt.ylabel('Frequency')
plt.title('Chirp F-T')

#IF信号函数参数
IF_angle_freq = fc*t+(slope*t*t)/2 - ((fc*(t-td) + (slope*(t-td)*(t-td))/2))
freqIF = slope*td
IFx = np.cos(-(2*np.pi*(fc*(t-td) + (slope*(t-td)*(t-td))/2))+(2*np.pi*angle_freq))
plt.subplot(4,2,4)
plt.plot(t[0:1024],IFx[0:1024])
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.title('IFx Signal')

#Range FFT
doppler = 10*np.log10(np.abs(np.fft.fft(IFx[0:1024])))
frequency = np.fft.fftfreq(1024, 1/Fs)
range = frequency*c/(2*slope)
plt.subplot(4,2,5)
plt.plot(range[0:512],doppler[0:512])
plt.xlabel('Frequency->Distance')
plt.ylabel('Amplitude')
plt.title('IF Signal FFT')

#2D plot
plt.subplot(4,2,6)
plt.specgram(IFx,1024,Fs)
plt.xlabel('Time')
plt.ylabel('Frequency')
plt.title('Spectogram')

plt.tight_layout(pad=3, w_pad=0.05, h_pad=0.05)
plt.show()

运行后全部的结果如下：

二、速度检测

(一)速度检测原理

微小范围的移动很难通过距离频率关系来检测，如下图所示，在帧期间频谱中没有发现明显的偏移。此外，如果有多个距离相同的目标，我们无法通过距离频率关系来区分它们，因为它们在频谱中具有相同的中频频率。

但是，我们可以通过测量IF信号的相位变化来检测这些微运目标，并区分不同的目标，通过相位变化估计速度的基本思想如下：

(二)速度检测Python仿真

一维速度估计：

步骤1：速维度数据提取

每个啁啾提取一个采样点，对于具有 128 个啁啾的帧，将有 128 个点的列表。

chirpamp = []
chirpnum = 1
while(chirpnum<=Nd):
    strat = (chirpnum-1)*1024
    end = chirpnum*1024
    chirpamp.append(IFx[(chirpnum-1)*1024])
    chirpnum = chirpnum + 1

步骤2：相位差和速度的速度维度 FFT

doppler = 10*np.log10(np.abs(np.fft.fft(chirpamp)))
FFTfrequency = np.fft.fftfreq(Nd,1/Fs)
velocity = 5*np.arange(0,Nd)/3

2D FFT 和速度-距离关系

Z_fft2 = abs(np.fft.fft2(mat2D))
Data_fft2 = Z_fft2[0:64,0:512]

本节完整代码如下，使用时将下面的代码续接在距离检测仿真代码之后，即可运行：

#Range Spectogram
plt.subplot(4,2,6)
plt.specgram(IFx,1024,Fs)
plt.xlabel('Time')
plt.ylabel('Frequency')
plt.title('Range-Frequency Spectogram')
plt.show()

#Speed Calculate
#速度维数据提取
chirpamp = []
chirpnum = 1
while(chirpnum<=Nd):
    strat = (chirpnum-1)*1024
    end = chirpnum*1024
    chirpamp.append(IFx[(chirpnum-1)*1024])
    chirpnum = chirpnum + 1
#速度维做FFT得到相位差
doppler = 10*np.log10(np.abs(np.fft.fft(chirpamp)))
FFTfrequency = np.fft.fftfreq(Nd,1/Fs)
velocity = 5*np.arange(0,Nd)/3
#plt.subplot(4,2,7)
plt.plot(velocity[0:int(Nd/2)],doppler[0:int(Nd/2)])
plt.xlabel('Velocity')
plt.ylabel('Amplitude')
plt.title('IF Velocity FFT')
plt.show()

#2D plot
mat2D=np.zeros((Nd, Nr))
i = 0
while(i<Nd):
    mat2D[i, :] = IFx[i*1024:(i+1)*1024]
    i = i + 1
#plt.matshow(mat2D)
#plt.title('Original data')
#图形绘制，需要修改一下才好看，这里就留个大家自行修改了。
Z_fft2 = abs(np.fft.fft2(mat2D))
Data_fft2 = Z_fft2[0:64,0:512]
#plt.subplot(4,2,8)
plt.imshow(Data_fft2)
plt.xlabel("Range")
plt.ylabel("Velocity")
plt.title('Velocity-Range 2D FFT')

plt.tight_layout(pad=3, w_pad=0.05, h_pad=0.05)
plt.show()

三、学习扩展

上述只是一个简单的示例程序，读者可以参考调皮哥之前分享的Matlab代码自行增加，例如静态杂波滤除、CFAR检测、DOA估计，以及其他的各种雷达信号处理算法仿真等，主动权完全在大家的手中，请按照自己的想法去创造吧。

学习雷达是一个理论结合实践的过程，尤其要重视实践的过程，在实践中的每一个问题及其解决方案，都是十分宝贵的学习经验。不论是做研究还是找工作，都要重视学习知识、提出问题、解决问题、复盘总结等能力。人生在世，学习时间尤其宝贵，望珍惜。

参考资料

【1】https://community./t5/Knowledge-Base-Articles/FMCW-radar-working-principle-simulation-based-on-python-Chapter-1-Distance/ta-p/366803
【2】https://community./t5/Knowledge-Base-Articles/FMCW-radar-working-principle-simulation-based-on-python-Chapter-2-Velocity/ta-p/367214