自適應(yīng)波束形成與Matlab程序代碼注解

1、1. 均勻線陣方向圖（1）matlab 程序clc;clear all;close all;imag=sqrt(-1);element_num=32;%陣元數(shù)為8d_lamda=1/2;%陣元間距d與波長(zhǎng)lamda的關(guān)系theta=linspace(-pi/2,pi/2,200);theta0=0;%來(lái)波方向w=exp(imag*2*pi*d_lamda*sin(theta0)*0:element_num-1');for j=1:length(theta) a=exp(imag*2*pi*d_lamda*sin(theta(j)*0:element_num-1'); p(j)=

2、w'*a;endpatternmag=abs(p);patternmagnorm=patternmag/max(max(patternmag);patterndB=20*log10(patternmag);patterndBnorm=20*log10(patternmagnorm);figure(1)plot(theta*180/pi,patternmag);grid on;xlabel('theta/radian')ylabel('amplitude/dB')title(num2str(element_num) '陣元均勻線陣方向圖',

3、'來(lái)波方向?yàn)?#39; num2str(theta0*180/pi) '度');hold on;figure(2)plot(theta,patterndBnorm,'r');grid on;xlabel('theta/radian')ylabel('amplitude/dB')title(num2str(element_num) '陣元均勻線陣方向圖','來(lái)波方向?yàn)?#39; num2str(theta0*180/pi) '度');axis(-1.5 1.5 -50 0); (2)

4、仿真結(jié)果A.來(lái)波方向?yàn)?° 不歸一化歸一化B.來(lái)波方向?yàn)?5°不歸一化歸一化C.隨著陣元數(shù)的增加，波束寬度變窄，分辨力提高，仿真圖如下： 非歸一化歸一化 不歸一化 歸一化2. 波束寬度與波達(dá)方向及陣元數(shù)的關(guān)系（1）matlab 程序clc;clear all;close all;imag=sqrt(-1);element_num1=16;element_num2=128;element_num3=1024;lambda=0.1;d=0.5*lambda;theta=0:0.5:90;for j=1:length(theta)fai(j)=theta(j)*pi/180-as

