基于MUSIC和NCMUSIC时延估计精度比较

 按理说NCMUSIC性能应该比MUSIC好，但是运行出来差距很大。不知道什么问题，求高手解答！

%时域MUSIC求TOA，
clc;
clear;
f0=0.3e6;                                 
f1=1.5e6;
B=f1-f0;                             
Fs=5e6;Ts=1/Fs;
T=200/Fs;
K=B/T;
nn=T/Ts;                                   %采样点数
tt=linspace(0,T,nn);
 nn=length(tt);
St=sin(2*pi*f0*tt+2*pi*K*tt.^2);          %产生线性调频信号

J=1;%试验数
SNR_db=-4:2:24;
tao=[5 11]*Ts;%两路延迟
d=length(tao);%信源数，信道数，延迟个数
tao1=tao(1)*ones(1,nn);
tao2=tao(2)*ones(1,nn);
x1=sin(2*pi*f0*tt+2*pi*K*tt.^2);
x2=sin(2*pi*f0*(tt-tao(1))+2*pi*K*((tt-tao(1)).^2));
x3=sin(2*pi*f0*(tt-tao2)+2*pi*K*((tt-tao2).^2));
A=[sin(2*pi*f0*(tt-tao1)+2*pi*K*((tt-tao1).^2))' sin(2*pi*f0*(tt-tao2)+2*pi*K*((tt-tao2).^2))'];%两路信号,LFM
MonteCarlo=100;
% M=rand(d,J);%幅度衰减矩阵
M=[1,1]';
%% MUSIC
for k = 1:length(SNR_db)  
    SNR = SNR_db(k); 
     SNR
     error1=0;
  for m=1:MonteCarlo

dPt1=sum(abs(x1).^2)/nn;
dPt2=sum(abs(x2).^2)/nn;
snr=10^(0.1*SNR);
dNoise1 = dPt1/snr;
dNoise2 = dPt2/snr;

u =sqrt(dNoise2)*randn(nn,J);
X=A*M+u;

Rx=X*X'/nn;
[EV,S]=eig(Rx);

EVA=diag(S)'; %将特征值矩阵对角线提取并转为一行
[EVA,I]=sort(EVA); %将特征值排序 从小到大
EVA=fliplr(EVA);
EV=fliplr(EV(:,I)); % 对应特征矢量排序
En = EV(:,d+1:end);

t=0:0.01*Ts:15*Ts;
aa=zeros(nn,length(t));
for i=1:length(t)
    aa(:,i)=2*pi*f0*(tt'-t(i)*ones(nn,1))+2*pi*K*(tt'-t(i)*ones(nn,1)).^2;
end
Aa=sin(aa);

den=diag(Aa'*En*En'*Aa);

P=10*log10(1./den);

[~, PeakLocation] = findpeaks(P,1:length(t),'SortStr', 'descend');
vEestimatedDelay = t(PeakLocation(1:d)).';
    vEestimatedDelay=sort(vEestimatedDelay);

 error1=error1+(vEestimatedDelay(1)-tao(1))^2+(vEestimatedDelay(2)-tao(2))^2;

  end
    RMSE_music(k)=sqrt(error1/MonteCarlo);
end
%% NCMUSIC
for k = 1:length(SNR_db)  
    SNR = SNR_db(k); 
     SNR
     error2=0;
  x2q=x2-mean(x2).*ones(1, nn);
x3q=x3-mean(x3).*ones(1, nn);
B=[x2q' x3q'];%两路信号,LFM
for m=1:MonteCarlo
     dPt=sum(abs(x1).^2)/nn;
     snr=10^(0.1*SNR);
     dNoise = dPt/snr;  % Noise power
     phi = [20 40 60 80];
     ns1 = x1 .* exp(1j * deg2rad(phi(1)));%将角从以度为单位转换为以弧度为单位
     ns2 = x2 .* exp(1j * deg2rad(phi(2)));
     ns = [ns1;ns2];
     u =sqrt(dNoise)*randn(nn,J);
     x = B * ns+u;
     y = [x;conj(x)];
     Ry = y * y' / nn;

     [EV,S]=eig(Ry);
     [EVx, Dx] = eig(Ry);                                  % 特征值分解 MUSIC-LIKE
     EVAx = diag(Dx)';
     [EVAx, I1x] = sort(EVAx);                                  % 特征值从小到大排列
     EVAx = fliplr(EVAx);                                     % 左右翻转，从大到小排列
     EVx = fliplr(EVx(:, I1x));                                     % 对应特征矢量排序
     En1 = EVx(1 : nn, d + 1 : 2 * nn);
     En2 = EVx(nn + 1 : 2 *nn, d + 1 : 2 * nn);

     t=0:0.01*Ts:10*Ts;
     aa=zeros(nn,length(t));
for i=1:length(t)
    aa=2*pi*f0*(tt'-t(i))+2*pi*K*(tt'-t(i)).^2;
    Aa=sin(aa);
    
    f1 = Aa' * En1 * En1' * Aa;
    f2 = Aa.' * En2 * En2' * conj(Aa);
    f3 = Aa' * En1 * En2' * conj(Aa);
    f4 = Aa' * En2 * En2' * Aa;
    Pmusic(i) = 1/(f1 - abs(f3));
end
P_MUSIC_dB=10*log10(abs(Pmusic));
[~, PeakLocation] = findpeaks(P_MUSIC_dB,1:length(t),'SortStr', 'descend');
vEestimatedDelay= sort(t(PeakLocation(1:d)));
 error2=error2+(vEestimatedDelay(1)-tao(1))^2+(vEestimatedDelay(2)-tao(2))^2;
end
    RMSE_ncmusic(k)=sqrt(error2/MonteCarlo);
 end

%% CRLB

for k = 1:length(SNR_db)  
    SNR = SNR_db(k); 
     SNR
aita1=1;
aita2=1;
     snr=10^(0.1*SNR);
dPt1=sum(abs(x1).^2)/nn;
dPt2=sum(abs(x2).^2)/nn;
dNoise1 = dPt1/snr;
dNoise2 = dPt2/snr;
     
sdao=((2*pi*f0+4*pi*K*tt).*cos(2*pi*f0*tt+2*pi*K*tt.^2))';
s_deltaD2dao=((2*pi*f0+4*pi*K*(tt-tao(2)).*cos(2*pi*f0*(tt-tao(2))+2*pi*K*(tt-tao(2)).^2)))';
s_deltaD1dao=((2*pi*f0+4*pi*K*(tt-tao(1)).*cos(2*pi*f0*(tt-tao(1))+2*pi*K*(tt-tao(1)).^2)))';
gamma=aita2*aita1*s_deltaD2dao'*s_deltaD1dao;
aitasigma=(aita1+aita2)^2;
CRLB(k)=((dNoise2+dNoise1*aitasigma)/2)*(1/(aita1^2*sdao'*sdao-gamma^2/(aita2^2*sdao'*sdao)));
end
% figure
% hold on
% plot(SNR_db,10*log10(RMSE_music),'-b');
% plot(SNR_db,10*log10(sqrt(CRLB)),'-r');
% 
% xlabel('SNR/dB');
% ylabel('RMSE');
% title('RMSE随信噪比变化');
% legend('MUSIC','CRLB');

figure
hold on
plot(SNR_db,RMSE_music*1e6,'-b^');
plot(SNR_db,RMSE_ncmusic*1e6,'-yo');
plot(SNR_db,sqrt(CRLB)*1e6,'-r');
xlabel('SNR/dB');
ylabel('RMSE/us');
title('RMSE随信噪比变化');
legend('MUSIC','NCMUSIC','CRLB');