线性约束最小方差准则

clear all; 
close all; 
clc; 
%--------------------------------- 
%信号参数设置 
Signal_No = 1; 
Interference_No = 2;%小于五 
S_No = Signal_No + Interference_No; %信号总数
%Xm(n)=S(n)+Nm(n)
sensor_No = 12; %阵元数 
%阵元M
azimuth = [0*pi/180 -40*pi/180 50*pi/180 30*pi/180 20*pi/180]; %方位角
%Ai=[A1.i，A2.i，....Am.i] 
Fs = 50e6; 
Fc = 10e6; 
F1 = [8e6 9e6 10e6 11e6]; %干扰功率
RB = 1e6; 
M = 2; %二进制 
wavelength = 3e8/Fc; 
d = wavelength/2; 
K = 300;  %快拍数 
Ps = 0; %dB 
SNR = -15; 
SIR = -90; 
M_No = 40; %码速率等于信息速率 
Data_No = M_No*Fs/RB; 
t = 1/Fs:1/Fs:Data_No/Fs; 
%--------------------------------- 
%功率转换 
Ps_l = 10^(Ps/10); 
Pi_l = 10^((Ps-SIR)/10); 
Pn_l = 10^((Ps-SNR)/10); 
%---------------------------------
 %信号生成 
bit = randint(1,M_No);          %产生信息序列 
bitstream = []; 
for i = 1:M_No 
 if bit(i) == 1 
      bitstream = [bitstream, ones(1, Fs/RB)];
   else 
      bitstream = [bitstream, -ones(1, Fs/RB)];  
  end 
end 
Carrier_R = cos(2*pi*Fc*t); 
S_R = Carrier_R .* bitstream; 
Carrier_I = sin(2*pi*Fc*t); 
S_I = Carrier_I .* bitstream; 
Signal_R = sqrt(Ps_l)*S_R/sqrt((S_R*S_R'/length(S_R))); 
Signal_I = sqrt(Ps_l)*S_I/sqrt((S_I*S_I'/length(S_I)));%生成BPSK复信号 
Signal(1,:) = complex(Signal_R, Signal_I); 
%--------------------------------- 
% %产生干扰 
for i = 2:Interference_No +1 
  Signal(i,:)=sqrt(2*Pi_l)*complex(sin(2*pi*F1(i-1)*t+pi/10*i).*sin(2*pi*F1(i-1)*t+pi/10*i));
end 
%--------------------------------- 
% 模拟天线接收 
%---------------------------------- 
%四次循环 
name = [' 3',' 6',' 9','12']; 
for mmm =1:4 
sensor_No = mmm*3; 
As = zeros(sensor_No,S_No); 
for i = 1:sensor_No
   for ii = 1:S_No 
    As(i,ii) = exp(-1i*2*pi*sin(azimuth(ii))*d/wavelength*(i-1)); 
   end 
end 
 x = As * Signal;
% 
% %--------------------------------- 
% %加噪声 
noise_R = randn(sensor_No,Data_No)*sqrt(Pn_l); 
noise_I = randn(sensor_No,Data_No)*sqrt(Pn_l); 
noise = complex(noise_R, noise_I); s = x + noise; 
%--------------------------------- 
%lcmv算法 
% Ad = [1 0 0 0 0 0 0 0 0 0 0 0 ]'; 
Ad = (As(:,1)); 
x1 = s(:,(1:K)); 
Rx=x1*x1'/K; 
 alfa=inv(Ad'*inv(Rx)*Ad); 
 Wopt=alfa*inv(Rx)*Ad; 
%------------------------ 
%---------------------------------
%mvdr 
% Ad = (As(:,1));
% x1 = s(:,(1:K)); 
% Rx=x1*x1'/K; 
% QR_Wopt1 = inv(Rx)*Ad; 
% Wopt=QR_Wopt1./QR_Wopt1(1); 
%--------------------------------- 
%波束图 
figure(1) 
whitebg('k'); 
all_azimuth = -pi/2:pi/180:pi/2; 
LMS_R = zeros(1,181); 
for i = 1:181 
  for m = 1:sensor_No
     LMS_R(i) = LMS_R(i) + conj(Wopt(m)) * exp(-1i*2*pi*sin(all_azimuth(i))*d/wavelength*(m-1)); 
  end 
  LMS_P(i) = (abs(LMS_R(i)))^2; 
end 
LMS_PdB = 10*log10(LMS_P/max(LMS_P));
 switch mmm 
   case 1 
    plot(all_azimuth*180/pi,LMS_PdB,'--g','LineWidth',2);hold on;
   case 2 
    plot(all_azimuth*180/pi,LMS_PdB,'.b','LineWidth',2);
   case 3 
    plot(all_azimuth*180/pi,LMS_PdB,'-.w','LineWidth',2);
   case 4 
    plot(all_azimuth*180/pi,LMS_PdB,':r','LineWidth',2);
 end 
grid on; 
axis([-90 90 -70 10]);
xlabel('方位角(度)') 
ylabel('阵列增益(dB)') 
%title(strcat('阵元数',name([2*mmm-1:2*mmm]),'个  ')); 
% figure(2); 
% polar ((-pi/2:pi/180:pi/2),LMS_PdB+109); figure(gcf);
 end 
hl=legend('阵元数3个','阵元数6个','阵元数9个','阵元数12个',1); 
%set(hl,'edgecolor','w','color',[1 1 1]); 
set(hl,'edgecolor','w','color',[0 0 0]); 
title('自适应波束形成'); 
drawnow;