雷达仿真：FMCW DDMA-MIMO 3D点云获取方法

1.DDMA-MIMO原理

        由于TDMA-MIMO采用不同单天线交替发射信号，没有更好的利用发射天线同时工作的发射资源，导致发射功率低以及损耗大，从而使得TDMA波形只能应用在近距离探测的低功率雷达场景。而DDMA波形则能很好的弥补TDMA上述缺点，与TDMA波形下单天线交替工作不同，DDMA波形下所有发射天线同时发射，但是每个发射天线的信号都加入一个特定的频率偏移(在速度维度体现为一个速度频移)，通过这个人为偏移的频率使得不同发射天线的信号在Doppler域上分离出来。使不同发射天线的信号在Doppler域上能分离开来。此外，由于DDMA相比于TDMA来说拥有更高的发射功率，因此DDMA波形也将成为往后车载毫米波雷达MIMO雷达的主流波形。有关于DDMA更多细节请参考DDMA原理，这里不在累赘。

2.信号建模

        有关于信号建模我在雷达仿真：FMCW TDMA-MIMO 3D点云获取方法已经介绍得非常详细，只需要做如下修改就可。我们基于TDMA信号模型做以下修改即可获得DDMA波形：

1) DDMA是Ntx个发射天线同时发射，因此我们需要在接收端分别叠加每个发射天线产生的回波信号。

2）为了区别不同发射天线发射的信号，每个发射天线施加一个固定的频率频移

上式子中，N通常由发射天线决定，但后续采用Empty Band解调通常满足下式：

式中Nt为发射天线数。在本案例中设置为4 Tx 4 Rx，其中发射天线的分布情况为线性排列，相邻发射天线距离为2λ，通过DDMA解调后，可以等价为1Tx16Rx虚拟阵列。由于固定的频率偏移在中频信号中没法体现出来，为了更直观的体现不同发射天线之间的固定频率频移，将这个频率偏移换算成固定的速度频移为：

 其中为发射脉冲间隔(chirp间隔)。因此DDMA信号模型建模代码如下：

for frameid = 1 : numframe
    for txid = 1 : numtx
        for rxid = 1 : numrx
            for chirpid = 1 : numchirp
                for targetid = 1 : numtarget 
                    for  sampleid = 1 : numsamples
                        R = target(targetid,1) + ((chirpid-1)*T)*(target(targetid,2) + (txid - 1)*lambda/2/Emptyband/T) + (frameid-1)*numchirp*T*(target(targetid,2)+ (txid - 1)*lambda/2/Emptyband/T); %计算目标随速度变化的距离,以忽略快时间对距离的影响
                        delay = 2 * R / parameter.c;
                        fixphase = exp(1j*2*pi*(f0*delay-0.5*slope*delay^2));         %中频信号的固定相位，末尾的相位为剩余视频相位，为了更真实还原毫米波雷达数据，本程中没有将其省略
                        fastsampdata(sampleid) = exp(1j*2*pi*slope*delay*t(sampleid)) * fixphase; %这里采用随机起伏   
                    end
                    azimuthphase = exp(1j*(2*pi/lambda)*sind(target(targetid,3)) * ( (rxid -1) * Detarx  + (txid - 1) * Detatx) );  %目标方位角信息
                    temprawdata(:,chirpid,(txid-1)*numrx + rxid,frameid) = temprawdata(:,chirpid,(txid-1)*numrx + rxid,frameid) + fastsampdata.* azimuthphase;
                end
            end
            %加入高斯噪声
            MixIQ = temprawdata(:,:,(txid-1)*numrx + rxid,frameid);
            xigema = std(MixIQ)/db2mag(snr);
            MixIQ = MixIQ + ((randn(size(MixIQ)) + 1i*randn(size(MixIQ)))).*xigema;
            temprawdata(:,:,(txid-1)*numrx + rxid,frameid) = MixIQ;   
        end
    end
    
    %%  叠加
    for i = 1 : numtx
        for k = 1 : numrx
            rawdata(:,:,i,frameid) = rawdata(:,:,i,frameid) + temprawdata(:,:,i + (k - 1)*numrx,frameid);
        end
    end

end

图 1 信号时域图

3.DDMA-MIMO仿真代码

        仿真环境：window 11 ，Matalb2021b

       4T4R两个空带，距离、速度、角度分别为(150,-12,30)、(100,4.4,-15)、(56, 25,13)的仿真效果为：

图 2 4T4R3个目标仿真效果图

