Music algorithm with KHATRI-RAO subspace for Doa estimation

[sanjoybasak1234]

Hi,
I am trying to implement the music algorithm for rank enhanced case with the following theory. Have anyone already tried this on Matlab? If so, please let me know how it worked? Please have a look in the code.
DOA ESTIMATION OF QUASI-STATIONARY SIGNALS VIA KHATRI-RAO SUBSPACE
by Wing-Kin Ma , Tsung-Han Hsieh , and Chong-Yung Chi
I am not really getting the performance as shown with the paper. I think some configuration went wrong, but not being able to figure it. For one source signal, the result is ok, however, as the theory is about rank enhancement. For sources> sensors, I am not getting the correct result. And interestingly, for 2 sensors, the algorithm shows reduction of dynamic range, which should not happen with music algorithm. Also, the detected azimuth is not exactly correct.

best regards
Sanjoy


az = [-30; 0 ; 40 ; 90 ]; % Azimuths
M = length(az); % Number of sources

L = 25000; % Number of data snapshots recorded by receiver
m_sig = randn(M,L); % Example: normally distributed random signals

% ========= (2) RECEIVED SIGNAL ========= %

% Wavenumber vectors (in units of wavelength/2)
k = pi*[sind(az)].';

N = 4; % Number of antennas
r = (-(N-1)/2N-1)/2)'; % Linear array

% Matrix of array response vectors
A = exp(-1j*r*k);

% Additive noise
sigma2 = 0.00; % Noise variance
n = sqrt(sigma2)*(randn(N,L) + 1j*randn(N,L))/sqrt(2);

% Received signal
signal_rx = A*m_sig + n;


T=length(signal_rx);
L=200; %length of each frame

M=T/L; % length of each segment


for m=1:M

signal_rx_new,:,m)=signal_rx,1:L);
Rxx_nw,:,m)=signal_rx_new,:,m)*(signal_rx_new,:,m))'/L;
RXX_newvec,m)=vec(Rxx_nw,:,m)); % vectorize
signal_rx,1:L)=[];

end

% P_1M construction
one_M=ones(1,M)';
I_M=eye(M);

P_1M=I_M-(one_M*one_M'/M);


Y_bar=RXX_newvec*P_1M; % noise covariance elimination

[E,D] = svd(Y_bar); % performing SVD

nb_sources=length(az);
Noise_subspace = E,1+nb_sources:end); % Noise subspace extraction


%% MUSIC search directions
AzSearch = (-90:1:90); % Azimuth values to search

% Corresponding points on array manifold to search

for i=1:length(AzSearch)

aznw=AzSearch(i);
kSearch = pi*[sind(aznw)].';
ASearch = exp(-1j*r*kSearch);

a_thetaprod=kron(conj(ASearch),ASearch); % taking Kronecker product

Z_testvalue=Noise_subspace'*a_thetaprod;

Z,i) = sum(abs(Z_testvalue).^2,1);

end


% Plot
figure();
plot(AzSearch,-10*log10(Z));
title('Simple 1D MUSIC Example');
xlabel('Azimuth (degrees)');
ylabel('MUSIC spectrum (dB)');
grid on; axis tight;