24小时热门版块排行榜

返回列表

影月水澜

新虫 (小有名气)

应助: 0 (幼儿园)
金币: 358.5
红花: 1
帖子: 53
在线: 61.1小时
虫号: 1826817
注册: 2012-05-20
专业: 敏感电子学与传感器

[求助] 程序运行结果不对，自己查不出来错，请各位大侠帮忙

下面是写的一段UCA ROOT-MUSIC程序，但是运行结果不对，不是设置的[20,50]。同样是本程序的噪声子空间，用谱峰搜索的MUSIC就是对的，但是用ROOT-MUSIC就不对。应该就是root-music那段程序有错，但是自己查也查不出来错误，请各位大侠帮帮忙，非常感谢！

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=dmod(x1,1,0.1,4,'psk',4);
y1=y1(1:sig_num);
y2=dmod(x2,1,0.1,4,'psk',4);
y2=y2(1:sig_num);
x=[y1';y2'];
s=x;
q=2;
p=12;%智能天线阵元数
fc=10.7*10^6;
lemda=3*10^8/fc;
banjing=20;
jh=pi./180;
k0=2*pi/lemda;
h=ceil(k0*banjing);
P=2*h+1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%均匀圆阵转化为均匀线阵
alpha=[90,90];
theta=[20,50];
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));
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
cv1(i+h+1)=j^(-abs(i));
end
cv=diag(cv1);
for m=-h:h
for i=0:p-1
      w(i+1,m+h+1)=exp(-j*m*gama(i+1));
end
end
V=1./sqrt(p)*w;
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;
rx=xn*xn'./sig_num;
ry=Fr*rx*Fr';
r=real(ry);
% r=Fe*rx*Fe';

[u,sig,v]=svd(r);%%svd奇异值分解，特征值降序排列
  signalspace=u(:,q);
noisespace=u(:,q+1:7);
% vv=conj(noisespace)*noisespace';
UU=noisespace*noisespace.';

%%%%%%root-music%%%%%%%%%%%%%%%%%

J1=besselj(1:h,k0*banjing*sin(90*jh));
J11=diag([fliplr(J1),besselj(0,k0*banjing*sin(90*jh)),J1]);

Q=J11*W*UU*W'*J11;

num1=1;
num2=1;
A1(1,1: P-1)=0;
A2(1,1: P-1)=0;
for l=P-1:-1:1
j=P;
jj=j-l;
for n=1:jj
      i=j-l;
      A1(1,num1)=A1(1,num1)+Q(i,j);
      A2(1,num2)=A2(1,num2)+Q(j,i);
      j=j-1;
      j;
end
num1=num1+1;
num2=num2+1;
end
A1;
A2;
A=[A1, trace(Q),fliplr(A2)];
ra=roots([A]);
rb=ra(abs(ra)<1);
[dumm,I]=sort(abs(abs(rb)-1));
w=angle(rb(I(1:q)));
fi=-w/jh

[ Last edited by 影月水澜 on 2013-7-2 at 21:53 ]

回复此楼