%classic music
clear all
clc
tic
%参数设定
M=8;
f=1000;c=1500;
lambda=c/f;d=lambda/2;
SNR=[5:2:25];
F=length(SNR);
N=128;%快拍数。
K=100;
doa=[-10 30];
P=length(doa);
for g=1:F %信噪比的循环。
for k=1:K
%阵列流型A
for i=1:P
A(:,i)=exp(-j*2*pi*d*[0:M-1]'/lambda*sin(doa(i)*pi/180));%均匀线性阵的阵列流型矢量。
end
%信源模型建立
for i=1:P
S(i,:)=sqrt(10^(SNR(g)/10))*(randn(1,N)+j*randn(1,N));%sqrt(10.^(snr/10))使信号符合该信噪比。
end
%通道幅度相位误差
for m=2:M
u1= unifrnd(-pi,pi,1,M-1);%1*7.
ph(1,1)=1;
ph(:,m)=exp(j*0.2*u1(1,m-1));%1*8相位误差。
u2=unifrnd(-1,1,1,M-1);%1*7。
al(1,1)=1;
al(:,m)=1+0.2*u2(1,m-1);%1*8幅度误差。
gamma(:,1)=1;
gamma(:,m)=ph(:,m)*al(:,m);%1*8幅相误差。
end
%接收信号模型建立
X=A*S+1/sqrt(2)*(randn(M,N)+j*randn(M,N));%1/sqrt(2)使高斯白噪声的能量为1。
%协方差矩阵特征值分解得到噪声子空间
R=X*X'/N;
[V,D]=eig(R);% [V,D]=eig(A)：求矩阵R的全部特征值，构成对角阵D，并求A的特征向量构成V的列向量。
[Y,I]=sort(diag(D));%diag是将对角阵的对角元素提取成一个向量.如果A是向量，sort(A) 对A中元素按照升序排列。
% 如果A是矩阵sort(A) 对A按每一列元素按照升序排列。Y为排列后矩阵，I是原矩阵中各元素位置所组成的新矩阵。
Un=V(:,I(1:M-P));%把原矩阵第1到M-P列中的数对应到V矩阵中的列数，建立新的矩阵赋值给Un。
Xe=diag(gamma)*A*S+1/sqrt(2)*(randn(M,N)+j*randn(M,N));%1/sqrt(2)使高斯白噪声的能量为1。
%协方差矩阵特征值分解得到噪声子空间
Re=Xe*Xe'/N;
[Ve,De]=eig(Re);% [V,D]=eig(A)：求矩阵R的全部特征值，构成对角阵D，并求A的特征向量构成V的列向量。
[Ye,Ie]=sort(diag(De));%diag是将对角阵的对角元素提取成一个向量.如果A是向量，sort(A) 对A中元素按照升序排列。
% 如果A是矩阵sort(A) 对A按每一列元素按照升序排列。Y为排列后矩阵，I是原矩阵中各元素位置所组成的新矩阵。
Une=Ve(:,Ie(1:M-P));%把原矩阵第1到M-P列中的数对应到V矩阵中的列数，建立新的矩阵赋值给Un。
%谱峰搜索部分
theta=-90:0.1:90;%线阵的搜索范围为-90~90度
for i=1:length(theta)
a_theta=exp(-j*(0:M-1)'*2*pi*d*sin(pi/180*theta(i))/lambda);
Pmusic(i)=1./abs((a_theta)'*Un*Un'*a_theta);
Pmusic2(i)=1./abs((a_theta)'*Une*Une'*a_theta);
P1(i)=10*log10(Pmusic(i)/min(Pmusic));
P2(i)=10*log10(Pmusic2(i)/min(Pmusic2));
end
[peaks,locals]=findpeaks(Pmusic);%峰值及其位置
b=sortrows( [peaks' locals']);%峰值按从小到大排列
bb=b';
B=bb(:,size(b',2)-1:size(b',2));
BB=B(2,:);%峰值角度对应的数字
BBB=sort(BB);%最大峰值角度对应的数字
M_theta=theta(BBB);%所测得入射角
mistake1(g,2*k-1:2*k)=M_theta-doa;%每次实际入射角与测得入射角得差值，11*200,11对应11个信噪比，200对应200个快拍。
[peaks2,locals2]=findpeaks(Pmusic2);%峰值及其位置
b2=sortrows( [peaks2' locals2']);%峰值按从小到大排列
bb2=b2';
B2=bb2(:,size(b2',2)-1:size(b2',2));
BB2=B2(2,:);%峰值角度对应的数字
BBB2=sort(BB2);%最大峰值角度对应的数字
M_theta2=theta(BBB2);%所测得入射角
mistake2(g,2*k-1:2*k)=M_theta2-doa;%每次实际入射角与测得入射角得差值
end
RMSE(g,:)=sqrt(sum(mistake1(g,:).^2,2)/2/K);
RMSE2(g,:)=sqrt(sum(mistake2(g,:).^2,2)/2/K);
end
figure(1)
plot(SNR,RMSE,'bd-',SNR,RMSE2,'ro-','MarkerFaceColor','c');
xlabel('SNR/dB');
ylabel('RMSE/°');
legend('无误差均方值RMSE1','有误差均方值RMSE2');
title('存在幅度相位误差与无幅相误差的角度误差均方值比较');
grid on;
figure(2)
plot(theta,P1,'b--',theta,P2,'r-');
axis([-60 60 0 60]);
xlabel('theta/°');
ylabel('功率谱/dB');
%h=legend('g=0','g=0.2',0);
%h=legend('g=0','g=0.2',0);
h=legend('理想情况下','存在幅相误差',0);
set(h,'box','off');
title('存在相位误差与理想情况下的功率谱比较');
%title('仅存在幅度误差与理想情况下的功率谱比较');
%title('存在幅相误差与理想情况下的功率谱比较');
grid on;
toc
