https://github.com/ekurtgl/FMCW-MIMO-Radar-Simulation是近几年比较新的雷达信号处理开源MATLAB仿真代码,涉及FMCW雷达与MIMO阵列体制的3d FFT、CA-CFAR、MUSIC算法、压缩感知等。在网络上也有一些分析该代码的博客,例如调皮连续波的开源代码 | FMCW-MIMO雷达仿真MATLAB_matlab 分析 mimo雷达分辨率-CSDN博客,还是很值得学习的。
有意思的是,本人在学习时发现,截止24年6月份该代码的版本里,都存在角度维FFT与MUSIC算法结果相反的问题,在一众分析文章中也没有发现该异常。两种算法得到的距离-角度谱罗列如下:
角度维FFT结果
MUSIC解角结果
目标初始位置在(50m,50°)(100m,10°)。可以看到,两种算法在角度维上的结果左右相反了。这是由什么导致的?本人详细分析了代码,最终找到了原因,现将原作者的部分代码放在下面,避免有人当时空警察。
首先来看一下角度维FFT部分的代码:
rangeFFT = fft(RDC(:,1:numChirps,:),N_range);
angleFFT = fftshift(fft(rangeFFT,length(ang_ax),3),3);
range_az = squeeze(sum(angleFFT,2)); % range-azimuth map
其中,RDC为雷达数据立方体,三个维度分别为距离维、多普勒维、角度维,从中取出1帧来成像;rangeFFT为距离维FFT结果,angleFFT为角度维FFT结果,角度维FFT长度为角度维坐标数量,且使用了fftshift变换,看起来似乎没有问题,可角度维的结果已经与目标角度相反了,这是为什么?再来看看作者标新立异的混频信号生成了:
for i = 1:numTX
for j = 1:numRX
delays_tar1{i,j} = (vecnorm(tar1_loc-repmat(rx_loc{j},N,1),2,2)+vecnorm(tar1_loc-repmat(tx_loc{i},N,1),2,2))/c;
delays_tar2{i,j} = (vecnorm(tar2_loc-repmat(rx_loc{j},N,1),2,2)+vecnorm(tar2_loc-repmat(tx_loc{i},N,1),2,2))/c;
end
end
phase = @(tx,fx) 2*pi*(fx.*tx+slope/2*tx.^2); % transmitted
mixed = cell(numTX,numRX);
for i = 1:numTX
for j = 1:numRX
disp(['Processing Channel: ' num2str(j) '/' num2str(numRX)]);
for k = 1:numChirps*numCPI
phase_t = phase(t_onePulse,fc);
phase_1 = phase(t_onePulse-delays_tar1{i,j}(k*numADC),fc); % received
phase_2 = phase(t_onePulse-delays_tar2{i,j}(k*numADC),fc);
signal_t((k-1)*numADC+1:k*numADC) = exp(1j*phase_t);
signal_1((k-1)*numADC+1:k*numADC) = exp(1j*(phase_t - phase_1));
signal_2((k-1)*numADC+1:k*numADC) = exp(1j*(phase_t - phase_2));
end
mixed{i,j} = signal_1 + signal_2;
end
end
与一般先根据公式计算回波信号再老老实实混频的仿真方式相比,该仿真直接进行相位的计算,设时刻的发射信号相位为:
则该时刻的混频信号为:
与计算回波信号再混频没有差别。值得一提的是,此处假设了每个chirp内目标位置不变,即同个chirp内时延相同。时延
则是针对每一组发射天线和接收天线组合分别求直线距离,也是相对精确的做法,那么问题究竟出在了哪里?我们再来分析一下雷达测角原理。
关于雷达测角原理的讲解,目前网上流传最广的是这一张图:
该图以向左为角度的正方向,而接收雷达向右排开;而本仿真中向右为角度的正方向,阵元也是向右排开的,在这种情况下,左侧阵元接收的信号应当比右侧阵元多出的路程(如果目标在原点左侧,
,这一句仍然成立),沿着正方向看,右侧阵元相较于左侧的相位差为:
换算得到:
这样看来,角度维FFT解出之后再计算
,结果似乎就应该是相差一个负号,即左右反转的,没有问题。这两个公式在CSDN的其他博客上都是没有负号的,大概就在于上面提到的方向问题,当然我也不清楚我的理解是否有问题,还有待进一步学习。
【雷达仿真 | FMCW TDMA-MIMO毫米波雷达信号处理仿真(可修改为DDMA-MIMO)】_利用所述曲线斜率得到的目标速度估计值v-r、2dfft的处理得到的模糊速度估计值v-a-CSDN博客中采用了另一种略显正常的回波信号仿真方式,先计算回波信号Sr再计算混频信号Sif。其中采用了近似的手法,即只计算一个阵元的相位,再利用相邻阵元的相位差计算其他阵元:
t = 0:1/fs:Tr-(1/fs); %chirp采样的时间序列
for chirpId = 1:chirps
for txId = 1:txNum
St = exp((1i*2*pi)*(centerFreq*(t+( chirpId-1)*Tr)+slope/2*t.^2)); %发射信号
for rxId = 1:rxNum
Sif = zeros(1,rangeBin);
for targetId = 1:targetNum
%%连续帧 目标设置,如果不需要连续帧,令Parameter.frame=0,即可。
if targetId==1
targetRange = target(targetId,1)-Parameter.frame;
targetSpeed = target(targetId,2);
targetAngle = target(targetId,3);
elseif targetId==2
targetRange = target(targetId,1)+0.5*Parameter.frame;
targetSpeed = target(targetId,2);
targetAngle = target(targetId,3);
elseif targetId==3
targetRange = target(targetId,1)+Parameter.frame;
targetSpeed = target(targetId,2);
targetAngle = target(targetId,3);
end
tau = 2 * (targetRange + targetSpeed * (txId - 1) * Tr) / c;
fd = 2 * targetSpeed / lambda;
wx = ((txId-1) * rxNum + rxId) / lambda * dx * sind(targetAngle);
Sr = 10*exp((1i*2*pi)*((centerFreq-fd)*(t-tau+( chirpId-1) * Tr)+slope/2*(t-tau).^2 + wx)); %回波信号
Sif = Sif + St .* conj(Sr);
%叠加20dB高斯白噪声
Sif = awgn(Sif,20);
end
rawData((txId-1) * rxNum + rxId,:,chirpId) = Sif;
end
end
end
由于它的相位差是正的,按照上文的分析就是阵元方向与角度正方向相反了,但是下面用复共轭信号进行混频(仿真使用复共轭是为了结果更清晰),相位差其实也加了个负号,所以阵元方向还是与角度正方向相同的,那它的角度维FFT结果怕不是也是反的了。
还有一点要说明的是MUSIC算法的结果怎么不是反的呢?这要归功于它的方向矢量。下面是开源代码中的MUSIC算法:
d = 0.5;
M = numCPI; % # of snapshots
for k=1:length(ang_ax)
a1(:,k)=exp(-1i*2*pi*(d*(0:numTX*numRX-1)'*sin(ang_ax(k).'*pi/180)));
end
rangeFFT = fft(RDC);
for i = 1:N_range
Rxx = zeros(numTX*numRX,numTX*numRX);
for m = 1:M
A = squeeze(sum(rangeFFT(i,(m-1)*numChirps+1:m*numChirps,:),2));
Rxx = Rxx + 1/M * (A*A');
end
% Rxx = Rxx + sqrt(noise_pow/2)*(randn(size(Rxx))+1j*randn(size(Rxx)));
[Q,D] = eig(Rxx); % Q: eigenvectors (columns), D: eigenvalues
[D, I] = sort(diag(D),'descend');
Q = Q(:,I); % Sort the eigenvectors to put signal eigenvectors first
Qs = Q(:,1); % Get the signal eigenvectors
Qn = Q(:,2:end); % Get the noise eigenvectors
for k=1:length(ang_ax)
music_spectrum2(k)=(a1(:,k)'*a1(:,k))/(a1(:,k)'*(Qn*Qn')*a1(:,k));
end
range_az_music(i,:) = music_spectrum2;
end
方向矢量a1代表信号从某个方向到达不同天线时产生的相位差。仍然是角度正方向和阵列延伸方向都朝右的情况,经过测试,该算法中方向矢量a1如果相位上添加负号,那么结果就是正确的;如果改为正的,结果就会相反。换言之,MUSIC算法的方向矢量似乎要与公式(3)对应。
正的,负的,个中弯弯绕绕实在叫人头疼 ,这个问题困扰了本人数天才想出了一个较好的解释,写这篇博客又花了本人半天,如果有更好的解释,欢迎讨论;如果略有帮助,还望各位佬不要吝惜自己的赞吧。