www.pudn.com > MUSCI_COPRIME.zip > MUSCI_COPRIME.m, change:2015-06-16,size:3481b


clc; 
clear all; 
close all; 
%% initialization 
fc = 1e9; %信号频率 1GHz 
fs = 2e9; %采样频率 2GHz 
lambda = 3e8/fc; %信号波长 0.3m 
d = lambda/2; %天线阵列间最小距离,即半波长,单位m 
snr = 0; %各信源在接收端统一的信噪比 
sampling_points = 500;  %快照数 
 
N=7;M=5; 
numbers_M = 20; %目标个数 
precision = 0.01; %精度,0.01度 
%% position 
%%%%%%%%%%%% 第一阵列%%%%%%%%%% 
% N=5 No_points = N = 5 
numbers_N_array_1 = N; %阵列中天线个数 
d_array_1 = M*d; %第一阵列天线间间距 
X_receivers_1 = (0:1:numbers_N_array_1-1).*d_array_1; %阵列中各天线的x轴坐标 
Y_receivers_1 = zeros(1,numbers_N_array_1); %阵列中各天线的y轴坐标 
%%%%%%%%%%%% 第二阵列%%%%%%%%%% 
% M=3 No_points = 2M-1 = 5 
numbers_N_array_2 = 2*M-1; %阵列中天线个数 
d_array_2 = N*d; %第一阵列天线间间距 
X_receivers_2 = (1:1:numbers_N_array_2).*d_array_2; %阵列中各天线的x轴坐标 
Y_receivers_2 = zeros(1,numbers_N_array_2); %阵列中各天线的y轴坐标 
%%%%%%%%%%%% 阵列重构%%%%%%%%%% 
X_array_co = [X_receivers_1 X_receivers_2]; 
X_array_co = sort(X_array_co); 
%%%%%%%%%%%% 信源%%%%%%%%%% 
% numbers_M = 10; %目标个数 
distance_SendPoint = 2.*(numbers_N_array_1.^2).*2.*d; %远场条件 
 
% theta = ( 2.*rand(1,numbers_M)-1 ).*90; %随机生成的角度 
gap = round (120/numbers_M); 
i = 1:1:numbers_M; 
theta = -60+i*gap; 
 
X_SendPoint = distance_SendPoint.*sin(theta./180.*pi); %随机生成的目标的x轴坐标 
Y_SendPoint = distance_SendPoint.*cos(theta./180.*pi); %随机生成的目标的y轴坐标 
i = 1:1:numbers_M; 
if( numbers_M/2 == 0) 
    fc_1(i) = fc+1e7*i-1e7*numbers_M/2+5e6; %各信源单频信号的频率点确定 
else 
    fc_1(i)=fc+1e7*i-1e7*(numbers_M+1)/2; 
end 
 
%figure % 散点图 
figure(1) 
scatter(X_receivers_1,Y_receivers_1,'k'); 
hold on; 
scatter(X_receivers_2,Y_receivers_2,'k','fill'); 
scatter(X_SendPoint,Y_SendPoint,'d'); 
 
%% sampling 
t = 0:1/fs:(sampling_points-1)/fs; % 500个采样点 
for i = 1:numbers_M 
    s(i,:) = exp(j*2*pi*fc_1(i).*t); %接收端的信号 
end 
for i = 1:numbers_M 
    phi = 2*pi*sin(theta(i)./180.*pi)/lambda;%固定相位差 
    % a = exp(-j.*2*pi*distance_SendPoint*sin(theta./180.*pi)./(2*d).*difference); %相位差矩阵 
    a_array(i,:) = exp(-j.*phi.*X_array_co);%相位差矩阵 
end 
a_array = a_array.';  %求转置.'  共轭转置' 
x = a_array*s; 
x = awgn(x,snr,'measured'); 
 
%%%%%%%%%%%%%%%%%%%test 
if 0 
% x = a_array; 
% sampling_points = 1; 
i = 0:1:N*M; 
array_continuous = i*d; 
for i = 1:numbers_M 
    phi = 2*pi*sin(theta(i)./180.*pi)/lambda;%固定相位差 
    a_array_continuous(i,:) = exp(-j.*phi.*array_continuous); 
end 
a_array_continuous = a_array_continuous.'; 
y = a_array_continuous*s; 
y1 = sum(y,2)/sampling_points; 
R_aa = y*y'/sampling_points; 
end 
%%%%%%%%%%%%%%%%%%%%% 
 
%% COPRIME_MUSIC 
R_xx = x*x'/sampling_points;  %x的自相关函数 
vec_R_xx = reconsitution(R_xx,X_array_co,d,N*M); 
R_xx = Spatially_smoothed(vec_R_xx);  %空间平滑,2MN+1 to MN+1 
 
%% MUSIC 
[V,D] = eig(R_xx); % 特征值对角矩阵D,特征向量V,满足R_xx*V=V*D 
[len,len] = size(R_xx); 
for i = 1:1:len; 
    d_eig(i) = D(i,i); 
end 
U_noise = V(:,1:1:(len-numbers_M)); 
difference = 0:1:N*M; %相位差矩阵差异 
for i = 1:1:ceil(180/precision+1) 
    theta0 = (i-1)*precision-90; 
    find_phi = pi*sin(theta0./180.*pi); 
    find_a = exp(-j*find_phi.*difference).'; 
    P_MUSIC(i) = 1/(find_a'*U_noise*U_noise'*find_a); 
    P_MUSIC(i) = abs(P_MUSIC(i)); 
end 
figure(2) 
P_MUSIC = P_MUSIC/max(P_MUSIC); 
plot(-90:precision:90,P_MUSIC); 
hold on; 
for i = 1:numbers_M 
    plot(theta(i),0:0.005:max(P_MUSIC)); 
end