关于克拉美罗下界(CRLB)-及不同DOA估计算法下的方差(性能)对比博文对应的代码

%% 计算线阵下：同一个角度、不同SNR下的CRLB、以及评估不同DOA估计算法下角度估计的方差 %% %线阵下&单目标下的CRLB公式：var(θ) ≥ 6/( (2*pi)^2*M*(M+1)/(M-1)*(L/λ)^2*cos^2(θ)*SNR) %式中，M为阵元个数，SNR为信噪比，L为阵列孔径大小，θ为目标角度(该角度是目标偏离阵列法向的角度)。 %本次仿真中，使用DBF、capon、music三种方法进行对比。 %author: https://blog.csdn.net/xhblair?spm=1000.2115.3001.5343% clear all; close all; clc; %------构造阵列------% arrayNum = 8; arrayPosition = (0:1:arrayNum-1); %均匀线阵 darray = 0.5; %两两阵元之间的间隔 lambda = 3e8/77e9; Darray = darray*(arrayNum-1)*lambda; %总的阵列长度 %------目标信息------% Target_theta = 5; %------回波信号的参数------% snap = 10; %快拍数 改变这里的快拍数，后面DBF测角时我们只抽取其中一个快拍？ SNR = 15:2:35; %SNR范围 S0 = exp(1j*2*pi*darray*sind(Target_theta)*arrayPosition); %------构造三种DOA估计算法下的必要参数------% thetaScan = linspace(-75,75,512); %------进行蒙特卡罗仿真------% MtklNum = 1000; %每个SNR下进行1000次的试验，并计算方差 CRLB = zeros(1,length(SNR)); VarDbf = zeros(1,length(SNR)); %设定一个装载方差的变量 VarCapon = zeros(1,length(SNR)); VarMusic = zeros(1,length(SNR)); for ii = 1:length(SNR) DbfResult = zeros(1,MtklNum); CaponResult = zeros(1,MtklNum); MusicResult = zeros(1,MtklNum); for jj = 1:MtklNum %得到信号 Amp = randn(snap,length(Target_theta)); %这里是给目标加随机的幅值！！不是给各个通道啊 % Amp = ones(snap,length(Target_theta)); S1 = Amp*S0; S = awgn(S1,SNR(ii),'measured'); %测角 %预处理 R1 = (S'*S./snap); %得到接收信号的协方差矩阵 R2 = inv(S'*S./snap); [EV,D] = eig(R1); %特征值分解,D为由INVR的特征值构成的对角矩阵，EV为其特征向量构成的矩阵 【这里是对协方差矩阵的分解！】 EVA = diag(D)'; %抽取D的对角线元素并转置 [EVA,I] = sort(EVA); %从小到大排序，I对应EVA中元素在原来EVA中的位置。 EVA = fliplr(EVA); %对特征值再从大到小排列 EV = fliplr(EV(:,I)); %对特征值对应的特征向量进行对应的从大到小的排序 EN = EV(:,length(Target_theta)+1:arrayNum); %把噪声子空间拿出来 %对每个角度 for pp = 1:length(thetaScan) a = exp(1j*2*pi*darray*sind(thetaScan(pp)).' *arrayPosition); % DbfResult_Tmp(pp) = abs(S*a'); %DBF DbfResult_Tmp(pp) = abs(S(1,:)*a'); %DBF(如果是多快拍数据，在给矩阵加噪声的过程中会影响相位，从而影响测角，这种情况下方差是很差的) CaponResult_Tmp(pp) = 1/abs(a*R2*a'); %capon % MusicResult_Tmp(pp) = abs((a*a')/(a*EN*EN'*a')); %Music MusicResult_Tmp(pp) = abs(1/(a*EN*EN'*a')); %Music end % figure(100); % subplot(131);plot(thetaScan,db(DbfResult_Tmp./max(DbfResult_Tmp)));title('dbf');grid on; % subplot(132);plot(thetaScan,db(CaponResult_Tmp./max(CaponResult_Tmp)));title('capon');grid on; % subplot(133);plot(thetaScan,db(MusicResult_Tmp./max(MusicResult_Tmp)));title('music');grid on; %获得各方法下的角度测量结果：（直接基于取最大的方法来得到测角的结果,这里还只是索引值） [~,DbfResult(jj)] = max(DbfResult_Tmp); [~,CaponResult(jj)] = max(CaponResult_Tmp); [~,MusicResult(jj)] = max(MusicResult_Tmp); end %计算该SNR下的CRLB，当然，也可以在外边一次性求 (需要把SNR从dB转成幅值) CRLB(ii) = 6 / ( (2*pi)^2 * arrayNum * (arrayNum+1)/(arrayNum-1) * (Darray/lambda)^2 * (cosd(Target_theta))^2 * 10^(SNR(ii)/10) ); %计算各方法下，在该SNR下的方差： %由索引值到角度值 DbfResult = (thetaScan(end) - thetaScan(1))/(length(thetaScan)-1).*(DbfResult-1) + thetaScan(1); CaponResult = (thetaScan(end) - thetaScan(1))/(length(thetaScan)-1).*(CaponResult-1) + thetaScan(1); MusicResult = (thetaScan(end) - thetaScan(1))/(length(thetaScan)-1).*(MusicResult-1) + thetaScan(1); for kk = 1:MtklNum %计算方差 [因为是无偏估计，这里求方差应该要用目标的实际值，而不是测量的均值？] % VarDbf(ii) = (DbfResult(kk) - sum(DbfResult)/MtklNum)^2 + VarDbf(ii); % VarCapon(ii) = (CaponResult(kk) - sum(CaponResult)/MtklNum)^2 + VarCapon(ii); % VarMusic(ii) = (MusicResult(kk) - sum(MusicResult)/MtklNum)^2 + VarMusic(ii); VarDbf(ii) = (DbfResult(kk) - Target_theta)^2 + VarDbf(ii); VarCapon(ii) = (CaponResult(kk) - Target_theta)^2 + VarCapon(ii); VarMusic(ii) = (MusicResult(kk) - Target_theta)^2 + VarMusic(ii); end VarDbf(ii) = VarDbf(ii) /MtklNum; VarCapon(ii) = VarCapon(ii)/MtklNum; VarMusic(ii) = VarMusic(ii)/MtklNum; currentSNR = SNR(ii) end % figure(100);hold off; %------画图------% figure(1); subplot(121); plot(SNR,CRLB,'r-p');hold on; plot(SNR,VarDbf,'b-*'); plot(SNR,VarMusic,'k-d'); % plot(SNR,VarCapon,'g-s'); hold off;grid on;legend('CRLB','DBF测角方差','Music测角方差','Capon测角方差'); xlabel('SNR/dB');ylabel('方差值');title('线阵-单目标-不同SNR下的各测角方法的性能(方差)对比'); subplot(122); plot(SNR,db(CRLB),'r-p');hold on; plot(SNR,db(VarDbf),'b-*'); plot(SNR,db(VarMusic),'k-d'); % plot(SNR,db(VarCapon),'g-s'); hold off;grid on;legend('CRLB','DBF测角方差','Music测角方差','Capon测角方差'); xlabel('SNR/dB');ylabel('方差值(dB)');title('线阵-单目标-不同SNR下的各测角方法的性能(方差)对比');