MUSIC.zip_coprime MUSIC_officerstz_互质MUSIC_互质阵 music_互质阵MUSIC、

%翟会2016.10.14 非圆互质线阵MUSIC %基于吕卫华的降维方法 clear,clc M = 3; N = 4; K = 3; L = 100; dm = N/2; dn = M/2; ds = 0.1; ddm = 0:dm:(M-1)*dm; ddn = 0:dn:(N-1)*dn; % theta = -90+180*rand(1,2) theta = [0,10,25]; alpha = [5 10 15]; snr = 20; Am = exp(-j*2*pi*ddm.'*sin(theta*pi/180)); An = exp(-j*2*pi*ddn.'*sin(theta*pi/180)); Alpha = diag(exp(-j*alpha*pi/180)); %NC S = randn(K,L); X0m = Am*Alpha*S;%NC X0n = An*Alpha*S;%NC Xm = awgn(X0m,snr,'measured'); Xn = awgn(X0n,snr,'measured'); Ym = [Xm;conj(Xm)]; Yn = [Xn;conj(Xn)]; Rym = Ym*Ym'/2*L; Ryn = Yn*Yn'/2*L; [EVm,Dm] = eig(Rym); [EVAm,Im] = sort(diag(Dm).'); EVm = fliplr(EVm(:,Im)); Enm = EVm(:,K+1:end); %噪声子空间 em=[1;0]; [EVn,Dn] = eig(Ryn); [EVAn,In] = sort(diag(Dn).'); EVn = fliplr(EVn(:,In)); Enn = EVn(:,K+1:end); en=[1;0]; c = 0; for tt=-90:ds:90 c = c+1; thet(c) = tt; am = exp(-j*2*pi*ddm.'*sin(tt*pi/180)); pm = [am,(linspace(0,0,M)).']; pm = [pm;fliplr(conj(pm))]; Qm = pm'*Enm*Enm'*pm; invQm = pinv(Qm); fm(c) = em'*invQm*em; an = exp(-j*2*pi*ddn.'*sin(tt*pi/180)); pn = [an,(linspace(0,0,N)).']; pn = [pn;fliplr(conj(pn))]; Qn = pn'*Enn*Enn'*pn; invQn = pinv(Qn); fn(c) = en'*invQn*en; end fm_am = abs(fm); %子阵M谱函数 fn_am = abs(fn); %子阵N谱函数 plot(thet,fm_am,'-b'); hold on plot(thet,fn_am,'-r'); title('NC CLA MUSIC'); xlabel('theta/度'); ylabel('幅值'); legend(['M=',num2str(M)],['N=',num2str(N)]); hold off; length_m = length(fm_am); length_n = length(fn_am); tm = 0; for mm=2:length_m-1 if fm_am(mm)>fm_am(mm-1) && fm_am(mm)>fm_am(mm+1) tm = tm+1; m_temp(tm,1) = -90+ds*(mm-1); %M峰值横坐标 m_temp(tm,2) = fm_am(mm); %M峰值纵坐标 end end tn = 0; for nn=2:length_n-1 if fn_am(nn)>fn_am(nn-1) && fn_am(nn)>fn_am(nn+1) tn = tn+1; n_temp(tn,1) = -90+ds*(nn-1); n_temp(tn,2) = fn_am(nn); end end [m_temp0,Im] = sort(m_temp(:,2)); %按行由小到大排列 m_temp = flipud(m_temp(Im,:)); [n_temp0,In] = sort(n_temp(:,2)); n_temp = flipud(n_temp(In,:)); ang_m = m_temp(1:K*N,:); ang_n = n_temp(1:K*M,:); % ang_m = pi*sin(ang_m*pi/180); % ang_n = pi*sin(ang_n*pi/180); [ang_m1,Im]= sort(ang_m(:,1)); %ang_m按照角度由小到大排序 ang_m = ang_m(Im,:); [ang_n1,In]= sort(ang_n(:,1)); %ang_n按照角度由小到大排序 ang_n = ang_n(In,:); %%%%%%寻找三个最接近的谱峰 for nn=1:K*N for mm=1:K*M c1(nn,mm) = abs(ang_m(nn)-ang_n(mm)); end end %%%寻找第一个幅度大于0.5的谱峰 while 1 [C_min11,I_rank1] = min(c1); [C_min12,I_colu1] = min(C_min11); ttm1 = ang_m(I_rank1(I_colu1),1); %%M阵列对应的角度 ttn1 = ang_n(I_colu1,1); %%N阵列对应的角度 famm1 = ang_m(I_rank1(I_colu1),2); %%%若对应的角度，谱峰的值过小，舍弃该角度，并重新搜索下一个更接近的角度 famn1 = ang_n(I_colu1,2); if((famm1>0.5)&&(famn1>0.5)) %当snr=-5时，正确的角度谱峰均远大于0.5 angle_doa1 = (ttm1+ttn1)/2; break; else %%%去掉该错误角度所在的行和列的影响 c1(:,I_colu1)=100; %将第一个角度所在的行列变为一个较大的数，这里用100，避免下次再搜索到它 c1(I_rank1(I_colu1),:)=100; end end %%%去掉第一个角度所在的行和列的影响 c2 = c1; c2(:,I_colu1)=100; %将第一个角度所在的行列变为一个较大的数，这里用100，避免下次再搜索到它 c2(I_rank1(I_colu1),:)=100; %%%寻找第二个幅度大于0.5的谱峰 while 1 [C_min21,I_rank2] = min(c2); [C_min22,I_colu2] = min(C_min21); ttm2 = ang_m(I_rank2(I_colu2)); ttn2 = ang_n(I_colu2); famm2 = ang_m(I_rank2(I_colu2),2); %%%若对应的角度，谱峰的值过小，舍弃该角度，并重新搜索下一个更接近的角度 famn2 = ang_n(I_colu2,2); if((famm2>0.5)&&(famn2>0.5)) %当snr=-5时，正确的角度谱峰均远大于0.5 angle_doa2 = (ttm2+ttn2)/2; break; else %%%去掉该错误角度所在的行和列的影响 c2(:,I_colu2)=100; %将第二个角度所在的行列变为一个较大的数，这里用100，避免下次再搜索到它 c2(I_rank2(I_colu2),:)=100; end end %%%去掉第二个角度所在的行和列的影响 c3 = c2; c3(:,I_colu2)=100; %将第二个角度所在的行列变为一个较大的数，这里用100，避免下次再搜索到它 c3(I_rank1(I_colu2),:)=100; %%%寻找第三个幅度大于0.5的谱峰 while 1 [C_min31,I_rank3] = min(c3); [C_min32,I_colu3] = min(C_min31); ttm3 = ang_m(I_rank3(I_colu3)); ttn3 = ang_n(I_colu3); famm3 = ang_m(I_rank3(I_colu3),2); %%%若对应的角度，谱峰的值过小，舍弃该角度，并重新搜索下一个更接近的角度 famn3 = ang_n(I_colu3,2); if((famm3>0.5)&&(famn3>0.5)) %当snr=-5时，正确的角度谱峰均远大于0.5 angle_doa3 = (ttm3+ttn3)/2; break; else %%%去掉该错误角度所在的行和列的影响 c3(:,I_colu3)=100; %将第一个角度所在的行列变为一个较大的数，这里用100，避免下次再搜索到它 c3(I_rank1(I_colu3),:)=100; end end doa = sort([angle_doa1 angle_doa2 angle_doa3])