ai5_pian0807.rar_MUSIC二维_二维music算法_信号超分辨_超分辨 MUSIC_超平面

%%%6 sensor. 4-8GHz%%% clc; clear; clear; sm=6; % 为天线得个数 snr=20; % 仿真时候用得信噪比 nd=200; % 采样个数 md=181; snr=10^(snr/20); f=5*10^3*10^6; % 载频码 dataf=10^6; wl=3*10^8/f; wl1=3*10^8/(f+dataf); diz=[25 0 0 0 0 0]/1000; dix=[-10 5.8 -10 -5.8 -10 4.1]/100; %参数都为距离不是角度 单位是厘米 diy=[10 10 5.8 10 -5.8 -4.1]/100; fw=[50 70]; %fangwei fy=[70 80]; %fuyang% rkx=fw*pi/180; %fangwei% rky=fy*pi/180; %fuyang% ll=length(rkx); p31=rand(1,1); for m=1:sm; for n=1:ll; %yxl a(m,n)=exp(-j*2*pi/wl*((dix(m)-dix(1))*sin(rkx(n))*cos(rky(n))+(diy(m)-diy(1))*sin(rkx(n))*sin(rky(n))+diz(m)*cos(rkx(n)))); a(m,n)=exp(-j*2*pi/wl*((dix(m))*cos(rkx(n))*cos(rky(n))+(diy(m))*sin(rkx(n))*cos(rky(n))+diz(m)*sin(rky(n)))); rand(1) end end for m=1:ll; for n=1:nd; p1=rand(1,1); p2=rand(1,1); sr(m,n)=sqrt(-2*snr*snr*log(p1))*cos(2*pi*p2); si(m,n)=sqrt(-2*snr*snr*log(p1))*sin(2*pi*p2); s(m,n)=sr(m,n)+j*si(m,n); end end for m=1:sm; for n=1:nd; p3=rand(1,1); p4=rand(1,1); wr(m,n)=sqrt(-2*log(p3))*cos(2*pi*p4); wi(m,n)=sqrt(-2*log(p3))*sin(2*pi*p4); w(m,n)=wr(m,n)+j*wi(m,n); end end %creat output signal cc1=cputime; x=a*s+w % X为6×100个数据 实际上工程整个数据由数字接收机得FIFO提供 x' r=x*x'/nd % r 为协方差矩阵 [z,d]=eig(r,'nobalance') d1=diag(d) ; %d1为列向量，各个元素为d的主对角线元素值，diag为建立或提取对角阵 [lambtal,k]=sort(d1); %将d1元素按照升序排列付给lambtal；k对应的是lambtal的元素对应在d1列向量中的位置 lambtal=(fliplr(lambtal'))' ; %fliplr:左右方向反转数组；此语句的功能是将特征值按照降幂排列 z=z(:,k) ; %特征向量按照升幂特征值的位置对应 z=fliplr(z) ; %特征向量按照降幂特征值的位置对应 for p=1:sm; c1=1; c2=0; for i=sm-p+1:sm; c1=c1*lambtal(i); c2=c2+lambtal(i); end c1=c1^(1.0/p); c2=c2/p; lr(p)=(c1/c2)^(nd*p/2); mdl(p)=(sm-p)*(sm+p+1)*log10(nd)/4-log10(lr(p)); end [lambta_min ,nx]=min(mdl); d=sm-nx; % MDL准则求信源个数 en=z(:,d+1:sm); % en为噪声子空间 x=linspace(0,360,md); rka=x*pi/180; %方位 y=linspace(0,90,md); rkb=y*pi/180; %俯仰 for k=1:sm; for m=1:md; for n=1:md; %yxl b(k,m,n)=exp(-j*2*pi/wl1*((dix(k)-dix(1))*sin(rka(m))*cos(rkb(n))+(diy(k)-diy(1))*sin(rka(m))*sin(rkb(n))+diz(k)*cos(rka(m)))); b(k,m,n)=exp(-j*2*pi/wl*((dix(k))*cos(rka(m))*cos(rkb(n))+(diy(k))*sin(rka(m))*cos(rkb(n))+diz(k)*sin(rkb(n)))); end end end for i=1:md; for j=1:md; pmu(i,j)=(b(:,i,j))'*en*(en)'*b(:,i,j); % pmu(i,j) 对应得是空间每一点得谱峰值 end end pmu=1./pmu; pmu=abs(pmu); pmu=20*log(pmu); plot3(x,y,pmu); mesh(x,y,pmu);%mesh可画出立体网状图，plot可画出立体曲面图。 hidden off; %meshc(x,y/4,pmu); qx=0; qy=0; if pmu(1,1)>=pmu(1,2); if pmu(1,1)>=pmu(2,1); qx=qx+1; qy=qy+1; pm(qx,qy)=pmu(i,j); pm(qx,qy)=pmu(1,1); pda(qx)=1; pdb(qy)=1; end end for k=2:md-1; if (pmu(1,k-1)<=pmu(1,k)); if (pmu(1,k)>=pmu(1,k+1)); if (pmu(1,k)>=pmu(2,k)); qx=qx+1; qy=qy+1; pm(qx,qy)=pmu(i,j); pda(qx)=1; pdb(qy)=j; end end end end for i=2:md-1; if (pmu(i,1)>=pmu(i-1,1)); if (pmu(i,1)>=pmu((i+1),1)); if (pmu(i,1)>=pmu(i,2)); qx=qx+1; qy=qy+1; pm(qx,qy)=pmu(i,1); pda(qx)=i; pdb(qy)=1; end end end for j=2:md-1; if(pmu(i-1,j)<=pmu(i,j)); if (pmu(i,j)>=pmu(i+1,j)); if(pmu(i,j-1)<=pmu(i,j)); if (pmu(i,j)>=pmu(i,j+1)); qx=qx+1; qy=qy+1; pm(qx,qy)=pmu(i,j); pda(qx)=i; pdb(qy)=j; end end end end end end aa=qx; bb=qy; cc2=cputime; % pm=diag(pm); %yxl pm=diag(pm) %pm1=diag(pm); hold on; % [qq]=find(pm>140); %yxl [qq]=find(pm>0) %qq为大于0的数在pm中的位置 pda=pda(qq); pdb=pdb(qq); pm1=pm(qq); hold on; [pm2,k]=sort(pm1); %将pm1元素按照升序排列付给pm2；k对应的是pm1的元素对应在pm2列向量中的位置 pm3=(fliplr(pm2'))' %fliplr:左右方向反转数组；此语句的功能是将特征值按照降幂排列 pm4=pm3(1:d); %取出前d个最大的数 pda=pda(:,k) ; %pda按照升幂值的位置对应 pda=fliplr(pda); pda=pda(1:(d)); pdb=pdb( :,k) ; %pda按照升幂值的位置对应 pdb=fliplr(pdb) ; pdb=pdb(1:(d)); % plot3((pdb-1)*2,(pda-1)/2,pm1,'*m'); % plot3((pdb-1)*2,(pda-1)/2,pm1,'*m'); % xlabel(' 空间方位（0\circ\sim 360\circ）') % ylabel(' 俯仰方位（0\circ\sim 90\circ）') %plot3((pdb-1),(pda-1)/4,pm4,'*m'); xlabel(' 方位（0\circ\sim 360\circ）') ylabel(' 俯仰（0\circ\sim 90\circ）') zlabel(' 谱密度') hold off; title('6 sensors random array') grid on; pm4 pda=(pda-1)*2 pdb=(pdb-1)/2 cc2=cputime; c=cc2-cc1;