模式空间法_模式空间算法_模式空间法_圆阵模式空间_DOA_DOA估计_

clc; clear; sig_num = 1000;%%信号的快拍数 snr=10; rou=0; x1=randint(sig_num,1,4);%产生sig_sum*1的矩阵，数值在0-3之间 x2=randint(sig_num,1,4); x3=randint(sig_num,1,4); x222=rou*x1+sqrt(1-(abs(rou))^2)*x3;%产生相关的两路信号 x2=round(x222); for i=1:sig_num if x2(i)>=4 x2(i)=3; elseif x2(i)<=0 x2(i)=0; end end y1=pskmod(x1,4,pi/4); y1=y1(1:sig_num); y2=pskmod(x2,4,pi/4); y2=y2(1:sig_num); x=[y1';y2']; s=x; q=2; p=8;%智能天线阵元数 h=3;%最大模式数(这里p>2*h即可） P=2*h+1; fc=2000*10^6;%*(10^6); %hz lemda=3*10^8/fc; banjing=lemda/(4*sin(pi/p)); jh=pi./180;%角度转弧度 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %均匀圆阵转化为均匀线阵 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% alpha=[10,30]; theta=[20,40]; kesei=2*pi*banjing./lemda;%% j=sqrt(-1);%虚数单位 n=[0:(p-1)]'; gama=2*pi*n/p; for i=1:q for n=0:(p-1) a(n+1,i)=exp(j*kesei*sin(theta(i)*jh)*cos(alpha(i)*jh-gama(n+1))); end end x11=a*s; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% vn=var(x11.')./(10^(snr/10));%var(),求样本方差无偏估计值 n11=randn(p,sig_num)+j*randn(p,sig_num); a1=(vn./var(n11.')); b1=((n11.'-ones(size(n11,2),1)*mean(n11.'))).'; for i=1:p c1(i,:)=a1(i)*b1(i,:); end xn=x11+c1; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%空间虚拟线阵%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for i=-h:h%-3到3，h=3 cv1(i+h+1)=j^(-abs(i));% end cv=diag(cv1);%%%%%构造矩阵J %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%p=8 for m=-h:h for i=0:p-1 w(i+1,m+h+1)=exp(-j*m*gama(i+1));%构造wq end end V=1./sqrt(p)*w;%wq Fe=cv*V'; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for i=-h:h alph(i+h+1)=2*pi*i/P; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for i=-h:h for m=-h:h v(m+h+1,i+h+1)=exp(j*2*pi*i*m./P); end end W=(1/sqrt(P))*v; Fr=W'*Fe; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %music算法实现过程 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% rx=xn*xn'./sig_num;%计算协方差矩阵 ry=Fr*rx*Fr'; r=real(ry); % r=Fe*rx*Fe'; [u,sig,v]=svd(r);%对矩阵进行分解，U和V代表二个相互正交矩阵，而Sig代表一对角矩阵 signalspace=u(:,1); noisespace=u(:,3:7); vv=conj(noisespace)*noisespace';%conj计算复数共轭值,共轭*转置 % [eigenvector,D]=eig(r,'nobalance'); % R_eigenvalue=diag(abs(D))'; % [eigenvalue,index]=sort(R_eigenvalue); % R_eigenvalue_decreasing=fliplr(eigenvalue); % eigenvector=eigenvector(:,index); % noisespace1=fliplr(eigenvector); % noisespace=noisespace1(:,q+1:2*h+1); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %music功率谱 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% bbb=1; for alpha1=0:100 kk=1; for theta1=0:90 for n=0:(p-1) aa(n+1)=exp(j*kesei*sin(theta1*jh)*cos(alpha1*jh-gama(n+1))); end bb=Fr*aa.'; ww=(bb.')*noisespace*(noisespace.')*bb; smusic(kk,bbb)=abs(1./ww); kk=kk+1; end bbb=bbb+1; end b=imregionalmax(smusic);%BW=imregionalmax(I);返回识别I中的区域最大值的二值图像BW [I,J]=find(b); alpha1=0:100; theta1=0:90; figure(1); contour(alpha1,theta1,10*log10(smusic),2,'bd'); xlabel('方位角/°'); ylabel('俯仰角/°'); grid on;