www.pudn.com > gene-selection.zip > GENE_AH.m, change:2012-08-10,size:2611b


%10th August 2012 
%===================================================================== 
% Programmer:  
% ArulMurugan Ambikapathi, E-mail: a.arulmurugan@gmail.com 
% ------------------------------------------------------- 
% Reference:  
% A.Ambikapathi, T.-H. Chan, C.-Y. Chi, and  K. Keizer, ``Hyperspectral 
% data geometry based estimation of number of endmembers using p-norm 
% based pure pixel identification algorithm,'' to appear 
% in IEEE Trans. Geoscience and Remote Sensing. 
%====================================================================== 
% An implementation of GENE-AH Algorithm 
% [MODELORDER] = GENE_AH(X,P_FA,N_max,CNN) 
%====================================================================== 
%  Input 
%  X is M-by-L data matrix where M is the spectral bands (or observations) and L is the number of pixels (data length).    
%  P_FA is the pre-assigned false alarm probability, say 10^{-6}. 
%  N_max is maximum bound on the number of endmembers, where N <= N_max <= M, say N_max=25. 
%  CNN is an estimate of the noise covariance matrix. 
%---------------------------------------------------------------------- 
%  Output 
%  MODELORDER is an estimate of the number of endmembers 
%---------------------------------------------------------------------- 
%%%%%P.S: the estimated number of endmembers using GENE-AH-MOD corresponds to 
%one less than that of GENE-AH (for proof, please see Lemma 1 and the discussions following Lemma 1 of the paper). 
%======================================================================== 
 
function [MODELORDER] = GENE_AH(X,P_FA,N_max,CNN) 
Rw=CNN; 
[M,L] = size(X);  
d = mean(X,2); 
U = X-d*ones(1,L); 
R = U*U'; 
Rx_true = (R./L) - Rw; 
[eV D] = eig(Rx_true); 
C = eV(:,end:-1:M-(N_max-2)); 
 
Xd_t = C'*(X-d*ones(1,L));  % Dimension reduced data vectors (Xd is (N-1)-by-L) 
Rw_small= C'*Rw*C; 
 
for N=3:N_max, 
[PP_Index]=TRIP(X,N,2,N_max,CNN); 
 
  
 ALPHA_est = Xd_t(:,PP_Index); 
 A_k_1 = ALPHA_est(:,[1:N-1]); 
 a_k= ALPHA_est(:,N); 
 clear yessame; 
if (N>3) 
    % Do the following optimization problem 
    cvx_quiet(true) 
    cvx_begin    
    variable theta_opt(N-1,1); 
    minimize norm((a_k-A_k_1*theta_opt)); 
    subject to              
         ones(1,N-1)*theta_opt==1; 
    cvx_end 
 
    b_vector= a_k-A_k_1*theta_opt; 
    zeta= theta_opt'*theta_opt ;   
    var_b = (1+zeta)*Rw_small; 
    r=(b_vector'*inv(var_b)*b_vector); 
    format long e 
    zero_to_r = gammainc(r/2,(N_max-1)/2,'lower'); 
r_to_inf=1- zero_to_r; 
 
if (r_to_inf>=P_FA) 
          MODELORDER=N-1; break; else MODELORDER=0; 
end 
end 
end