On Robust Capon Beamformingand Diagonal Loading

"On Robust Capon Beamformingand Diagonal Loading" 是阵列信号处理领域被引用很广泛的一篇论文，这里贴出这篇论文的MATLAB代码。

发这篇博文的目的是希望感兴趣以及有问题的网友在评论中留言交流。

% FUNCTION [a_est, phiRCB, phiB, phiC, thetaRange ] = RCB_JianLi( R, L, epsilon, d )
% 
% The function calculates the robust Capon beamformer. This implementation 
% follows (please cite if used): 
%
% J. Li, P. Stoica, and Z. Wang, "On Robust Capon Beamforming and Diagonal
% Loading", IEEE Transactions on Signal Processing, Vol 51, No. 7, pp.
% 1702-1715, July 2003. 
%
% Parameters:
%   R           Input data,sample covariance matrix
%   L           Number of grid points to evaluate over.
%   epsilon     Size of the allowed uncertainty sphere (default 1.0).
%   d           distance of antenna (default 0.5 lambda)
%   a_est       Corrected steering vector 
%   PRCB        The power of robust Capon beamformer.
%   phiB        The power of classical beamformer.
%   PCapon      The power of Capon beamformer.
%   thetaRange  The considered frequency grid. 
%
% By Zhisong Wang. Last modified by Andreas Jakobsson, 160104.
%
function [a_est, PRCB, phiB, PCapon, thetaRange ] = rcb_JianLi( R, L, epsilon, d )
   
    if nargin < 3
        epsilon = 1;
    end
    thetaRange = linspace( -90, 90, L );
    
    if nargin < 4
        d = 0.5;
    end

    [U,Gamma,V] = svd(R); %对采样协方差矩阵进行特征分解
    Gamma       = real(Gamma);  % Gamma是一个对角矩阵，特征值实部
    Gamma_vec   = diag(Gamma);  % 取出Gamma的对角，得到一个向量
    lambda_max  = max(diag(Gamma));
    lambda_min  = min(diag(Gamma));
    num = 1; 
    m   = length(R);
    phiB = zeros(L,1);    
    PCapon = phiB;    
    PRCB = phiB; 

    for assumed_theta=thetaRange
        a_bar=exp(-2*pi*1j*d*sin(assumed_theta*pi/180)*[0:m-1].');       % Assumed steering vector (without calibration errors).
        
%         PCapon(num)=real(1/(a_bar'*inv(R)*a_bar));                    % Capon estimate.
        Delay_sum_weight=a_bar/m;
        phiB(num)=real(Delay_sum_weight'*R*Delay_sum_weight);       % Classical beamformer estimate.

        % Compute RCB using Newton's method 
        z=U'*a_bar;
        
        error=10;       % Initially set a large number for the error

        % Lower bound for the optimal solution of lambda
        LowerBound=(norm(a_bar)/sqrt(epsilon)-1)/lambda_max;
        disp(norm(a_bar));
        disp(LowerBound);

        lambda_n=LowerBound;	

        while (abs(error)>1e-6) 
                 [lambda_m,error]=RCB_Newton(lambda_n,z,Gamma_vec,m,epsilon);
                 lambda_n=lambda_m;	
        end
        lambda_Newton=lambda_n;     % Optimal solution for lambda
        PCapon(num)=real(1/(a_bar'*inv(R+1/lambda_Newton*eye(8))*a_bar));
        
        for qq=1:m 
          gamma_i=Gamma(qq,qq);
          QQ(qq,qq)=1/(1/lambda_Newton^2 + 2/lambda_Newton*gamma_i + gamma_i^2)*gamma_i;
          QQ_scale(qq,qq)=QQ(qq,qq)*gamma_i;
        end	
        a_est = a_bar - inv(eye(m) + lambda_Newton * R) * a_bar; % 公式(20)
        sigema_power = 1/(z'*QQ*z);	% 公式(30)
        a_norm = z'*QQ_scale*z;      % 优化之后的导向矢量的范数
        PRCB(num) = real(sigema_power * a_norm / m); % 公式(28),消除模糊度
        
        num=num+1;
    end
end

% 因为f(lambda)是单调递减函数，所以利用迭代法逼近来得到lambda的结果
function [lambda_new,error]=RCB_Newton(lambda,z,Gamma_vec,m,epsilon)
    f=sum(abs(z).^2./(lambda*Gamma_vec + ones(m,1)).^2)-epsilon; %论文公式(24)
    f_derivative= -2*sum(abs(z).^2.*Gamma_vec./(lambda*Gamma_vec + ones(m,1)).^3);
    lambda_new=lambda-(f/f_derivative);
    error=sum(abs(z).^2./(lambda_new*Gamma_vec + ones(m,1)).^2)-epsilon;
end