clear; close all;

da = readmatrix('./ULAData/ULA70.xlsx');
[r,~]=size(da);
angs = zeros(r/79,1);

c = physconst('LightSpeed');
    fble = 2.4e9; %240GHz
    lambda = c/fble;

dd = 0.04/lambda  ;                                                                 %阵元间距d = lambda/2
d = 0:dd:3*dd;

% for i =1:r/79
    % phase shift between two consecutive samples with a 1-μs difference is around 90°.
    % Taking the offset frequency into account as well, it might vary between 80°-100°.
    % For 1- μs slots, 2 times the value should be used.
    t= (10-1)*79;   
    mr = da(1+t:7+t,1)+ 1i * da(1+t:7+t,2); %reference samples

    %calculate phrase diffrence
    ra = angle(mr);
    pd = ra(1:6)-ra(2:7);
    mpd = mod(2*pi+pd,2*pi);
    %average or median
    pd = mean(mpd);% average
    %pd = median(mdr);% median
    dpd = rad2deg(pd);

    %phrase conpensation
    m = da(8+t:79+t,1)+ 1i * da(8+t:79+t,2);
    %amplitude
    amp=abs(m);
    rm = angle(m);

    cpd = 2*pd*(0:71)';
    cpdm = rm + cpd;
    crm = mod(cpdm,2*pi);

    dm = rad2deg(crm);
    

    %i=cos,q=sin
    cm = amp.*(cos(crm)+1i*sin(crm));

    cmx = reshape(cm,[4,18]);


Rxx = cmx*cmx'/18;   %构建协方差矩阵
% doa = musicdoa(Rxx,1);
[EV,D] = eig(Rxx);                                                          %拿到特向量EV + 特征值D  新版本matlab已经从小到大排序好了
% figure(7)
% bar3(D);title('特征值矩阵排列')
%   循环搜索特征向量正交的时候
idx = 1;
[SP,SP_inv] = deal(zeros(181,1));
scale = -90:90;
for angle_degree = scale
    a =exp(1j.*2*pi*d*sind(angle_degree)).';            %构建信号导向矢量，用共轭转至全部加负号
    En = EV(:,1:end-1);                                 %用前面的几个小特征值的特征向量
    SP(idx) = (a'*En)*(En'*a);                          %利用前面讲的正交来判断结果
    SP_inv(idx) = 1/abs((a'*En)*(En'*a));               %用倒数翻转一下 变成峰值
    idx = idx + 1;
end
SP_db = db(SP);                                         %转换为dB
SP_inv_db = db(abs(SP_inv));                            %转换为dB

figure(2)
subplot(211);plot(scale,SP_inv);
xlabel('入射角/(degree)');ylabel('空间谱/(dB)');
grid on;title('正交表示目标信号来向')
subplot(212);plot(scale,SP_inv_db);
xlabel('入射角/(degree)');ylabel('空间谱/(dB)');
grid on;title('用伪谱表示目标信号来向')

%      angs(i) = doa; 
%  end