通过Emptyband算法解调，采用OMP-CS算法进行方位角估计（关于OMP-CS算法原理可参看我另一篇文章基于压缩感知(CS)的DOA估计方法-OMP-CS算法），这里为什么选用OMP-CS算法进行到DOA估计，主要原因是通过DDMA波形解调后，比较容易解调出每个发射天线和接收天线的单快拍数据，而数据采样点减少时，由于传统的基于协方差估计(Capon、MUSIC、ESPRIT)的方法失效，因此这里采用OMP-CS算法。如果各位头铁硬要使用上诉传统的方法，可以积累多个帧进行操作。方位角估计效果为：

图 3 目标的方位向

将图2图3中目标的距离、速度、角度进行匹配，得到目标的3D点云图像：

图 4 目标点云图

DDMA仿真的部分代码结构如下：

%%  Author : Poulen
%%  Data :  2023.9.11

%%  内容介绍：DDMA-MIMO demo,该案例用于学习DDMA原理以及接收端解调。DDMA主要原理为：在发射端每个发射天线同时发射，但是相邻发射天线之间加入一个人为的固定相移；
%%  这使得不同的发射天线发射的chirp在doppler域得以区分出来，在接收端按照这种差异解调出来。
%%  缺点：使得最大不模糊速度区间变为原来的1/N(N为发射天线数量)，即最大不模糊速度区间减小了
clear 
close all;
clc;
%%  DDMA波形参数设置
%%  参数设计
parameter = RadarSystemDesign();
maxv = parameter.dopplerindex;
%%  获取原始数据
rawdata = RawDataCapture(parameter);

%%  展示时域信号波形
timedomian = squeeze(rawdata(:,10,1,1));
figure;
subplot(2,1,1);
plot(real(timedomian),'Color',[0 0 0],'LineWidth',1.2);
title('原始信号时域波形(实部)');
ylabel('信号幅度');
xlabel('时域采样点');
xlim([0 384]);
subplot(2,1,2);
plot(imag(timedomian),'Color',[0 0 0],'LineWidth',1.2);
title('原始信号时域波形(实部)');
ylabel('信号幅度');
xlabel('时域采样点');
xlim([0 384]);

clear timedomian;

%% 2DFFT(距离-速度估计)
rangedoppler = RangeDopplerEst(parameter,rawdata);

%% 多通道积累
mulaccumulata = MulAccumulataEst(parameter,rangedoppler);

figure;
mesh(parameter.dopplerindex,parameter.rangeindex(1:end/2-38),db(mulaccumulata(1:end/2-38,:,1)));
xlabel('速度(m/s)');
ylabel('距离(m)');
title('RD 图(4T4R Empty Band = 6)');

%% 解调
[CFAROUT,tagetestindex,targeest] = CA_2DCFAR(parameter,mulaccumulata);     %二维cfar
% figure;
% mesh(parameter.dopplerindex,parameter.rangeindex(1:end/2-38),CFAROUT((1:end/2-38),:,1));
% xlabel('速度(m/s)');
% ylabel('距离(m)');
% title('RD CA-CFAR 图');
% % clear mulaccumulata;

[target,AzimuthData] = EmptyBand(parameter,tagetestindex,mulaccumulata,rangedoppler);

%% doa 估计(OMP-CS算法)
[Azimuspectral,pointcloud] = OMPCS(parameter,target,AzimuthData);

figure;
plot(parameter.scale,abs(Azimuspectral));
xlabel('方位向(°)');
ylabel('幅度');
title('OMPCS 角度谱');
%%  绘制3D点云散点图（在该散点图中仅展示第一帧的图像信息）
figure;
X = pointcloud(1,:);
Y = pointcloud(2,:);
Z = pointcloud(3,:);
S = 50;
C = [0 1 0; 1 0 0; 0.5 0.5 0.5];                                            %注意：这是三个点目标是的RGB原组，若要改变目标数量，该原组需要改变
get = scatter3(X,Y,Z,S,C);
xlim([0 parameter.rangeindex(end/2-38)]);
xlabel('距离(m)');
ylim([parameter.dopplerindex(1) parameter.dopplerindex(end)]);
ylabel('速度(m/s)');
zlim([parameter.scale(1) parameter.scale(end)]);
zlabel('角度(°)');
title('3D点云图');
view(30,35);

总结

        文中难免有笔误之处，望各位海涵。另外，制作不易望多多支持。 

代码提取链：TDMA-DDMA 3D点云仿真代码合集