关于雷达仿真中角度维左右反转的问题解决

最新推荐文章于 2024-09-13 23:25:03 发布

爱莉希雅9527

最新推荐文章于 2024-09-13 23:25:03 发布

阅读量829

点赞数 24

文章标签： matlab 信号处理 bug

本文链接：https://blog.csdn.net/m0_60497281/article/details/140700907

版权

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

与一般先根据公式计算回波信号再老老实实混频的仿真方式相比，该仿真直接进行相位的计算，设 $t$ 时刻的发射信号相位为：

$\phi (t)=2\pi (f_{c}t+\frac{u}{2}t^{2})\textup{(1)}$

则该时刻的混频信号为：

$s (t)=e^{j2\pi *(f_{c}t+\frac{u}{2}t^{2}-f_{c}(t-\tau )-\frac{u}{2}(t-\tau )^{2})}\\=e^{j2\pi *(u\tau t+f\tau -\frac{u}{2}\tau ^{2})}\\\approx e^{j2\pi *(u\tau t+f\tau)}\textup{(2)}$

与计算回波信号再混频没有差别。值得一提的是，此处假设了每个chirp内目标位置不变，即同个chirp内时延 $\tau$ 相同。时延 $\tau$ 则是针对每一组发射天线和接收天线组合分别求直线距离，也是相对精确的做法，那么问题究竟出在了哪里？我们再来分析一下雷达测角原理。

关于雷达测角原理的讲解，目前网上流传最广的是这一张图：

该图以向左为角度的正方向，而接收雷达向右排开；而本仿真中向右为角度的正方向，阵元也是向右排开的，在这种情况下，左侧阵元接收的信号应当比右侧阵元多出 $d\sin \theta$ 的路程（如果目标在原点左侧， $d\sin \theta< 0$ ，这一句仍然成立），沿着正方向看，右侧阵元相较于左侧的相位差为：

$\Delta \phi =-2\pi \frac{d\sin \theta }{\lambda }\textup{(3)}$

换算得到：

$\theta =-\arcsin \frac{\lambda \sin \Delta \phi }{2\pi d}\textup{(4)}$

这样看来，角度维FFT解出 $\Delta \phi$ 之后再计算 $\theta$ ，结果似乎就应该是相差一个负号，即左右反转的，没有问题。这两个公式在CSDN的其他博客上都是没有负号的，大概就在于上面提到的方向问题，当然我也不清楚我的理解是否有问题，还有待进一步学习。

【雷达仿真 | FMCW TDMA-MIMO毫米波雷达信号处理仿真（可修改为DDMA-MIMO）】_利用所述曲线斜率得到的目标速度估计值v-r、2dfft的处理得到的模糊速度估计值v-a-CSDN博客中采用了另一种略显正常的回波信号仿真方式，先计算回波信号Sr再计算混频信号Sif。其中采用了近似的手法，即只计算一个阵元的相位，再利用相邻阵元的相位差 $2\pi \frac{d\sin \theta }{\lambda }$ 计算其他阵元：

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）对应。

正的，负的，个中弯弯绕绕实在叫人头疼，这个问题困扰了本人数天才想出了一个较好的解释，写这篇博客又花了本人半天，如果有更好的解释，欢迎讨论；如果略有帮助，还望各位佬不要吝惜自己的赞吧。