本文首發于公眾號【調皮連續波】，其他平臺為自動同步，內容若不全或亂碼，請前往公眾號閱讀。保持關注調皮哥，和1.5W雷達er一起學習雷達技術！

2023年度會員內容更新公告（05.04）

序號	類別	內容	文件路徑
1	雷達代碼	本文內容	根目錄雷達代碼庫

【正文】

編輯|雷達小助理 審核|調皮哥

---

代碼運行平臺：MATLAB2022a

操作系統：Windows 10 專業版

文件目錄：

1、參數設置

雷達的參數按照原始數據采集的參數來設置，本文的MATLAB仿真中設置的雷達參數如下所示，同時設定了距離維和速度維的坐標軸。

%% 參數設置n_chirps=128;                   %每幀脈沖數n_samples=256;                  %采樣點數/脈沖N=256;                          %距離向FFT點數M=128;                          %多普勒向FFT點數
fs=10e6;                        %采樣率c=3.0e8;                        %采樣率f0=77e9;                        %初始頻率lambda=c/f0;                    %雷達信號波長d=lambda/2;                     %天線陣列間距
Tc=50e-6;                       %單chirp周期Tf=40e-3;                       %幀周期B=1344.19e6;                    %信號帶寬HzrangeRes=c/2/B*Tc/2*fs/N;       %距離刻度Ts=0127*Tc;                 %快時間軸Ta=0255*Tf;                 %總時間軸
Rr=0(N-1)*rangeRes;   %距離軸Vr=(-M/2:M/2-1)*lambda/Tc/M/2;  %速度軸

2、數據讀取

數據讀取采用TI官方提供的函數實現，在本文的仿真程序中采用readDCA1443()函數，得到4個接收通道的數據。

%% 數據讀取fname='adc_data.bin'; %文件路徑名稱adcdata =readDCA1443(fname);n_frame=floor(length(adcdata)/n_chirps/n_samples); %數據總幀數
data_rx1= reshape(adcdata(1,:),n_samples,n_chirps,n_frame); %RX1數據data_rx2= reshape(adcdata(2,:),n_samples,n_chirps,n_frame); %RX2數據data_rx3= reshape(adcdata(3,:),n_samples,n_chirps,n_frame); %RX3數據data_rx4= reshape(adcdata(4,:),n_samples,n_chirps,n_frame); %RX4數據
%微多普勒數據profile=zeros(n_frame,n_chirps);

3、距離速度譜

利用距離維FFT和速度維FFT實現，其中在距離維FFT之前采用相量均值相消去除了零頻，窗函數選擇海明窗。仿真結果如下所示：

4、非相參積累

疊加四個通道的信號幅值，如下所示：

 rx_2dfft=abs(rx1_2dfft)+abs(rx2_2dfft)+abs(rx3_2dfft)+abs(rx4_2dfft);

結果如下所示：

5、二維CFAR

CFAR參數根據雷達的分辨率等參數設定，本文仿真的參數設置如下所示，其中，虛警概率pfa =10^-6。

%通過完整的距離多普勒圖滑動窗口%在兩個維度中選擇參考單元的數量Tr = 8; Td = 4;%選擇被測單元（CUT）周圍兩個維度的保護單元數量，以進行準確CFAR檢測Gr = 4;Gd = 2;

結果如下所示：

2D CFAR的算法執行效率較低，讀者可以根據需求選擇先進行速度維CFAR，然后再進行距離維CFAR。

6、峰值搜索

首先建立檢測矩陣，擴展補零，便于檢測邊界目標。

rx_2dcfar_temp=padarray(rx_2dcfar,[1,1],0);   %建立檢測矩陣，擴展補零，檢測邊界目標[r,c]=size(rx_2dcfar);                        %矩陣大小

建立3*3的滑窗，通過找出滑窗內最大的值及其坐標，用于后續提取峰值點。

for i=1:r                                     %創建3*3滑窗，找出較滑窗內數據最大值，剩余用0覆蓋    for j=1:c        if  rx_2dcfar_temp(i+1,j+1)>0                               a=rx_2dcfar_temp(i:i+2,j:j+2);    %創建3*3滑窗            b=max(max(a));                    %找出較滑窗內數據最大值            [x,y]=find(a==max(max(a)));       %找出較滑窗內數據最大值在滑窗內坐標            temp=zeros(3,3);                  %創建3*3全0矩陣            rx_2dcfar_temp(i:i+2,j:j+2)=temp; %用3*3全0矩陣覆蓋檢測矩陣內數據，代表檢測矩陣數據檢測完成            temp(x,y)=b;                      %創建3*3全0保留最大值的矩陣            rx_2dcfar_temp(i:i+2,j:j+2)=temp; %用上述矩陣覆蓋數據矩陣內數據，        end    endend
rx_2dcfar_plots=rx_2dcfar_temp(2:r+1,2:c+1);

效果如下所示：

7、微多普勒

微多普勒提取沒有采用STFT法，而是直接提取RD譜中的低速部分速度，作為微多普勒。這種方法比較簡單，也能在一定程度上近似微多普勒，工程上比較實用。

但是這種方法需要精細確定目標檢測點的距離門分布范圍，否則會漏掉一些目標點。例如，下面的代碼中，j代表距離門范圍。直接將每一幀的目標速度疊加起來，最后就能夠得到近似的微多普勒譜。

%微動多普勒forj=3:6    
   profile(n,:)=profile(n,:)+ rx_2dcfar(j,:);
end

仿真結果如下圖所示：

8、角度估計

采用角度維FFT，角度估計的其他方法，比如DBF、Capon、MUSIC、壓縮感知以及其他超分辨算法等，需要讀者自行探索，公眾號以往的一些文章也提到過一些，可以參考。

9、微多普勒速度時域包絡

利用微多普勒頻率fd乘以速度v ，可以得到fd*v =2(v^2)/λ，即得到了以2/λ為幅度系數的微多普勒速度包絡v^2,也就是微多普勒速度的功率，如下圖所示。

包絡檢波得到的是帶通信號的基帶部分，在本文中也就是微多普勒速度的時域信號。

10、微多普勒調制頻譜圖

微多普勒速度頻域可以利用對微多普勒速度時域信號做FFT求得，而實信號FFT得到雙邊譜，如下圖所示。

% 微動多普勒繪圖profile=profile';figure(8)colormap(jet);imagesc(Ta,Vr,(profile));title('微動多普勒圖');xlabel('時間 s');ylabel('速度m/s');
profile_x=Vr*profile;%微動多普勒與速度乘積累加figure(9)plot(Ta,profile_x);xlabel('時間 s');ylabel('幅度');title("微多普勒時域包絡圖")
profile_x=fftshift(fft(fftshift(profile_x,256)));  %FFTFs=1/(Tf*256);F=(-128:127)*Fs; %頻率軸figure(10)plot(F,(abs(profile_x)));title("微多普勒頻域圖");xlabel("頻率 Hz");ylabel("幅值");

這個譜圖的峰值表示，微多普勒頻率，可以看到上圖，峰值較強的有4個。

本文到此結束，年度會員直接查看公告目錄，非會員可私信博主。

【點擊以下鏈接可直達各個業務模塊】

加入雷達群	加入年度會員
雷達項目交流	付費咨詢
商業推廣合作	文章投稿指南

【本期結束】

---

本文是空閑時個人的心得體會，僅供參考。目前我還有很多內容需要學習，如果還有沒有說到或者不全面的地方，還請指正，感謝大家。

喜歡本文，可以轉發朋友圈。歡迎關注【調皮連續波】和備用號【跳頻連續波】。