0 1 2 3 4 5 6 7
>> for n=1:length(a)
A1(n)=[exp(-j*pi*cos(q1-(pi*a(n)/4))/lam)];
A2(n)=[exp(-j*pi*cos(q2-(pi*a(n)/4))/lam)];
A3(n)=[exp(-j*pi*cos(q3-(pi*a(n)/4))/lam)];
A4(n)=[exp(-j*pi*cos(q4-(pi*a(n)/4))/lam)];
end
>> A=[A1',A2',A3',A4']; %得出A矩阵
>> n=1:1900;
>> v1=.06;
>> v2=.02;
>> v3=.03;
>> v4=.073;
>> D=[1*cos(v1*n);1*sin(v2*n);1*sin(v3*n);1*square(v4*n)]; %四个输入信号
>> U=A*D; %总的输入信号
>> U1=(U)';
>> c=cov(U*U1); %总输入信号的协方差矩阵
>> [s,z]=eig(c); %求协方差矩阵的特征矢量及特征值
>> Vn=s(:,[5:8]); %取出与零特征值对应的特征矢量
>> ci=inv(c);
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.942991e-19.
>> bb=[1 0 0 0]';
>> Wopte=A'\bb; %求解线性方程组(求Wopte)
>> q1b=[2*pi:-2*pi/180:2*pi/180];
>> b=0:7;
>> for t=1:length(q1b)
h(t)=q1b(t);
for m=1:length(b)
Aa(m)=[exp(-j*pi*cos(h(t)-(pi*a(m)/4))/lam)];
end
A1a=Aa';
Pmusic(t)=(A1a)'*A1a*(inv ((A1a)'*Vn*(Vn)'*A1a));
Pcap(t)=inv((A1a)'*ci*(A1a));
T(t)=q1b(t);
P1=abs(Pmusic); %MUSIC算法估计
P2=abs(Pcap); %Capon算法估计
Ye(t)=Wopte'*A1a;
P=abs(Ye);
end
>> figure(1)
>> T1=T*180/pi;
>> semilogy(T1,P1);grid %MUSIC算法波达方向估计