一、均勻圓陣(UCA, Uniform Circular Array)的MUSIC算法
假設一個半徑為R的M元均勻圓陣的所有陣元均位於坐標系X-Y平面內，第k-1個陣元坐標為
，第i個窄帶信號波長為
，來波方向為
，如圖1，則第k-1個陣元到圓心(即原點)的波程差
為：
均勻圓陣
存在P個入射信號均勻圓陣的接收模型可表示為：
其他步驟與基於ULA的MUSIC算法一致。
令任意兩陣元間的波程差
為：
當
時，即產生相位模糊。將均勻圓陣各陣元投影到入射方向，得到一個隨入射方向變動的非均勻線陣。需要保證在任意入射方向上投影出的非均勻線陣，其最小間隔總是小於信號波長，模糊譜峰對測向結果影響較小，即：
在方向角相同時，水平入射(
)信號的波程差最長，且投影出的非均勻線陣隨方向角不同周期變化，因此只需要討論水平入射信號對應投影線陣的不同情況。
在MUSIC算法中，陣元的最小間隔越大模糊譜峰峰值就越大。但在均勻圓陣中，陣元間隔隨着入射波方向變化，因此算法性能受到最小間隔最大值的影響。根據來波方向不同，入射方向上的第k個陣元投影間隔分別為：
當M為奇數時，對着陣元入射，投影點重合為${\rm{(M + 1)/2}}$個；當M為偶數時，對着相鄰陣元連線中點入射，投影點重合為M/2個，此時投影線陣的非零最小間隔的值最大，且取得該最大值k=1時。進一步可求得半徑的選取關系：
選取半徑時，按上式等比例縮放，即能使對應的奇數陣與偶數陣有近似的抗相位模糊的性能。
每一路接收的結構圖：
二、MUSIC
1、
clear all;
%產生三信源，角度分別為-40°、30°、45°，采用8PSK調制，滾降系數為0.5的平方根升余弦濾波
Nsym=500;%符號個數
Fsym=1;%符號速率
M=3;%一個符號對應的比特數
Fbit=M*Fsym;%比特速率
Nsour=3;%信源數
Angle=[5,15,35];%信源的來波方向
Fc=10;%載波頻率
Fs=100;%抽樣頻率
R=0.5;%滾降因子
Del=5;%群延遲因子
% Nsamp=50;%采樣點數或者快拍數
S1=randint(Nsym,1,2^M);
S2=randint(Nsym,1,2^M);
S3=randint(Nsym,1,2^M);
PM1=pmmod(S1,Fc,Fs,pi/8,pi/4);
PM2=pmmod(S2,Fc,Fs,pi/8,pi/4);
PM3=pmmod(S3,Fc,Fs,pi/8,pi/4);
Rcos11=rcosflt(PM1,Fsym,Fs,'fir/sqrt/Fs',R,Del);
Rcos21=rcosflt(PM2,Fsym,Fs,'fir/sqrt/Fs',R,Del);
Rcos31=rcosflt(PM3,Fsym,Fs,'fir/sqrt/Fs',R,Del);
Rcos1=0.99*Rcos11+Rcos21+1.02*Rcos31;%構造相干信源--信源1、信源2與信源3
Rcos2=Rcos11+Rcos21+Rcos31;%構造相干信源--信源1、信源2與信源3
Rcos3=Rcos11+1.03*Rcos21+1.05*Rcos31;%構造相干信源--信源1、信源2與信源3
save xyc3 Rcos1 Rcos2 Rcos3
%產生三信源，角度分別為-40°、30°、45°，采用8PSK調制，滾降系數為0.5的平方根升余弦濾波
Nsamp=1024;%采樣點數或者快拍數
i=sqrt(-1);
j=i;
Ntx=8;%陣列數
SNR=[2,2,2];%三信源的信噪比
% sn=10; %----單信號源
Lamda=2;%信號波長
D=Lamda/2;%線性陣列的距離
p=3;%子陣個數
L=Ntx-p+1;%子陣陣元數
nr=randn(Ntx,Nsamp);
ni=randn(Ntx,Nsamp);
n=nr+j*ni;%產生背景噪聲
load xyc3;
t=1:Nsamp;
% s1=[Rcos1(t).'];%接收信號的采樣點數%----單信號源
s1=[Rcos1(t).';Rcos2(t).';Rcos3(t).'];%矩陣維數=信源數*抽樣點數
ps=diag((s1*s1')/Nsamp);%無噪聲信號功率--%矩陣維數=信源數*1
delta1=(1./(2*10.^(SNR/10)))*ps;%矩陣維數=1*1
% delta1=ps./(2*10.^(sn/10)); %----單信號源
delta2=diag(delta1);%矩陣維數=1*1
delta=sqrt(delta2);%噪聲幅度值--%矩陣維數=1*1
Rev_s1=(1./delta')*s1;%SNR條件下的信號幅度--%矩陣維數=信源數*抽樣點數
%計算各信源SNR比條件下，陣列接收到的信號幅度%
Pn=zeros(Nsamp,1);
pn=zeros(Ntx,Nsamp);
Pn=diag(n'*n);
for h=1:Nsamp
pn(:,h)=n(:,h)./sqrt(Pn(h,:));
end
Rev_n=pn;
%計算各陣列接收到的背景噪聲下的信號幅度%
tmp=-j*2*pi*D*sin(Angle*pi/180)/Lamda;%---%矩陣維數=1*信源數
% tmp=-j*2*pi*D*sin(1*pi/180)/Lamda; %----單信號源
tmp1=[0:Ntx-1]';%矩陣維數=陣元數*1
tmp4=[0:L-1]';%子矩陣維數=子矩陣陣元數*1
a1=tmp1*tmp;%矩陣維數=陣元數*信源數
A=exp(a1);%方向矩陣--%矩陣維數=陣元數*信源數
X=A*Rev_s1+Rev_n;%陣列接收到的信號幅度--%矩陣維數=陣元數*抽樣點數
Rxx=(X*X')/Nsamp;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%空間平滑算法%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sub_FRxx=zeros(L,L);
sub_BRxx=zeros(L,L);
for i=1:p
sub_FR=zeros(L,Nsamp);
sub_BR=zeros(L,Nsamp);
sub_FR=X(i:1:i+L-1,:);
last=Ntx+1-i;
first=last-L+1;
sub_BR=conj(X(last:-1:first,:));
sub_FRxx=sub_FRxx+((sub_FR*sub_FR')./Nsamp);
sub_BRxx=sub_BRxx+((sub_BR*sub_BR')./Nsamp);
end
sub_FRxx=sub_FRxx./p;
sub_BRxx=sub_BRxx./p;
sub_Rxx=(sub_FRxx+sub_BRxx)./2;
[VFB,HFB]=eig(sub_Rxx);
[HFB,IFB]=sort(diag(HFB),1);
VFB=VFB(:,IFB);
VnFB=VFB(:,1:L-Nsour);
VsFB=VFB(:,L-Nsour+1:L);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%空間平滑算法%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ScanAng=[-90:1:90];
for i=1:length(ScanAng)
tmp2=-j*2*pi*D*sin(ScanAng(i)*pi/180)/Lamda;
tmp3=tmp2*tmp1;
tmp5=tmp2*tmp4;
A_Sita=exp(tmp3);
Sub_Sita=exp(tmp5);
Sub_FBsita(i)=(Sub_Sita'*Sub_Sita)/(Sub_Sita'*VnFB*VnFB'*Sub_Sita);
end
figure(1);
semilogy(ScanAng,real(Sub_FBsita),'bo-');
axis([-60 60 0.1 1e7]);
xlabel('M_Angle(deg)');
ylabel('M_Spectrum');
grid on
2、
clear all;
%產生三信源，角度分別為-40°、30°、45°，采用8PSK調制，滾降系數為0.5的平方根升余弦濾波
Nsym=500;%符號個數
Fsym=1;%符號速率
M=3;%一個符號對應的比特數
Fbit=M*Fsym;%比特速率
Nsour=3;%信源數
Angle=[10,40,80];%信源的來波方向
Fc=10;%載波頻率
Fs=100;%抽樣頻率
R=0.5;%滾降因子
Del=5;%群延遲因子
% Nsamp=50;%采樣點數或者快拍數
S1=randint(Nsym,1,2^M);
S2=randint(Nsym,1,2^M);
S3=randint(Nsym,1,2^M);
PM1=pmmod(S1,Fc,Fs,pi/8,pi/4);
PM2=pmmod(S2,Fc,Fs,pi/8,pi/4);
PM3=pmmod(S3,Fc,Fs,pi/8,pi/4);
Rcos11=rcosflt(PM1,Fsym,Fs,'fir/sqrt/Fs',R,Del);
Rcos21=rcosflt(PM2,Fsym,Fs,'fir/sqrt/Fs',R,Del);
Rcos31=rcosflt(PM3,Fsym,Fs,'fir/sqrt/Fs',R,Del);
Rcos1=Rcos11;%構造相干信源--信源1、信源2與信源3
Rcos2=Rcos21;%構造相干信源--信源1、信源2與信源3
Rcos3=Rcos31;%構造相干信源--信源1、信源2與信源3
save xyc3 Rcos1 Rcos2 Rcos3
%產生三信源，角度分別為-40°、30°、45°，采用8PSK調制，滾降系數為0.5的平方根升余弦濾波
Nsamp=1024;%采樣點數或者快拍數
i=sqrt(-1);
j=i;
Ntx=8;%陣列數
SNR=[10,10,10];%三信源的信噪比
% sn=10; %----單信號源
Lamda=2;%信號波長
D=Lamda/2;%線性陣列的距離
p=3;%子陣個數
L=Ntx-p+1;%子陣陣元數
nr=randn(Ntx,Nsamp);
ni=randn(Ntx,Nsamp);
n=nr+j*ni;%產生背景噪聲
load xyc3;
t=1:Nsamp;
% s1=[Rcos1(t).'];%接收信號的采樣點數%----單信號源
s1=[Rcos1(t).';Rcos2(t).';Rcos3(t).'];%矩陣維數=信源數*抽樣點數
ps=diag((s1*s1')/Nsamp);%無噪聲信號功率--%矩陣維數=信源數*1
delta1=(1./(2*10.^(SNR/10)))*ps;%矩陣維數=1*1
% delta1=ps./(2*10.^(sn/10)); %----單信號源
delta2=diag(delta1);%矩陣維數=1*1
delta=sqrt(delta2);%噪聲幅度值--%矩陣維數=1*1
Rev_s1=(1./delta')*s1;%SNR條件下的信號幅度--%矩陣維數=信源數*抽樣點數
%計算各信源SNR比條件下，陣列接收到的信號幅度%
Pn=zeros(Nsamp,1);
pn=zeros(Ntx,Nsamp);
Pn=diag(n'*n);
for h=1:Nsamp
pn(:,h)=n(:,h)./sqrt(Pn(h,:));
end
Rev_n=pn;
%計算各陣列接收到的背景噪聲下的信號幅度%
tmp=-j*2*pi*D*sin(Angle*pi/180)/Lamda;%---%矩陣維數=1*信源數
% tmp=-j*2*pi*D*sin(1*pi/180)/Lamda; %----單信號源
tmp1=[0:Ntx-1]';%矩陣維數=陣元數*1
tmp4=[0:L-1]';%子矩陣維數=子矩陣陣元數*1
a1=tmp1*tmp;%矩陣維數=陣元數*信源數
A=exp(a1);%方向矩陣--%矩陣維數=陣元數*信源數
X=A*Rev_s1+Rev_n;%陣列接收到的信號幅度--%矩陣維數=陣元數*抽樣點數
Rxx=(X*X')/Nsamp;
[V,H]=eig(Rxx);%MUSIC算法---MUltiSIgnal Classification
[H,I]=sort(diag(H),1);%特征值按照升序排列
V=V(:,I);%特征值對應的特征向量也按照相應特征值的升序排列
Vn=V(:,1:Ntx-Nsour);%噪聲子空間---協方差的特征向量--最小特征值對應的特征向量
Vs=V(:,Ntx-Nsour+1:Ntx);%信號子空間---協方差的特征向量--最大特征值對應的特征向量
ScanAng=[-90:1:90];
for i=1:length(ScanAng)
tmp2=-j*2*pi*D*sin(ScanAng(i)*pi/180)/Lamda;
tmp3=tmp2*tmp1;
tmp5=tmp2*tmp4;
A_Sita=exp(tmp3);
MUSIC_Sita(i)=(A_Sita'*A_Sita)/(A_Sita'*Vn*Vn'*A_Sita);
end
figure(1);
semilogy(ScanAng,real(MUSIC_Sita),'g*-');
axis([-90 90 0.1 1e7]);
xlabel('M_Angle(deg)');
ylabel('M_Spectrum');
grid on
3、
clear all;
%產生三信源，角度分別為-40°、30°、45°，采用8PSK調制，滾降系數為0.5的平方根升余弦濾波
Nsym=500;%符號個數
Fsym=1;%符號速率
M=3;%一個符號對應的比特數
Fbit=M*Fsym;%比特速率
Nsour=3;%信源數
Angle=[5,8,35];%信源的來波方向
Fc=10;%載波頻率
Fs=100;%抽樣頻率
R=0.5;%滾降因子
Del=5;%群延遲因子
% Nsamp=50;%采樣點數或者快拍數
S1=randint(Nsym,1,2^M);
S2=randint(Nsym,1,2^M);
S3=randint(Nsym,1,2^M);
PM1=pmmod(S1,Fc,Fs,pi/8,pi/4);
PM2=pmmod(S2,Fc,Fs,pi/8,pi/4);
PM3=pmmod(S3,Fc,Fs,pi/8,pi/4);
Rcos11=rcosflt(PM1,Fsym,Fs,'fir/sqrt/Fs',R,Del);
Rcos21=rcosflt(PM2,Fsym,Fs,'fir/sqrt/Fs',R,Del);
Rcos31=rcosflt(PM3,Fsym,Fs,'fir/sqrt/Fs',R,Del);
Rcos1=0.99*Rcos11+Rcos21+1.02*Rcos31;%構造相干信源--信源1、信源2與信源3
Rcos2=Rcos11+Rcos21+Rcos31;%構造相干信源--信源1、信源2與信源3
Rcos3=Rcos11+1.03*Rcos21+1.05*Rcos31;%構造相干信源--信源1、信源2與信源3
save xyc3 Rcos1 Rcos2 Rcos3
%產生三信源，角度分別為-40°、30°、45°，采用8PSK調制，滾降系數為0.5的平方根升余弦濾波
Nsamp=512;%采樣點數或者快拍數
i=sqrt(-1);
j=i;
Ntx=8;%陣列數
SNR=[2,2,2];%三信源的信噪比
% sn=10; %----單信號源
Lamda=2;%信號波長
D=Lamda/2;%線性陣列的距離
p=3;%子陣個數
L=Ntx-p+1;%子陣陣元數
nr=randn(Ntx,Nsamp);
ni=randn(Ntx,Nsamp);
n=nr+j*ni;%產生背景噪聲
load xyc3;
t=1:Nsamp;
% s1=[Rcos1(t).'];%接收信號的采樣點數%----單信號源
s1=[Rcos1(t).';Rcos2(t).';Rcos3(t).'];%矩陣維數=信源數*抽樣點數
ps=diag((s1*s1')/Nsamp);%無噪聲信號功率--%矩陣維數=信源數*1
delta1=(1./(2*10.^(SNR/10)))*ps;%矩陣維數=1*1
% delta1=ps./(2*10.^(sn/10)); %----單信號源
delta2=diag(delta1);%矩陣維數=1*1
delta=sqrt(delta2);%噪聲幅度值--%矩陣維數=1*1
Rev_s1=(1./delta')*s1;%SNR條件下的信號幅度--%矩陣維數=信源數*抽樣點數
%計算各信源SNR比條件下，陣列接收到的信號幅度%
Pn=zeros(Nsamp,1);
pn=zeros(Ntx,Nsamp);
Pn=diag(n'*n);
for h=1:Nsamp
pn(:,h)=n(:,h)./sqrt(Pn(h,:));
end
Rev_n=pn;
%計算各陣列接收到的背景噪聲下的信號幅度%
tmp=-j*2*pi*D*sin(Angle*pi/180)/Lamda;%---%矩陣維數=1*信源數
% tmp=-j*2*pi*D*sin(1*pi/180)/Lamda; %----單信號源
tmp1=[0:Ntx-1]';%矩陣維數=陣元數*1
tmp4=[0:L-1]';%子矩陣維數=子矩陣陣元數*1
a1=tmp1*tmp;%矩陣維數=陣元數*信源數
A=exp(a1);%方向矩陣--%矩陣維數=陣元數*信源數
X=A*Rev_s1+Rev_n;%陣列接收到的信號幅度--%矩陣維數=陣元數*抽樣點數
Rxx=(X*X')/Nsamp;
Rxx_fb=zeros(L,L);
Rxx_f=zeros(L,L);
Rxx_b=zeros(L,L);
J=fliplr(eye(L));
for m=1:p
for k=1:p
Rxx_f=Rxx_f+Rxx(m:1:m+L-1,k:1:k+L-1)*Rxx(k:1:k+L-1,m:1:m+L-1);
Rxx_b=Rxx_b+J*conj(Rxx(m:1:m+L-1,k:1:k+L-1))*conj(Rxx(k:1:k+L-1,m:1:m+L-1))*J;
end
end
Rxx_f=Rxx_f./p;
Rxx_b=Rxx_b./p;
Rxx_fb=(Rxx_f+Rxx_b)./p;
[V_fb,H_fb]=eig(Rxx_fb);%特征分解---MUltiSIgnal Classification
[H_fb,I_fb]=sort(diag(H_fb),1);%特征值按照升序排列
V_fb=V_fb(:,I_fb);%特征值對應的特征向量也按照相應特征值的升序排列
Vn_fb=V_fb(:,1:L-Nsour);
Vs_fb=V_fb(:,L-Nsour+1:L);
ScanAng=[-90:1:90];
for i=1:length(ScanAng)
tmp2=-j*2*pi*D*sin(ScanAng(i)*pi/180)/Lamda;
tmp3=tmp2*tmp1;
tmp5=tmp2*tmp4;
A_Sita=exp(tmp3);
Sub_Sita=exp(tmp5);
fb_sita(i)=(Sub_Sita'*Sub_Sita)/(Sub_Sita'*Vn_fb*Vn_fb'*Sub_Sita);
end
figure(1);
semilogy(ScanAng,real(fb_sita),'r-');
axis([-60 60 0.1 1e7]);
xlabel('M_Angle(deg)');
ylabel('M_Spectrum');
grid on
4、
%=========================================================================
% UCA_multi_in_2D
%
%=========================================================================
clc;
clear all;
close all;
%------------------------常數表-------------------------------
c = 3e8;
namda = c/18e9;
est_num = 1;
iteration = 100;
sr_array = [-50,-47.5,-45,-42.5,-40,-35,-27.5];
phi = 60;
%% ----------------入射信號模型-------------------------------
N_x = 2^5; %快拍點數
F0 = 18e9; %中心頻率
B = 20e6; %帶寬
Fs = 2*B; %采樣頻率
Ts = 1/Fs; %采樣時間
T = (N_x-1)*Ts; %快拍持續時間
u = B/T; %頻率變化率
t = -T/2:Ts:T/2; %時間軸點
l = c/18e9;
st = exp(1j*2*pi*(F0*t+.5*u*t.^2));
dir7 =(46:.25:56)*pi/180; %(-50:.25:-40)*pi/180;
dir8 =(62.5:.25:72.5)*pi/180;
dir9 =(35:.25:45)*pi/180;%(-51:.25:-41)*pi/180;
dir10=(49:.25:59)*pi/180;
ang=(50:.25:70)*pi/180;
e_dir7 = zeros(1,length(sr_array));
e_dir8 = zeros(1,length(sr_array));
e_dir9 = zeros(1,length(sr_array));
e_dir10= zeros(1,length(sr_array));
for ss = 1:length(sr_array)
snr = sr_array(ss);
%--------------------7陣元---------------------------------------
sensor_num = 7;
d_angle = (0:sensor_num-1)'*2*pi/sensor_num;
R = namda/(1-cos(d_angle(2)));
x = R*cos(d_angle);y = R*sin(d_angle);
theta = d_angle(2)*180/pi;
for it = 1:iteration
n = (randn(sensor_num,N_x)+1j*randn(sensor_num,N_x))/sqrt(2)*sqrt(10^(-45/10));
tao = x*(sin(phi*pi/180).*cos(theta*pi/180))+y*(sin(phi*pi/180).*sin(theta*pi/180));
A = exp(-1j*2*pi*tao./l);
xt = A*(sqrt(10^(snr/10))*st)+n;
% -------------------2D-MUSIC算法-----------------------
Rxx = xt*xt'/N_x;
[U,S] = svd(Rxx);
[~,index] = sort(diag(S));
Un = U(:,index(1:sensor_num-est_num));
Gn = Un*Un';
Pmusic = zeros(length(dir7),length(ang));
for i=1:length(dir7)
for k=1:length(ang)
a_tao = sin(ang(k))*(x*cos(dir7(i))+y*sin(dir7(i)));
a_theta = exp(-1j*2*pi*a_tao/l);
Pmusic(i,k)=1./abs((a_theta)'*Gn*a_theta);
end
end
[a,b]=find(Pmusic==max(max(Pmusic)));
aa = dir7(a)*180/pi;
e_dir7(ss) = e_dir7(ss)+(aa-theta)^2;
end
%--------------------8陣元-------------------------------------
sensor_num = 8;
d_angle = (0:sensor_num-1)'*2*pi/sensor_num;
R = namda/(2*sin(d_angle(2))*sin(0.5*d_angle(2)));
x = R*cos(d_angle);y = R*sin(d_angle);
theta = 1.5*d_angle(2)*180/pi;
for it = 1:iteration
tao = x*(sin(phi*pi/180).*cos(theta*pi/180))+y*(sin(phi*pi/180).*sin(theta*pi/180));
A = exp(-1j*2*pi*tao./l);
n = (randn(sensor_num,N_x)+1j*randn(sensor_num,N_x))/sqrt(2)*sqrt(10^(-45/10));
xt = A*(sqrt(10^(snr/10))*st)+n;
% -------------------2D-MUSIC算法-----------------------
Rxx = xt*xt'/N_x;
[U,S] = svd(Rxx);
[~,index] = sort(diag(S));
Un = U(:,index(1:sensor_num-est_num));
Gn = Un*Un';
Pmusic = zeros(length(dir8),length(ang));
for i=1:length(dir8)
for k=1:length(ang)
a_tao = sin(ang(k))*(x*cos(dir8(i))+y*sin(dir8(i)));
a_theta = exp(-1j*2*pi*a_tao/l);
Pmusic(i,k)=1./abs((a_theta)'*Gn*a_theta);
end
end
[a,b]=find(Pmusic==max(max(Pmusic)));
aa = dir8(a)*180/pi;
e_dir8(ss) = e_dir8(ss)+(aa-theta)^2;
end
%---------------9陣元---------------------------------------------------
sensor_num = 9;
d_angle = (0:sensor_num-1)'*2*pi/sensor_num;
R = namda/(1-cos(d_angle(2)));
x = R*cos(d_angle);y = R*sin(d_angle);
theta = d_angle(2)*180/pi;
for it = 1:iteration
n = (randn(sensor_num,N_x)+1j*randn(sensor_num,N_x))/sqrt(2)*sqrt(10^(-45/10));
tao = x*(sin(phi*pi/180).*cos(theta*pi/180))+y*(sin(phi*pi/180).*sin(theta*pi/180));
A = exp(-1j*2*pi*tao./l);
xt = A*(sqrt(10^(snr/10))*st)+n;
% -------------------2D-MUSIC算法-----------------------
Rxx = xt*xt'/N_x;
[U,S] = svd(Rxx);
[~,index] = sort(diag(S));
Un = U(:,index(1:sensor_num-est_num));
Gn = Un*Un';
Pmusic = zeros(length(dir9),length(ang));
for i=1:length(dir9)
for k=1:length(ang)
a_tao = sin(ang(k))*(x*cos(dir9(i))+y*sin(dir9(i)));
a_theta = exp(-1j*2*pi*a_tao/l);
Pmusic(i,k)=1./abs((a_theta)'*Gn*a_theta);
end
end
[a,b]=find(Pmusic==max(max(Pmusic)));
aa = dir9(a)*180/pi;
e_dir9(ss) = e_dir9(ss)+(aa-theta)^2;
end
%---------------10陣元---------------------------------------------------
sensor_num = 10;
d_angle = (0:sensor_num-1)'*2*pi/sensor_num;
R = namda/(2*sin(d_angle(2))*sin(0.5*d_angle(2)));
x = R*cos(d_angle);y = R*sin(d_angle);
theta = 1.5*d_angle(2)*180/pi;
for it = 1:iteration
tao = x*(sin(phi*pi/180).*cos(theta*pi/180))+y*(sin(phi*pi/180).*sin(theta*pi/180));
A = exp(-1j*2*pi*tao./l);
n = (randn(sensor_num,N_x)+1j*randn(sensor_num,N_x))/sqrt(2)*sqrt(10^(-45/10));
xt = A*(sqrt(10^(snr/10))*st)+n;
% -------------------2D-MUSIC算法-----------------------
Rxx = xt*xt'/N_x;
[U,S] = svd(Rxx);
[~,index] = sort(diag(S));
Un = U(:,index(1:sensor_num-est_num));
Gn = Un*Un';
Pmusic = zeros(length(dir10),length(ang));
for i=1:length(dir10)
for k=1:length(ang)
a_tao = sin(ang(k))*(x*cos(dir10(i))+y*sin(dir10(i)));
a_theta = exp(-1j*2*pi*a_tao/l);
Pmusic(i,k)=1./abs((a_theta)'*Gn*a_theta);
end
end
[a,b]=find(Pmusic==max(max(Pmusic)));
aa = dir10(a)*180/pi;
e_dir10(ss) = e_dir10(ss)+(aa-theta)^2;
end
end
figure;
% subplot(121);
% plot(sr_array+45,e_dir10(1,:)/iteration,'-^k',sr_array+45,e_dir9(1,:)/iteration,'-*k');
% legend('8元UCA','9元UCA');
% grid on;%axis([-10,20,-.01,.45]);
% colormap gray;
% xlabel('信噪比/dB');ylabel('均方誤差/°');
% title('半徑相同精度試驗');
%
%
% subplot(122);
plot(sr_array+45,e_dir7/iteration,'-^k',sr_array+45,e_dir8/iteration,'-*k');hold on;
plot(sr_array+45,e_dir9/iteration,'-sk',sr_array+45,e_dir10/iteration,'-dk');
legend('7陣元','8陣元','9陣元','10陣元');
grid on;axis([-5.5,18,-.02,.65]);
colormap gray;
xlabel('信噪比/dB');ylabel('均方誤差/°');
title('不同陣元最大半徑測向試驗');
5、
%=========================================================================
% UCA_multi_in_2D
%
%=========================================================================
clc;
clear all;
close all;
%------------------------常數表-------------------------------
c = 3e8;
namda = c/18e9;
est_num = 1;
iteration = 100;
sr_array = [-50,-47.5,-45,-42.5,-40,-35,-27.5];
phi = 60;
%% ----------------入射信號模型-------------------------------
N_x = 2^5; %快拍點數
F0 = 18e9; %中心頻率
B = 20e6; %帶寬
Fs = 2*B; %采樣頻率
Ts = 1/Fs; %采樣時間
T = (N_x-1)*Ts; %快拍持續時間
u = B/T; %頻率變化率
t = -T/2:Ts:T/2; %時間軸點
l = c/18e9;
st = exp(1j*2*pi*(F0*t+.5*u*t.^2));
R = 0.01;
dir9(1,:)=(-50:.25:-40)*pi/180;
dir9(2,:)=(75:.25:85)*pi/180;
dir10(1,:)=(-51:.25:-41)*pi/180;
dir10(2,:)=(107:.25:117)*pi/180;
ang=(50:.25:70)*pi/180;
e_dir9= zeros(2,length(sr_array));
e_ang9= zeros(2,length(sr_array));
e_dir10= zeros(2,length(sr_array));
e_ang10= zeros(2,length(sr_array));
for ss = 1:length(sr_array)
snr = sr_array(ss);
%---------------7陣元---------------------------------------------------
sensor_num = 9;
d_angle = (0:sensor_num-1)'*2*pi/sensor_num;
x = R*cos(d_angle);y = R*sin(d_angle);
theta_array = [-9*d_angle(2)/8,d_angle(3)]*180/pi;
for it = 1:iteration
n = (randn(sensor_num,N_x)+1j*randn(sensor_num,N_x))/sqrt(2)*sqrt(10^(-45/10));
for dd = 1:2
theta = theta_array(dd);
tao = x*(sin(phi*pi/180).*cos(theta*pi/180))+y*(sin(phi*pi/180).*sin(theta*pi/180));
A = exp(-1j*2*pi*tao./l);
xt = A*(sqrt(10^(snr/10))*st)+n;
% -------------------2D-MUSIC算法-----------------------
Rxx = xt*xt'/N_x;
[U,S] = svd(Rxx);
[~,index] = sort(diag(S));
Un = U(:,index(1:sensor_num-est_num));
Gn = Un*Un';
Pmusic = zeros(length(dir9(dd,:)),length(ang));
for i=1:length(dir9(dd,:))
for k=1:length(ang)
a_tao = sin(ang(k))*(x*cos(dir9(dd,i))+y*sin(dir9(dd,i)));
a_theta = exp(-1j*2*pi*a_tao/l);
Pmusic(i,k)=1./abs((a_theta)'*Gn*a_theta);
end
end
[a,b]=find(Pmusic==max(max(Pmusic)));
aa = dir9(dd,a)*180/pi;
bb = ang(b)*180/pi;
e_dir9(dd,ss) = e_dir9(dd,ss)+(aa-theta)^2;
e_ang9(dd,ss) = e_ang9(dd,ss)+(bb-phi)^2;
end
end
%---------------8陣元---------------------------------------------------
sensor_num = 8;
d_angle = (0:sensor_num-1)'*2*pi/sensor_num;
x = R*cos(d_angle);y = R*sin(d_angle);
theta_array = [-9*d_angle(2)/8,2.5*d_angle(2)]*180/pi;
for it = 1:iteration
for dd = 1:2
theta = theta_array(dd);
tao = x*(sin(phi*pi/180).*cos(theta*pi/180))+y*(sin(phi*pi/180).*sin(theta*pi/180));
A = exp(-1j*2*pi*tao./l);
n = (randn(sensor_num,N_x)+1j*randn(sensor_num,N_x))/sqrt(2)*sqrt(10^(-45/10));
xt = A*(sqrt(10^(snr/10))*st)+n;
% -------------------2D-MUSIC算法-----------------------
Rxx = xt*xt'/N_x;
[U,S] = svd(Rxx);
[~,index] = sort(diag(S));
Un = U(:,index(1:sensor_num-est_num));
Gn = Un*Un';
Pmusic = zeros(length(dir10(dd,:)),length(ang));
for i=1:length(dir10(dd,:))
for k=1:length(ang)
a_tao = sin(ang(k))*(x*cos(dir10(dd,i))+y*sin(dir10(dd,i)));
a_theta = exp(-1j*2*pi*a_tao/l);
Pmusic(i,k)=1./abs((a_theta)'*Gn*a_theta);
end
end
[a,b]=find(Pmusic==max(max(Pmusic)));
aa = dir10(dd,a)*180/pi;
bb = ang(b)*180/pi;
e_dir10(dd,ss) = e_dir10(dd,ss)+(aa-theta)^2;
e_ang10(dd,ss) = e_ang10(dd,ss)+(bb-phi)^2;
end
end
end
figure;
subplot(121);
plot(sr_array+45,e_dir9(1,:)/iteration,'-^k',sr_array+45,e_dir9(2,:)/iteration,'-*k');
legend('-9/8 360/M','3 360/M');
grid on;%axis([-10,20,-.01,.45]);
colormap gray;
xlabel('信噪比/dB');ylabel('均方誤差/°');
title('9元陣方向角誤差');
%
% subplot(222);
% plot(sr_array+45,e_ang9(1,:)/iteration,'-^k',sr_array+45,e_ang9(2,:)/iteration,'-*k');
% legend('-9/8 360/M','3 360/M');
% grid on;axis([-10,20,-.1,1.45]);
% colormap gray;
% xlabel('信噪比/dB');ylabel('均方誤差/°');
% title('7元陣俯仰角誤差');
subplot(122);
plot(sr_array+45,e_dir10(1,:)/iteration,'-^k',sr_array+45,e_dir10(2,:)/iteration,'-*k');
legend('-5/4 360/M','5/2 360/M');
grid on;%axis([-10,20,-.01,.45]);
colormap gray;
xlabel('信噪比/dB');ylabel('均方誤差/°');
title('8元陣方向角誤差');
%
%
% subplot(224);
% plot(sr_array+45,e_ang10(1,:)/iteration,'-^k',sr_array+45,e_ang10(2,:)/iteration,'-*k');
% legend('-5/4 360/M','3 360/M');
% grid on;axis([-10,20,-.1,1.45]);
% colormap gray;
% xlabel('信噪比/dB');ylabel('均方誤差/°');
% title('8元陣俯仰角誤差');
6、
%=========================================================================
% UCA_multi_in_2D
%
%=========================================================================
clc;
clear all;
close all;
%------------------------常數表-------------------------------
c = 3e8;
namda = c/18e9;
est_num = 1;
iteration = 100;
sr_array = [-50,-47.5,-45,-42.5,-40,-35,-27.5];
phi = 60;
%% ----------------入射信號模型-------------------------------
N_x = 2^5; %快拍點數
F0 = 18e9; %中心頻率
B = 20e6; %帶寬
Fs = 2*B; %采樣頻率
Ts = 1/Fs; %采樣時間
T = (N_x-1)*Ts; %快拍持續時間
u = B/T; %頻率變化率
t = -T/2:Ts:T/2; %時間軸點
l = c/18e9;
st = exp(1j*2*pi*(F0*t+.5*u*t.^2));
dir7 =(46:.25:56)*pi/180; %(-50:.25:-40)*pi/180;
dir8 =(62.5:.25:72.5)*pi/180;
dir9 =(35:.25:45)*pi/180;%(-51:.25:-41)*pi/180;
dir10=(49:.25:59)*pi/180;
ang=(50:.25:70)*pi/180;
e_dir7 = zeros(1,length(sr_array));
e_dir8 = zeros(1,length(sr_array));
e_dir9 = zeros(1,length(sr_array));
e_dir10= zeros(1,length(sr_array));
R = 0.1;
for ss = 1:length(sr_array)
snr = sr_array(ss);
%--------------------7陣元---------------------------------------
sensor_num = 7;
d_angle = (0:sensor_num-1)'*2*pi/sensor_num;
x = R*cos(d_angle);y = R*sin(d_angle);
theta = d_angle(2)*180/pi;
for it = 1:iteration
n = (randn(sensor_num,N_x)+1j*randn(sensor_num,N_x))/sqrt(2)*sqrt(10^(-45/10));
tao = x*(sin(phi*pi/180).*cos(theta*pi/180))+y*(sin(phi*pi/180).*sin(theta*pi/180));
A = exp(-1j*2*pi*tao./l);
xt = A*(sqrt(10^(snr/10))*st)+n;
% -------------------2D-MUSIC算法-----------------------
Rxx = xt*xt'/N_x;
[U,S] = svd(Rxx);
[~,index] = sort(diag(S));
Un = U(:,index(1:sensor_num-est_num));
Gn = Un*Un';
Pmusic = zeros(length(dir7),length(ang));
for i=1:length(dir7)
for k=1:length(ang)
a_tao = sin(ang(k))*(x*cos(dir7(i))+y*sin(dir7(i)));
a_theta = exp(-1j*2*pi*a_tao/l);
Pmusic(i,k)=1./abs((a_theta)'*Gn*a_theta);
end
end
[a,b]=find(Pmusic==max(max(Pmusic)));
aa = dir7(a)*180/pi;
e_dir7(ss) = e_dir7(ss)+(aa-theta)^2;
end
%--------------------8陣元-------------------------------------
sensor_num = 8;
d_angle = (0:sensor_num-1)'*2*pi/sensor_num;
x = R*cos(d_angle);y = R*sin(d_angle);
theta = 1.5*d_angle(2)*180/pi;
for it = 1:iteration
tao = x*(sin(phi*pi/180).*cos(theta*pi/180))+y*(sin(phi*pi/180).*sin(theta*pi/180));
A = exp(-1j*2*pi*tao./l);
n = (randn(sensor_num,N_x)+1j*randn(sensor_num,N_x))/sqrt(2)*sqrt(10^(-45/10));
xt = A*(sqrt(10^(snr/10))*st)+n;
% -------------------2D-MUSIC算法-----------------------
Rxx = xt*xt'/N_x;
[U,S] = svd(Rxx);
[~,index] = sort(diag(S));
Un = U(:,index(1:sensor_num-est_num));
Gn = Un*Un';
Pmusic = zeros(length(dir8),length(ang));
for i=1:length(dir8)
for k=1:length(ang)
a_tao = sin(ang(k))*(x*cos(dir8(i))+y*sin(dir8(i)));
a_theta = exp(-1j*2*pi*a_tao/l);
Pmusic(i,k)=1./abs((a_theta)'*Gn*a_theta);
end
end
[a,b]=find(Pmusic==max(max(Pmusic)));
aa = dir8(a)*180/pi;
e_dir8(ss) = e_dir8(ss)+(aa-theta)^2;
end
%---------------9陣元---------------------------------------------------
sensor_num = 9;
d_angle = (0:sensor_num-1)'*2*pi/sensor_num;
x = R*cos(d_angle);y = R*sin(d_angle);
theta = d_angle(2)*180/pi;
for it = 1:iteration
n = (randn(sensor_num,N_x)+1j*randn(sensor_num,N_x))/sqrt(2)*sqrt(10^(-45/10));
tao = x*(sin(phi*pi/180).*cos(theta*pi/180))+y*(sin(phi*pi/180).*sin(theta*pi/180));
A = exp(-1j*2*pi*tao./l);
xt = A*(sqrt(10^(snr/10))*st)+n;
% -------------------2D-MUSIC算法-----------------------
Rxx = xt*xt'/N_x;
[U,S] = svd(Rxx);
[~,index] = sort(diag(S));
Un = U(:,index(1:sensor_num-est_num));
Gn = Un*Un';
Pmusic = zeros(length(dir9),length(ang));
for i=1:length(dir9)
for k=1:length(ang)
a_tao = sin(ang(k))*(x*cos(dir9(i))+y*sin(dir9(i)));
a_theta = exp(-1j*2*pi*a_tao/l);
Pmusic(i,k)=1./abs((a_theta)'*Gn*a_theta);
end
end
[a,b]=find(Pmusic==max(max(Pmusic)));
aa = dir9(a)*180/pi;
e_dir9(ss) = e_dir9(ss)+(aa-theta)^2;
end
%---------------10陣元---------------------------------------------------
sensor_num = 10;
d_angle = (0:sensor_num-1)'*2*pi/sensor_num;
x = R*cos(d_angle);y = R*sin(d_angle);
theta = 1.5*d_angle(2)*180/pi;
for it = 1:iteration
tao = x*(sin(phi*pi/180).*cos(theta*pi/180))+y*(sin(phi*pi/180).*sin(theta*pi/180));
A = exp(-1j*2*pi*tao./l);
n = (randn(sensor_num,N_x)+1j*randn(sensor_num,N_x))/sqrt(2)*sqrt(10^(-45/10));
xt = A*(sqrt(10^(snr/10))*st)+n;
% -------------------2D-MUSIC算法-----------------------
Rxx = xt*xt'/N_x;
[U,S] = svd(Rxx);
[~,index] = sort(diag(S));
Un = U(:,index(1:sensor_num-est_num));
Gn = Un*Un';
Pmusic = zeros(length(dir10),length(ang));
for i=1:length(dir10)
for k=1:length(ang)
a_tao = sin(ang(k))*(x*cos(dir10(i))+y*sin(dir10(i)));
a_theta = exp(-1j*2*pi*a_tao/l);
Pmusic(i,k)=1./abs((a_theta)'*Gn*a_theta);
end
end
[a,b]=find(Pmusic==max(max(Pmusic)));
aa = dir10(a)*180/pi;
e_dir10(ss) = e_dir10(ss)+(aa-theta)^2;
end
end
figure;
plot(sr_array+45,e_dir7/iteration,'-^k',sr_array+45,e_dir8/iteration,'-*k');hold on;
plot(sr_array+45,e_dir9/iteration,'-sk',sr_array+45,e_dir10/iteration,'-dk');
legend('7陣元','8陣元','9陣元','10陣元');
grid on;axis([-5.5,18,-.02,.25]);
colormap gray;
xlabel('信噪比/dB');ylabel('均方誤差/°');
title('不同陣元相同半徑測向試驗');
7、
%=========================================================================
% Circular Array Classical-Music
%
%=========================================================================
clc;
clear all;
close all;
c = 3e8;
phi = 60;
namda = c/18e9;
R = 7.5/100;
R = 10/100;
snr = -35; %信噪比
N_x = 2^5; %快拍點數
F0 = 18e9; %中心頻率
B = 20e6; %帶寬
Fs = 40e6; %采樣頻率
Ts = 1/Fs; %采樣時間
T = (N_x-1)*Ts; %快拍持續時間
u = B/T; %頻率變化率
t = -T/2:Ts:T/2; %時間軸點
l = c/F0;
st = sqrt(10^(snr/10))*exp(1j*2*pi*(F0*t+.5*u*t.^2));
%% -----------------9------------------------
sensor_num = 9;
R = 7.1239/100;
d_angle = (0:sensor_num-1)'*2*pi/sensor_num;
theta =d_angle(2)*180/pi;
x = R*cos(d_angle);y = R*sin(d_angle);
tao = sin(phi*pi/180)*(x*cos(theta*pi/180)+y*sin(theta*pi/180));
A = exp(-1j*2*pi*tao/l);
n = randn(sensor_num,N_x)*sqrt(10^(-45/10));
xt = A*st+n;
% -------------------2D-MUSIC算法-----------------------
Rx = xt*xt'/N_x;
[U,S] = eig(Rx);
est_sour = 1;
[~,index] = sort(diag(S));
Un = U(:,index(1:sensor_num-est_sour));%*diag([0.05,50,3,1,0.001,1000,777]);
Gn = Un*Un';
dir=(-180:.25:179.8)*pi/180;
ang=(20:.25:91)*pi/180;
Pmusic9 = zeros(length(dir),length(ang));
for i=1:length(dir)
for k=1:length(ang)
a_tao = sin(ang(k))*(x*cos(dir(i))+y*sin(dir(i)));
a_theta = exp(-1j*2*pi*a_tao/l);
Pmusic9(i,k)=1./abs((a_theta)'*Gn*a_theta);
end
end
P_music9 = 10*log10(Pmusic9/min(min(Pmusic9)));
%% -----------------9------------------------
sensor_num = 10;
R = 4.5879/100;
%R = 10/100;
d_angle = (0:sensor_num-1)'*2*pi/sensor_num;
theta =1.6*d_angle(2)*180/pi;
x = R*cos(d_angle);y = R*sin(d_angle);
tao = sin(phi*pi/180)*(x*cos(theta*pi/180)+y*sin(theta*pi/180));
A = exp(-1j*2*pi*tao/l);
n = randn(sensor_num,N_x)*sqrt(10^(-45/10));
xt = A*st+n;
% -------------------2D-MUSIC算法-----------------------
Rx = xt*xt'/N_x;
[U,S] = eig(Rx);
disp(est_sour);
[~,index] = sort(diag(S));
Un = U(:,index(1:sensor_num-est_sour));%*diag([0.05,50,3,1,0.001,1000,777]);
Gn = Un*Un';
dir=(-180:.25:179.8)*pi/180;
ang=(20:.25:91)*pi/180;
Pmusic10 = zeros(length(dir),length(ang));
for i=1:length(dir)
for k=1:length(ang)
a_tao = sin(ang(k))*(x*cos(dir(i))+y*sin(dir(i)));
a_theta = exp(-1j*2*pi*a_tao/l);
Pmusic10(i,k)=1./abs((a_theta)'*Gn*a_theta);
end
end
P_music10 = 10*log10(Pmusic10/min(min(Pmusic10)));
figure;
% subplot(221);
% [xx,yy] = meshgrid(ang*180/pi,dir*180/pi);
% mesh(xx,yy,P_music9);
% title('9元陣二維空間譜');
% xlabel('俯仰角/°');ylabel('方向角/°');zlabel('空間譜/dB');
% axis([20,91,-180,180,0,24]);%colormap gray;
subplot(121);
[xx,yy] = meshgrid(ang*180/pi,dir*180/pi);
mesh(xx,yy,P_music9);
title('9元陣方向角空間譜');
xlabel('俯仰角/°');ylabel('方向角/°');zlabel('空間譜/dB');
axis([20,91,-180,180,0,24]);%colormap gray;
% subplot(222);
% [xx,yy] = meshgrid(ang*180/pi,dir*180/pi);
% mesh(xx,yy,P_music10);
% title('10元陣二維空間譜');
% xlabel('俯仰角/°');ylabel('方向角/°');zlabel('空間譜/dB');
% axis([20,91,-180,180,0,24]);%colormap gray;
subplot(122);
[xx,yy] = meshgrid(ang*180/pi,dir*180/pi);
mesh(xx,yy,P_music10);
title('10元陣方向角空間譜');
xlabel('俯仰角/°');ylabel('方向角/°');zlabel('空間譜/dB');
axis([20,91,-180,180,0,24]);%colormap gray;
8、
% 波程差圖
clear all;
close all;
clc;
sensor_num = 9;
c = 3e8;
d_angle = (0:sensor_num-1)'*2*pi/sensor_num;
phi = 90;
namda = c/18e9;
if mod(sensor_num,2);
R = namda/(1-cos(d_angle(2)));%相鄰最小間距最大值
%R = namda/(cos(floor(sensor_num/4)*d_angle(2))+cos(d_angle(2)*floor((sensor_num+1)/4)-.5));%相鄰最大間距最大值
%R = namda/(1+cos(d_angle(2)/2));%任意陣元間距最大
else
R = namda/(2*sin(d_angle(2))*sin(0.5*d_angle(2)));%相鄰最小間距最大值
%R = namda/(2*sin(d_angle(2)/2));%相鄰最大間距最大值
%R = namda/2;%任意陣元間距最大值
end
R = 1;
dir = (0:.4:d_angle(2)*180/pi)*pi/180;
%idx=nchoosek(1:sensor_num,2); % 16取2的組合
idx1 = [1,2;2,9;9,3;3,8;8,4;4,7;7,5;5,6];
idx2 = [2,1;1,3;3,9;9,4;4,8;8,5;5,7;7,6];
r = zeros(length(dir),size(idx1,1));
x = R*cos(d_angle);y = R*sin(d_angle);
for i = 1:50
tao = (x*cos(dir(i))*sin(90*pi/180)+y*sin(dir(i))*sin(90*pi/180));
r(i,:) = -diff(tao(idx1),1,2);
end
for i = 51:101
tao = (x*cos(dir(i))*sin(90*pi/180)+y*sin(dir(i))*sin(90*pi/180));
r(i,:) = -diff(tao(idx2),1,2);
end
figure;
subplot(211);
plot(dir(1:5:length(dir))*180/pi,0.234*ones(1,length(dir(1:5:length(dir)))),'.k',dir*180/pi,r(:,end),'-.k',dir*180/pi,r(:,1:end-1),'k');grid on;
legend('0.234','k=1','k=其他');
title('(a)9元陣投影間隔')
xlabel('入射方向\°');ylabel('投影間隔/R');
axis([0,40,-0.01,0.7]);
% ---------------------------------------------------------------------
sensor_num = 10;
d_angle = (0:sensor_num-1)'*2*pi/sensor_num;
dir = linspace(0,d_angle(2),101);
%idx=nchoosek(1:sensor_num,2); % 16取2的組合
idx1 = [1,2;2,10;10,3;3,9;9,4;4,8;8,5;5,7;7,6];
idx2 = [2,1;1,3;3,10;10,4;4,9;9,5;5,8;8,6;6,7];
r = zeros(length(dir),size(idx1,1));
x = R*cos(d_angle);y = R*sin(d_angle);
for i = 1:50
tao = (x*cos(dir(i))*sin(90*pi/180)+y*sin(dir(i))*sin(90*pi/180));
r(i,:) = -diff(tao(idx1),1,2);
end
for i = 51:101
tao = (x*cos(dir(i))*sin(90*pi/180)+y*sin(dir(i))*sin(90*pi/180));
r(i,:) = -diff(tao(idx2),1,2);
end
subplot(212);
plot(dir(1:5:length(dir))*180/pi,0.3633*ones(1,length(dir(1:5:length(dir)))),'.k',dir*180/pi,r(:,2),'-.k',dir*180/pi,r(:,[1,3:6,7]),'k');grid on;
legend('0.3633','k=1','k=其他');
title('(b)10元陣投影間隔');
xlabel('入射方向/°');ylabel('投影間隔/R');
axis([0,36,-0.01,0.7]);