基于干扰协方差矩阵重构的鲁棒波束形成算法.zip

%% 功能描述：存在DOA偏差情况下对比不同（鲁棒）波束形成方法形成的波束图 %% -- 开始编码 --------------------------------------------------- clc; clear; % close all; %% -- 信号参数设置 ------------------------------------------------ f_SigThetaPi = 80*pi/180; % 信号方位角度 f_SNRdB = 0; % 信噪比 vf_IntThetaPi = [60, 105]*pi/180; % 干扰方位角度 vf_INRdB = [30, 40]; % 干噪比 n_IntNum = length(vf_IntThetaPi); % 加入的干扰个数 f_NoisePower = 1; % 噪声功率为 1 %% -- 阵列参数设置 ------------------------------------------------ n_EleNum = 16; % 阵元个数 n_SnapNum = 50; % 快拍数目 f_DToLamda = 0.5; % 归一化波长（阵元间距比波长，半波长间距） vf_ElePosLamd = [0:n_EleNum-1]'*f_DToLamda; % 阵元归一化位置序列 vf_SigSteerVec = exp(-1j*2*pi*vf_ElePosLamd*cos(f_SigThetaPi)); % 信号导向矢量 af_IntSteerVec = exp(-1j*2*pi*vf_ElePosLamd*cos(vf_IntThetaPi));% 信号导向矢量 f_DOAError = 2/180*pi; f_DesireThetaPi = f_SigThetaPi-f_DOAError; % 波束形成期望角度 vf_DesireSteerVec = exp(-1j*2*pi*vf_ElePosLamd*cos(f_DesireThetaPi)); %% -- 信号、干扰、噪声产生 ----------------------------------------- %---------产生噪声------------------------- f_NoiseAmpl = sqrt(f_NoisePower/2); % 噪声幅度 af_RcvNoise = f_NoiseAmpl*(randn(n_EleNum, n_SnapNum)+1j*randn(n_EleNum, n_SnapNum)); %---------产生信号------------------------- f_SigPower = f_NoisePower*10^(f_SNRdB./10); % 信号功率 f_SigAmpl = sqrt(f_SigPower/2); % 信号幅度 % vf_SigData = randi([0,1],[n_SnapNum,1])*2-1; vf_SigData = f_SigAmpl*(randn(n_SnapNum,1)+1j*randn(n_SnapNum,1)); af_RcvSig = vf_SigSteerVec*vf_SigData.'; % 阵列接收的信号 %---------产生干扰-------------------------- vf_IntAmpl = sqrt((f_NoisePower*10.^(vf_INRdB./10))/2); % 干扰幅度 vf_IntData = kron(vf_IntAmpl,ones(n_SnapNum, 1)).*(randn(n_SnapNum,n_IntNum)+1j*randn(n_SnapNum,n_IntNum)); af_RcvInt = af_IntSteerVec*vf_IntData.'; % 阵列接收的干扰 %---------阵列接收数据---------------------- af_RcvDataAll = af_RcvNoise+af_RcvSig+af_RcvInt; % 阵列接收数据总和 %---------计算协方差矩阵逆------------------ af_SmiR = af_RcvDataAll*af_RcvDataAll'/n_SnapNum; % 协方差矩阵 af_InvSmiR = inv(af_SmiR); % 协方差矩阵求逆 af_Rcv_i_n = af_RcvNoise+af_RcvInt; af_R_IN = af_Rcv_i_n*af_Rcv_i_n'/n_SnapNum; % 干扰噪声协方差矩阵 af_InvR_IN = inv(af_R_IN); % 协方差矩阵求逆 %% -- 空间能量分布计算 --------------------------------------------- vf_SearchThetaPi = [40:0.5:140].'/180*pi; % 搜索角度序列 n_ThetaSNum = length(vf_SearchThetaPi); % 搜索角度个数 vf_SpaceSpec = zeros(n_ThetaSNum, 1); % 空间谱存储变量 for i_loop = 1:n_ThetaSNum vf_SearchSteerVec = exp(-1j*2*pi*vf_ElePosLamd*cos(vf_SearchThetaPi(i_loop))); % 扫描导向矢量 vf_SpaceSpec(i_loop) = 1./abs(vf_SearchSteerVec'*af_InvSmiR*vf_SearchSteerVec); end figure();plot(vf_SearchThetaPi/pi*180, 10*log10(vf_SpaceSpec)); grid on; xlabel('Angle (^o)'); ylabel('Spatial Spectrum (dB)') %% -- 方向图计算 -------------------------------------------------- vf_BeamPatternCBF = zeros(n_ThetaSNum, 1); % CBF波束图存储变量 vf_BeamPatternSMI = zeros(n_ThetaSNum, 1); % Capon（采样协方差矩阵）波束图存储变量 vf_BeamPatternLSMI = zeros(n_ThetaSNum, 1); % 对角加载Capon（采样协方差矩阵）波束图存储变量 vf_BeamPatternEIG = zeros(n_ThetaSNum, 1); % 特征子空间鲁棒波束图存储变量 vf_BeamPatternWorst = zeros(n_ThetaSNum, 1); % 最差性能波束图存储变量 vf_BeamPatternRFix = zeros(n_ThetaSNum, 1); % 协方差矩阵重构束图存储变量 vf_BeamPatternOpt = zeros(n_ThetaSNum, 1); % 最优（干扰噪声协方差矩阵）波束图存储变量 f_DLdB = 10; % 对角加载量，噪声功率的10dB af_InvSmiR_Dl = inv(af_SmiR + 10^(f_DLdB/10)*eye(n_EleNum));% 协方差矩阵求逆 %---------CBF波束形成权向量计算-------------- vf_CBFWeight = vf_DesireSteerVec/n_EleNum; % CBF 波束形成权向量 %---------SMI波束形成权向量计算-------------- vf_SMIWeight = af_InvSmiR*vf_DesireSteerVec/(vf_DesireSteerVec'*af_InvSmiR*vf_DesireSteerVec); % SMI-Capon波束形成权向量 %---------LMI波束形成权向量计算-------------- vf_LSMIWeight = af_InvSmiR_Dl*vf_DesireSteerVec/(vf_DesireSteerVec'*af_InvSmiR_Dl*vf_DesireSteerVec); % SMI-Dl-Capon波束形成权向量 vf_LSMIWeight = vf_LSMIWeight./norm(vf_LSMIWeight)/sqrt(n_EleNum); %---------最优波束形成权向量计算------------- vf_OptWeight = af_InvR_IN*vf_DesireSteerVec/(vf_DesireSteerVec'*af_InvR_IN*vf_DesireSteerVec); % SMI-Capon波束形成权向量 %---------特征子空间方法权向量计算------------- [V,S,U] = svd(af_SmiR); % 特征值分解，R=VSU'; E = V(:,1:3); % D = S(1:3,1:3); % vf_EigWeight = E*inv(D)*E'*vf_DesireSteerVec; % vf_EigWeight = vf_EigWeight./norm(vf_EigWeight)/sqrt(n_EleNum); %---------最差性能波束形成权向量计算------------- epsilon = 2; cvx_quiet(true) % cvx_begin % variable W(n_EleNum) complex % minimize real(W'*af_SmiR*W) % subject to % 1+epsilon*norm(W) <= real(W'*vf_DesireSteerVec); % imag(W'*vf_DesireSteerVec) == 0; % cvx_end % vf_WorstWeight = W./norm(W)/sqrt(n_EleNum); % get weight by CVX，并归一化 %---------协方差矩阵重构波束形成权向量计算------------- % 重构协方差矩阵 f_DeltaTheta = 5*pi/180; af_RIN_Fix = zeros(n_EleNum, n_EleNum); for i=1:n_ThetaSNum if abs(vf_SearchThetaPi(i)-f_DesireThetaPi) > f_DeltaTheta; vf_SearchSteerVec = exp(-1j*2*pi*vf_ElePosLamd*cos(vf_SearchThetaPi(i)));% 扫描导向矢量 af_RIN_Fix = af_RIN_Fix + vf_SearchSteerVec*vf_SearchSteerVec'... /abs(vf_SearchSteerVec'*af_InvSmiR_Dl*vf_SearchSteerVec); end end % 协方差矩阵加权 af_InvRIN_Fix = pinv(af_RIN_Fix); cvx_quiet(true) % cvx_begin % variable E_V(n_EleNum) complex % minimize quad_form(vf_DesireSteerVec+E_V, af_InvRIN_Fix) % subject to % quad_form(vf_DesireSteerVec+E_V, af_RIN_Fix) <= quad_form(vf_DesireSteerVec, af_RIN_Fix) E_V'*vf_DesireSteerVec == 0; % cvx_end % vf_SteerFix = sqrt(n_EleNum)*(vf_DesireSteerVec+E_V)/norm(vf_DesireSteerVec+E_V); % get weight by CVX vf_RFixWeight = af_InvRIN_Fix*vf_SteerFix/abs(vf_SteerFix'*af_InvRIN_Fix*vf_SteerFix); % get weight by CVX vf_RFixWeight = vf_RFixWeight./norm(vf_RFixWeight)/sqrt(n_EleNum); %---------方向图计算（各个角度扫描）---------- for i_loop = 1:n_ThetaSNum f_ThetaSNow = vf_SearchThetaPi(i_loop); % 当前扫描角度 vf_SearchSteerVec = exp(-1j*2*pi*vf_ElePosLamd*cos(f_ThetaSNow));% 扫描导向矢量 vf_BeamPatternCBF(i_loop) = abs(vf_CBFWeight'*vf_SearchSteerVec).^2; vf_BeamPatternSMI(i_loop) = abs(vf_SMIW