tic                                 %程序計時開始
clear;
clc;

n=0:1/255/1000:1/1000;              %時間長度1ms,采樣頻率256KHz

w1=2000*pi
s1=sin(w1*n);
subplot(3,1,1)
plot(s1)		                     %輸入正弦信號1,頻率f1=1KHz

w2=4000*pi
s2=sin(w2*n)
subplot(3,1,2)
plot(s2)		                     %輸入正弦信號2,頻率f2=2KHz

S=[s1;s2];                           %輸入信號矩陣

cita1=pi/6;                          %信號1方位角30度
cita2=3*pi/8;                        %信號2方位角135度

M=[0:5]			                     %陣元數M=6
ratio=1/2                            %波長間距比

a1=exp(j*2*pi*ratio*M*cos(cita1))
a2=exp(j*2*pi*ratio*M*cos(cita2))
 A=[a1.' a2.'];		                 %方位角矩陣

N=randn(length(M),length(n))    	 %白噪聲矩陣

X=A*S+1.5*N			                 %陣元輸出信號

Rx=X*X'/length(s1)		             %求自相關矩陣
    
[V,D]=eig(Rx)			             %自相關矩陣的特征分解,,得到特征值矩陣和與特征值一一對應的特征向量矩陣

for i=1:length(M)                    %對特征值排序
    P(i)=D(i,i)
end

[Q,index]=sort(P)                    %將與特征值對應的特征向量也按照特征值的順序排列

for i=1:length(M)-2                  %取排了序的特征向量的前4個來構造噪聲子空間En
    En(:,i)=V(:,index(i))
end

a=[0:180]                           
m=M'
a_cita=exp(j*2*pi*m*ratio*cos(pi*a/180))

temp=En'*a_cita

for b=1:180;
    fanshu=norm(temp(:,b));
    out(b)=1/fanshu^2;		        %進行譜峰搜索
end;

subplot(3,1,3)
plot(out)

toc                                 %%程序計時結束