5、in(sin(theta(j)*pi/180)-lambda/(element_num1*d);psi(j)=theta(j)*pi/180-asin(sin(theta(j)*pi/180)-lambda/(element_num2*d);beta(j)=theta(j)*pi/180-asin(sin(theta(j)*pi/180)-lambda/(element_num3*d);endfigureplot(theta,fai,'r',theta,psi,'b',theta,beta,'g');grid on;xlabel('the

6、ta');ylabel('width in radians');title('波束寬度與達(dá)波方向及陣元數(shù)目的關(guān)系');legend('N=16','N=128','N=1024');(2) 仿真結(jié)果 結(jié)果3. 當(dāng)陣元間距時(shí)，會(huì)出現(xiàn)柵瓣，導(dǎo)致空間模糊（1）仿真結(jié)果 非歸一化歸一化4. 類(lèi)似于時(shí)域?yàn)V波，天線方向圖是最優(yōu)權(quán)的傅立葉變換（1）matlab 程序clc;clear all;close all;imag=sqrt(-1);element_num=32;source_num=1;d_lambda=0.5;

7、theta=linspace(-pi/2,pi/2,200);theta0=0;w=exp(imag*2*pi*d_lambda*sin(theta0)*0:element_num-1');for j=1:length(theta) a=exp(imag*2*pi*d_lambda*sin(theta(j)*0:element_num-1'); p(j)=w'*a;endpatternmag=abs(p);patternmagnorm=patternmag/max(max(patternmag);patterndB=20*log10(patternmag);patter

8、ndBnorm=20*log10(patternmagnorm);figure(1)subplot(1,2,1);plot(theta,patterndBnorm);grid on;xlabel('theta/radian');ylabel('amplitude/dB');axis(-2.0 2.0 -50 0);subplot(1,2,2);pfft=fftshift(fft(w,256);pfftmag=abs(pfft);pfftmagnorm=pfftmag/max(max(pfftmag);pfftdB=20*log10(pfftmagnorm);pf

9、ftdBnorm=20*log10(pfftmagnorm);plot(linspace(-pi/2,pi/2,256),pfftdBnorm);grid on;xlabel('theta/radian');ylabel('FFT_amplitude/dB');axis(-2.0 2.0 -50 0);（2）仿真結(jié)果5. 最大信噪比準(zhǔn)則方向圖和功率譜（1） matlab 程序clc;clear all;close all;imag=sqrt(-1);element_num=8;%陣元數(shù)為8d_lambda=0.5;%間距為半波長(zhǎng)theta=-90:0.5:90;

10、%掃描范圍theta0=0;%來(lái)波方位theta1=20;%干擾方向L=512;%采樣點(diǎn)數(shù)for i=1:L amp0=10*randn(1); amp1=200*randn(1); ampn=1; s(:,i)=amp0*exp(imag*2*pi*0.5*sin(theta0*pi/180)*0:element_num-1'); j(:,i)=amp1*exp(imag*2*pi*0.5*sin(theta1*pi/180)*0:element_num-1'); n(:,i)=ampn*exp(randn(element_num,1)+imag*randn(element_

11、num,1);endRs=1/L*s*s'%信號(hào)自相關(guān)矩陣Rnj=1/L*(j*j'+n*n'); %干擾+噪聲的自相關(guān)矩陣V,D=eig(Rs,Rnj); %（Rs,Rnj）的廣義特征值和特征向量D,I=sort(diag(D); %特征向量排序Wopt=V(:,I(8);%最優(yōu)權(quán)矢量for j=1:length(theta)a=exp(imag*2*pi*d_lambda*sin(theta(j)*pi/180)*0:element_num-1'); f(j)=Wopt'*a; p(j)=a'*Rs*a+a'*Rnj*a;endF=2

12、0*log10(abs(f)/max(max(abs(f);P=20*log10(abs(p)/max(max(abs(p);subplot(1,2,1)plot(theta,F);grid on;hold on;plot(theta0,-80:0,'.');plot(theta1,-80:0,'.');xlabel('theta/0');ylabel('F in dB');title('max-SNR 方向圖');axis(-90 90 -80 0);hold on;subplot(1,2,2);plot(the

13、ta,P,'r');grid on;xlabel('theta/0');ylabel('功率 in dB');title('max-SNR 功率譜');grid on;axis(-90 90 -80 0);（2）仿真結(jié)果6. ASC旁瓣相消-MSE準(zhǔn)則(1) matlab 程序clc;close all;clear all;imag=sqrt(-1);M=32;%輔助天線數(shù)目d_lambda=0.5;%陣元間距theta0=-30;%來(lái)波方向theta1=60;%干擾方向L=512;%采樣單元數(shù)s=zeros(1,512); %預(yù)

14、劃分一個(gè)區(qū)域for ii=1:L amp0=1*randn(1);%信號(hào)的幅度隨機(jī)產(chǎn)生，保證信號(hào)之間是不相關(guān)的 amp1=200*randn(1);ampn=1; jam(:,ii)=amp1*exp(imag*2*pi*d_lambda*sin(theta1*pi/180)*0:M-1')+ampn*(randn(M,1)+imag*randn(M,1); %干擾+噪聲 s(ii)=amp0*exp(imag*2*pi*d_lambda*sin(theta0*pi/180)+amp1*exp(imag*2*pi*d_lambda*sin(theta1*pi/180)+ampn*(ra

15、ndn(1,1)+imag*randn(1,1);%接收信號(hào)（信號(hào)+干擾+噪聲）s0(ii)=amp0*exp(imag*2*pi*d_lambda*sin(theta0*pi/180);endRx=1/L*jam*jam'r_xd=1/L*jam*s'Wopt=pinv(Rx)*r_xd;delta=s0-(s-Wopt'*jam);delta1=abs(mean(delta.2)-(mean(delta).2);theta=linspace(-pi/2,pi/2,200);for jj=1:length(theta) a=exp(imag*2*pi*d_lambda

16、*sin(theta(jj)*0:M-1'); f(jj)=Wopt'*a;endF=20*log10(abs(f)/max(max(abs(f);figure(1)plot(theta*180/pi,F);grid on;hold on;plot(theta0,-50:0,'.');plot(theta1,-50:0,'.');xlabel('theta/°');ylabel('F/dB');title('MSE準(zhǔn)則下的方向圖');axis(-90 90 -50 0);（2）仿真結(jié)果7.

17、 線性約束最小方差(LCMV)準(zhǔn)則 (1) matlab 程序clc;clear all;close all;imag=sqrt(-1);element_num=8;%陣元數(shù)d_lambda=0.5;%陣元間距與波長(zhǎng)的關(guān)系theta=-90:0.5:90; %搜索范圍theta0=0; %三個(gè)信號(hào)源的來(lái)波方向theta1=30;theta2=60;L=512;%采樣單元數(shù)for i=1:L amp0=10*randn(1); amp1=100*randn(1); amp2=10*randn(1); ampn=10;x(:,i)=amp0*exp(imag*2*pi*d_lambda*sin(t

18、heta0*pi/180)*0:element_num-1')+amp1*exp(imag*2*pi*d_lambda*sin(theta1*pi/180)*0:element_num-1')+amp2*exp(imag*2*pi*d_lambda*sin(theta2*pi/180)*0:element_num-1')+ampn*(randn(element_num,1)+imag*randn(element_num,1);endRx=1/L*x*x'steer1=exp(imag*2*pi*d_lambda*sin(theta0*pi/180)*0:elem

19、ent_num-1');steer2=exp(imag*2*pi*d_lambda*sin(theta1*pi/180)*0:element_num-1');steer3=exp(imag*2*pi*d_lambda*sin(theta2*pi/180)*0:element_num-1');C=steer1 steer2 steer3;F=1 0 1'%把三個(gè)方向都作為來(lái)波方向w=inv(Rx)*C*(inv(C'*inv(Rx)*C)*F;for j=1:length(theta) a=exp(imag*2*pi*d_lambda*sin(theta(

20、j)*pi/180)*0:element_num-1'); f(j)=w'*a; p(j)=1/(a'*inv(Rx)*a);endF=20*log10(abs(f)/(max(max(abs(f);subplot(1,2,1)plot(theta,F);grid on;hold on;plot(theta0,-20:0,'.');plot(theta1,-20:0,'.');plot(theta2,-20:0,'.');xlabel('theta/°');ylabel('F/dB'

21、;);title('Capon beamforming 方向圖');axis(-90 90 -20 0);P=20*log10(abs(p)/(max(max(abs(p);subplot(1,2,2)plot(theta,P);grid on;hold on;plot(theta0,-20:0,'.');plot(theta1,-20:0,'.');plot(theta2,-20:0,'.');xlabel('theta/°');ylabel('P/dB');title('Cap

22、on beamforming 功率譜');axis(-90 90 -20 0);（2） 仿真結(jié)果8. Capon beamforming（1） matlab 程序clc;clear all;close all;imag=sqrt(-1);element_num=8;%陣元數(shù)d_lambda=0.5;%陣元間距與波長(zhǎng)的關(guān)系theta=-90:0.5:90; %搜索范圍theta0=0; %三個(gè)信號(hào)源的來(lái)波方向theta1=20;theta2=60;L=1000;%采樣單元數(shù)for i=1:L amp0=10*randn(1); amp1=200*randn(1); amp2=200*ra

23、ndn(1); ampn=3; x(:,i)=amp0*exp(imag*2*pi*d_lambda*sin(theta0*pi/180)*0:element_num-1')+amp1*exp(imag*2*pi*d_lambda*sin(theta1*pi/180)*0:element_num-1')+amp2*exp(imag*2*pi*d_lambda*sin(theta2*pi/180)*0:element_num-1')+ampn*(randn(element_num,1)+imag*randn(element_num,1);endRx=1/L*x*x'

24、;R=inv(Rx);steer=exp(imag*2*pi*d_lambda*sin(theta0*pi/180)*0:element_num-1');w=R*steer/(steer'*R*steer);%最優(yōu)權(quán)矢量for j=1:length(theta) a=exp(imag*2*pi*d_lambda*sin(theta(j)*pi/180)*0:element_num-1'); f(j)=w'*a; p(j)=1/(a'*R*a);endF=20*log10(abs(f)/(max(max(abs(f);subplot(1,2,1)plot(

25、theta,F);grid on;hold on;plot(theta0,-50:0,'.');plot(theta1,-50:0,'.');plot(theta2,-50:0,'.');xlabel('theta/°');ylabel('F/dB');title('Capon beamforming 方向圖');axis(-90 90 -50 0);P=20*log10(abs(p)/(max(max(abs(p);subplot(1,2,2)plot(theta,P);grid on;

26、hold on;xlabel('theta/°');ylabel('P/dB');title('Capon beamforming 功率譜');axis(-90 90 -90 0);（2） 仿真結(jié)果9. 不同方法估計(jì)協(xié)方差矩陣的Capon波束形成（1）matlab 程序clc;clear all;close all;imag=sqrt(-1);element_num=8;%陣元數(shù)為8d_lambda=0.5;%間距為半波長(zhǎng)theta=-90:0.5:90;%掃描范圍theta0=0;%來(lái)波方向theta1=50;%干擾方向L=1024;

27、%采樣單元數(shù)for i=1:L amp0=10*randn(1); amp1=50*randn(1); ampn=0.5; s(:,i)=amp0*exp(imag*2*pi*d_lambda*sin(theta0*pi/180)*0:element_num-1'); j(:,i)=amp1*exp(imag*2*pi*d_lambda*sin(theta1*pi/180)*0:element_num-1'); n(:,i)=ampn*exp(imag*2*pi*randn(1)*0:element_num-1');endRx=1/L*(s+j+n)*(s+j+n)&#

28、39;%接收信號(hào)自相關(guān)矩陣Rnj=1/L*(j+n)*(j+n)'%干攏+噪聲的自相關(guān)矩陣e=exp(imag*2*pi*d_lambda*sin(theta0*pi/180)*0:element_num-1');Wopt_Rx=inv(Rx)*e/(e'*inv(Rx)*e);%采用接收信號(hào)的權(quán)矢量Wopt_Rnj=inv(Rnj)*e/(e'*inv(Rnj)*e);%采用干攏+噪聲信號(hào)的權(quán)矢量for j=1:length(theta) a=exp(imag*2*pi*d_lambda*sin(theta(j)*pi/180)*0:element_num-1

29、'); f1(j)=Wopt_Rx'*a; f2(j)=Wopt_Rnj'*a;endF1=20*log10(abs(f1)/max(max(abs(f1);F2=20*log10(abs(f2/max(max(abs(f2);figure;plot(theta,F1,theta,F2,'r');grid on;hold on;plot(theta0,-50:0,'.');plot(theta1,-50:0,'.');xlabel('theta/°');ylabel('F(1,2)/dB&

30、#39;);title('不同方法估計(jì)協(xié)方差矩陣的Capon波束形成');axis(-90 90 -60 0);（2）仿真結(jié)果 10 多點(diǎn)約束的Capon波束形成和方向圖（1）matlab 程序clc;clear all;close all;imag=sqrt(-1);element_num=8;d_lambda=0.5;theta=-90:0.3:90;theta0=0;theta1=20;theta2=50;L=512;Rx=zeros(element_num,element_num);%產(chǎn)生協(xié)方差矩陣for i=1:L amp0=10*randn(1); amp1=10*

31、randn(1); amp2=50*randn(1);ampn=0.5*randn(1);%噪聲的幅度隨機(jī)產(chǎn)生，保證噪聲與信號(hào)之間是不相關(guān)的 j(:,i)=amp1*exp(imag*2*pi*d_lambda*sin(theta1*pi/180)*0:element_num-1')+amp2*exp(imag*2*pi*d_lambda*sin(theta2*pi/180)*0:element_num-1')+ampn*exp(imag*2*pi*randn(1)*0:element_num-1'); x(:,i)=amp0*exp(imag*2*pi*d_lambd

32、a*sin(theta0*pi/180)*0:element_num-1')+j(:,i);%表示接收信號(hào)endRx=1/L*x*x'R=inv(Rx);w=amp0*exp(imag*2*pi*d_lambda*sin(theta0*pi/180)*0:element_num-1')+amp1*exp(imag*2*pi*d_lambda*sin(theta1*pi/180)*0:element_num-1')+amp2*exp(imag*2*pi*d_lambda*sin(theta2*pi/180)*0:element_num-1');for j=

33、1:length(theta) a=exp(imag*2*pi*d_lambda*sin(theta(j)*pi/180)*0:element_num-1'); f(j)=w'*a; p(j)=1/(a'*R*a);endF=20*log10(abs(f)/max(max(abs(f);P=20*log10(abs(p)/max(max(abs(p);figure;subplot(1,2,1);plot(theta,F);grid on;hold on;plot(theta0,-50:0,'.');plot(theta1,-50:0,'.'

34、;);plot(theta2,-50:0,'.');xlabel('theta/°');ylabel('F/dB');title('Capon beamforming方向圖');axis(-90 90 -50 0);subplot(1,2,2);plot(theta,P);hold on;grid on;plot(theta0,-90:0,'.');plot(theta1,-90:0,'.');plot(theta2,-90:0,'.');xlabel('theta

35、/°');ylabel('P/dB');title('Capon beamforming功率譜');（2） 仿真結(jié)果11 自適應(yīng)波束形成方向圖（1）matlab 程序clc;clear all;close all;imag=sqrt(-1);element_num=8;c=3e8;f=5e8;lambda=c/f;d_lambda=0.5;theta=-90:0.5:90;theta0=0;theta1=45;theta2=60;L=2048;for i=1:L amp0=10*randn(1); amp1=100*randn(1); amp2=100*randn(1); ampn=1