matlab 循环程序的问题

matlab 循环程序的问题
这段程序的问题是,只要加入C33矩阵的部分,谱估计结果就不对,就算特征分解不用C33也不行,可是从数值上看C33和C3是一样的,问题到底出在哪里?
%经典MUSIC算法
%信号不相干
clear all;
clc;
%初始数据
sensor_number=5; %阵元数
N_x=1024; %快拍数
f=8*10^9; %信号频率
l=3*10^8/f; %波长
d=0.5*l; %阵元间距
snr=20;
snr1=10.^(snr/20); %信噪比
theta=[0 20 40]; %信号角度
rtheta=theta*pi/180;
ll=length(theta);
%生成信号
for m=1:sensor_number;
for n=1:ll;
A(m,n)=exp(-j*2*pi*(m-1)*d*sin(rtheta(n))/l);
end
end
%解析信号虚部是实部的希尔伯特变换,采用复信号的表示方法,模拟入射信号
for m=1:ll;
for n=1:N_x;
p1=rand(1,1);
p2=rand(1,1);
sr(m,n)=sqrt(-2*snr1*snr1*log(p1))*cos(2*pi*p2);
si(m,n)=sqrt(-2*snr1*snr1*log(p1))*sin(2*pi*p2);
s(m,n)=sr(m,n)+j*si(m,n);
end
end
%高斯白噪声
for m=1:sensor_number;
for n=1:N_x;
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
x=A*s+w;
%构造协方差矩阵,特征分解
C2=x*x'/N_x;
C33=zeros(M,M);
%for j=1:M
%for i=1:M-j+1
%C33(i,i+j-1)=C2(1,j);
% end
%end
%for j=2:M
%for i=1:M-j+1
%C33(i+j-1,i)=conj(C2(1,j));
%end
%end
C3=zeros(M,M);
for i=1:M
C3(i,i)=C2(1,1);
end
for i=1:M-1
C3(i,i+1)=C2(1,2);
end
for i=1:M-1
C3(i+1,i)=conj(C2(1,2));
end
for i=1:M-2
C3(i,i+2)=C2(1,3);
end
for i=1:M-2
C3(i+2,i)=conj(C2(1,3));
end
for i=1:M-3
C3(i,i+3)=C2(1,4);
end
for i=1:M-3
C3(i+3,i)=conj(C2(1,4));
end
for i=1:M-4
C3(i,i+4)=C2(1,5);
end
for i=1:M-4
C3(i+4,i)=conj(C2(1,5));
end
[V D]=eig(C3);
[lambda index]=sort(diag(D),'descend');
lambda=lambda';
%AIC准则判断信源数
for k=0:sensor_number-1
temp0=sum(lambda(k+1:sensor_number))/(sensor_number-k);
temp1=prod(lambda(k+1:sensor_number))^(1/(sensor_number-k));
Lambda=temp0/temp1;
AIC(k+1)=2*N_x*(sensor_number-k)*log(Lambda)+2*k*(2*sensor_number-k);
end
[aic_min b]=min(AIC);
num=b-1;%估计出的信源数
disp('source_number=');
disp(num);
%噪声子空间
Un=V(:,index(num+1:sensor_number));
Gn=Un*Un';
%谱峰搜索
searching_doa=-90:0.1:90;
for i=1:length(searching_doa)
a_theta=exp(-j*(0:sensor_number-1)'*2*pi*d*sin(pi*searching_doa(i)/180)/l);
Pmusic(i)=1./abs((a_theta)'*Gn*a_theta);
end
plot(searching_doa,10*log10(Pmusic),'r');
xlabel('入射角/degree');
ylabel('空间谱/dB');
legend('Music Spectrum');
title('经典MUSIC估计');
grid on;

改了点东西,你看你是不是这个意思?
clear all;
clc;
%初始数据
sensor_number=5; %阵元数
N_x=1024; %快拍数
f=8*10^9; %信号频率
l=3*10^8/f; %波长
d=0.5*l; %阵元间距
snr=20;
snr1=10.^(snr/20); %信噪比
theta=[0 20 40]; %信号角度
rtheta=theta*pi/180;
ll=length(theta);
%生成信号
for m=1:sensor_number;
for n=1:ll;
A(m,n)=exp(-1i*2*pi*(m-1)*d*sin(rtheta(n))/l);
end
end
%解析信号虚部是实部的希尔伯特变换,采用复信号的表示方法,模拟入射信号
for m=1:ll;
for n=1:N_x;
p1=rand(1,1);
p2=rand(1,1);
sr(m,n)=sqrt(-2*snr1*snr1*log(p1))*cos(2*pi*p2);
si(m,n)=sqrt(-2*snr1*snr1*log(p1))*sin(2*pi*p2);
s(m,n)=sr(m,n)+1i*si(m,n);
end
end
%高斯白噪声
for m=1:sensor_number;
for n=1:N_x;
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)+1i*wi(m,n);
end
end
x=A*s+w;
%构造协方差矩阵,特征分解
C2=x*x'/N_x;
M=m;
C33=zeros(M,M);
%for j=1:M
%for i=1:M-j+1
%C33(i,i+j-1)=C2(1,j);
% end
%end
%for j=2:M
%for i=1:M-j+1
%C33(i+j-1,i)=conj(C2(1,j));
%end
%end
C3=zeros(M,M);
for i=1:M
C3(i,i)=C2(1,1);
end
for i=1:M-1
C3(i,i+1)=C2(1,2);
end
for i=1:M-1
C3(i+1,i)=conj(C2(1,2));
end
for i=1:M-2
C3(i,i+2)=C2(1,3);
end
for i=1:M-2
C3(i+2,i)=conj(C2(1,3));
end
for i=1:M-3
C3(i,i+3)=C2(1,4);
end
for i=1:M-3
C3(i+3,i)=conj(C2(1,4));
end
for i=1:M-4
C3(i,i+4)=C2(1,5);
end
for i=1:M-4
C3(i+4,i)=conj(C2(1,5));
end
[V D]=eig(C3);
[lambda index]=sort(diag(D),'descend');
lambda=lambda';
%AIC准则判断信源数
for k=0:sensor_number-1
temp0=sum(lambda(k+1:sensor_number))/(sensor_number-k);
temp1=prod(lambda(k+1:sensor_number))^(1/(sensor_number-k));
Lambda=temp0/temp1;
AIC(k+1)=2*N_x*(sensor_number-k)*log(Lambda)+2*k*(2*sensor_number-k);
end
[aic_min b]=min(AIC);
num=b-1;%估计出的信源数
disp('source_number=');
disp(num);
%噪声子空间
Un=V(:,index(num+1:sensor_number));
Gn=Un*Un';
%谱峰搜索
searching_doa=-90:0.1:90;
for i=1:length(searching_doa)
a_theta=exp(-1i*(0:sensor_number-1)'*2*pi*d*sin(pi*searching_doa(i)/180)/l);
Pmusic(i)=1./abs((a_theta)'*Gn*a_theta);
end
plot(searching_doa,10*log10(Pmusic),'r');
xlabel('入射角/degree');
ylabel('空间谱/dB');
legend('Music Spectrum');
title('经典MUSIC估计');
grid on;