music_doa_3.m

阵列信号处理

close all
clear all
clc
%模拟信号
N=14;     %阵元数
K=200;    %快拍数
M=3;      %信源数
theta=[10,30,80]; %第一个是期望信号，后面的都是干扰
Xt=zeros(N,K);
Snr=0;     %信噪比
Sir1=-10;  %信干比1
Sir2=-20;  %信干比2
Pd=1;      %给出信号功率
Pn=Pd/(10^(Snr/10));    %计算噪声功率
Pi1=Pd/(10^(Sir1/10));  %计算干扰1功率
Pi2=Pd/(10^(Sir2/10));  %计算干扰2功率
theta=theta*pi/180;
ai1=exp(-j*[0:N-1]'*sin(theta(2))*pi);       %阵列导向矢量（干扰信号1）
ai2=exp(-j*[0:N-1]'*sin(theta(3))*pi);       %阵列导向矢量（干扰信号2）
ad=exp(-j*[0:N-1]'*sin(theta(1))*pi);        %阵列导向矢量（期望信号）
Sd=sqrt(Pd/2)*(randn(1,K)+j*randn(1,K));     %构造期望信号包络
Si1=sqrt(Pi1/2)*(randn(1,K)+j*randn(1,K));   %构造干扰信号1包络
Si2=sqrt(Pi2/2)*(randn(1,K)+j*randn(1,K));   %构造干扰信号2包络
Xt=ad*Sd+ai1*Si1+ai2*Si2;                    %构造输入信号（期望信号+干扰）
Xt=Xt+sqrt(Pn/2)*(randn(N,K)+j*randn(N,K));  %构造输入信号（期望信号+干扰+噪声）
%以上是模拟信号产生
%Music DOA
Rx=Xt*Xt'/K;
[a,b,c]=svd(Rx);
Un=a(:,M+1:N);
Pn=Un*Un';
thetas=-90:90;
for m=1:length(thetas)
    tm=thetas(m)*pi/180;
    am=exp(-j*pi*sin(tm)*[0:N-1]');
    S_Music(m)=1/norm(Pn*am);
end
%Capon Beamformming
Rx=Xt*Xt'/K;
Wopt=inv(Rx)*ad;
thetas=-90:90;
for m=1:length(thetas)
    tm=thetas(m)*pi/180;
    am=exp(-j*pi*[0:N-1]'*sin(tm));
    S_Capon(m)=10*log(abs(Wopt'*am));
end
plot(thetas,S_Capon,'r--','linewidth',2)
grid
hold on
plot(thetas,S_Music,'b-','linewidth',2)
legend('Capon','Music',2)