%2017.6.22
%PM算法，不需要svd，快速获得噪声子空间。

clc
close all;
clear all

M=10;                                                          %x轴阵元数
d_lamda=1/2;                                                   %阵元间距
K=2;                                                           %信号数
theta=[45,80];                                                 %x轴方位角
fc=[200*10^6,500*10^6];                                        %信号频率
fs=1000*10^6;                                                  %采样频率
L=1000;                                                        %快拍数
snr=20;  %信噪比

for k=1:L
    S(:,k)=sqrt(10.^(snr/10))*randn(1)*exp(1j*2*pi*fc'/fs*(k-1));
end
for kk=1:K
    A(:,kk)=exp(-1j*2*pi*[0:M-1]'*d_lamda*sin(theta(kk)/180*pi));   %导向矢量
end

A1=A(1:K,:);
A2=A(K+1:M,:);

P=pinv(A1')*A2';                                                 %传播算子
Q=[P;-eye(M-K)];                                                 %构造噪声子空间
Gn=Q*Q';
searching_doa=-90:0.1:90;                                        %搜索范围为-90~90度
 for i=1:length(searching_doa)
   a_theta=exp(-1j*(0:M-1)'*2*pi*d_lamda*sin(pi*searching_doa(i)/180));
   Pmusic(i)=1./abs((a_theta)'*Gn*a_theta);
 end
plot(searching_doa,10*log(Pmusic),'r');