ucaesprit_0123.m

? bpsk信号的产生可用与仿真的输入信号

% UCA ESPRIT
% Two dimensional DOA estimation through Uniform Circular Array  
% UCA ESPRIT proposed by Cherian P. Mathews and Michael D. Zoltowski,
% REFERENCE:
% C. P. Mathews and M. D. Zoltowski, 
% "Eigenstructure techniques for 2D angle estimation with uniform circular arrays",
% IEEE Trans. Signal Processing, vol. 42, no. 9, pp. 2395-2407, Sept. 1994
% The source code was completed by Zhijie Feng, Shenzhen University, January 2008.

clear all;
close all;

% ======= elevation and azimuth in degree ======== %
azimuth   = [-40 -15  13  19  -32];          % from -180 - 180
elevation = [ 60  41  28  17  2];         % from 0 - 90
% ================================================ %
snapshot = 512;
n =  size( azimuth, 2 ); % number of  signals
SNR = 20;
P = 10.^(SNR/10);
Ps = diag( ones( 1, n )*P ) ;        % power of singal
Pn = 1;         % power of noise    

N = 13;         %  number of sensors N > 2n , test for N = 5, 6 ,8 , 10 ,...
Border = 6;ceil(N/2);%5;%6;     % order of first kind Bessel function 
if Border > N/2
    Border = Border - 1;
end
M = -Border : Border;
M1 = 2*Border + 1;
                %  X is an observation data
Ntry = 200; 
AZIMUTH_mean = zeros(Ntry,length(azimuth));
ELEVATION_mean = zeros(Ntry,length(azimuth));
for mm = 1:Ntry                
    X = ucadata_huang( azimuth ,elevation ,Ps,snapshot,Pn,N);
    Rx = X*X'/snapshot;     %   Covariance of observation datas
    
    temp1 = ones( N,1 )*M;
    R = 2*pi*j*[ 0: N-1 ]'/N ; 
    temp2 = R*ones( 1, M1 );
    Wm = exp( temp1.*temp2 )/N; %   Weight matrix wm
    
    Cv = diag(j.^( -abs( M ) ) );
    V = sqrt( N )*Wm;           %   Weight matrix V
    FeH = Cv*V';                %   Eq.( 9 ) in reference
    
    alpha = 2*pi*M/M1;
    alpha = ones( M1, 1 )*alpha;
    temp3 = M'*ones( 1,M1 );
    W = exp( j*temp3.*alpha )/sqrt( M1 );   % Weight matrix
    FrH = W'*FeH;               % Real Beam Former 
    
    Ry = FrH*Rx*FrH';
    Rr = real( Ry );            % Real Beam space covariance matrix
    
%    [ U, D, Vec ] = svd( Rr );
    [Vec, Deig] = eig_sort(Rr);
    S = Vec( :, 1: n );         % Real Beam space signal subspace
    % G = Vec( :, n + 1 : end );  % Real Beam space noise subspace
    
    Co = diag( [( -1 ).^( [ Border:-1:1 ]) ones( 1, M1 - Border ) ] );
    % FoH = Co*W*FeH;
    So_hat = Co*W*S;
    % ---------------------------------- %
    Me = M1 - 2;
    S_1 = So_hat( 1: Me,: );
    S0 = So_hat( 2: Me + 1,: );
    temp = abs( 2 - Border : Border );
    D = diag( (-1).^temp );
    In = zeros( Me,Me );
    
    for i = 1:Me
        In(i, Me +1-i) = 1;
    end
    
    E2 = D*In*( conj( S_1 ) ); % 
    E = [ S_1,E2 ];  
    T = diag( [ -( Border - 1 ):( Border - 1 ) ] )/pi;
    Ls = E\( T*S0 );
    y = Ls( 1:end/2 , : )';
    [ Ui, Vi ] = eigs( y );
    EV = diag( Vi );
    AZIMUTH(mm,:) = ( angle( EV ) * 180 /pi )';
    ELEVATION(mm,:) = ( asin( abs( EV ) ) * 180/pi )';
end
AzimuthMean = mean(AZIMUTH);
ElevationMean = mean(ELEVATION);



clc
fprintf( ' @----------------------------- UCA ESPRIT -----------------------------@\n ' );
Info = sprintf('Number of elements = %d \nNumber of snapshots = %d \nNumber of signals = %d \nSNR = %d dB \nBessel = %d ',N ,snapshot,n ,SNR,Border)
fprintf( 'Original Angles are as follows: ( in degree )\n' ),
azimuth,elevation

fprintf( ' @---------------------------- Simulation Result ---------------------------@ \n' );
fprintf( 'Estimated Angles are as follows: ( in degree )\n' ),
AzimuthMean, ElevationMean
fprintf( 'Maximum Estimation Error are as follows: ( in degree )\n' ),
AzimuthMaxError = max(abs(sort(AzimuthMean) - sort(azimuth))),
ElevationMaxError = max(abs(sort(ElevationMean) - sort(elevation)))
fprintf( ' @---------------------------------- END  -----------------------------------@ ' );