MVDR波束形成Matlab代码

% Multiple beamforming. clear;clc;close all; starttime=now; %%%%%%%%%%%%%%%%%%The initial condition%%%%%%%%%%%%%%%%%%%% N=8; % the array elments no. d=0.5; % in wavelength unit. Var=1; % white noise variance. p=1; % no. of desired signal SNR=[10 40]'; % source SNR in dB. Angs1=[0.2 -0.4]'; % the incident angle vector in radians. L=length(Angs1); % sources no. snapshot=200; % Sampling data number. sir_data=200; % Iteration numbers of LMS algorithm experimentno=20; % Averaging experimentno. g=ones(p,1); ave_sir=zeros(1,sir_data); sir_end=0; step=10^(-7); % Stepsize %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% j=sqrt(-1); n=[0:1:N-1]'; pj=zeros(L,1); for i=1:L pj(i)=Var*10^(SNR(i)/10); end amp=sqrt(pj); %%%%%%%% Far-Field %%%%%%%% angf=2*pi*d*Angs1; aaf=n*angf'; D=exp(j*aaf); % the Direction Matrix. %%%%%%%%%%%%%%%% generate the autocorrelation matrix %%%%% ny=zeros(N,1); for ite=1:experimentno ite x=zeros(N,snapshot); %%%%% Generate the BPSK source %%%%% % sam_t=ones(L,1); % the normalized sampling interval. % pha=zeros(L,1); % sources initial phase. % chg_pha=ones(L,1); % Data length to exchage BPSK. % factor=ones(L,1); % for k=1:snapshot % for kk=1:L % if rem(k,chg_pha(kk)) == 0 % b_level=round(rand); % factor(kk)=round(b_level-0.5)*factor(kk); % end % end % bpha=pi/2*factor; % source(:,k)=amp.*exp(j*(2*pi*(k-1)*sam_t+bpha)); % end %%%%%%% Genrate project's source %%%%%% source=kron(amp,ones(1,snapshot)); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% x=D*source; samn=randn(N,snapshot)+j*randn(N,snapshot); samn=samn/sqrt(2)*Var; y=x+samn; Rxx=y*y'/snapshot; %%%%%%%%% MVDR %%%%%%%%% w=inv(Rxx)*D(:,1:p)*inv(D(:,1:p)'*inv(Rxx)*D(:,1:p))*g; % MVDR weighting vector mvdr(ite,:)=pattern(n*d,w,ny); sir_end=sir_end+Sir(D,pj,w,L,p); %%%%%%% MVDR-LMS %%%%%%% w_int=D(:,1:p)*inv(D(:,1:p)'*D(:,1:p))*g; sir_int(1)=Sir(D,pj,w_int,L,p); % wq=w_int; % P=eye(N)-D(:,1:p)*inv(D(:,1:p)'*D(:,1:p))*D(:,1:p)'; for data=2:sir_data error=source(1,data)-w_int'*y(:,data); w_int=w_int+step*conj(error)*y(data); % Haykin % w_int=P*(w_int-satep*y(:,data)*conj(w_int'*y(:,data)))+wq; % Frost sir_int(data)=Sir(D,pj,w_int,L,p); end mvdr_lms(ite,:)=pattern(n*d,w_int,ny); ave_sir=ave_sir+sir_int; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% mvdr=abs(mvdr).^2; pat=Normalization(mvdr,experimentno); mvdr_lms=abs(mvdr_lms).^2; pat2=Normalization(mvdr_lms,experimentno); ave_sir=10*log10(ave_sir/experimentno); sir_end=10*log10(sir_end/experimentno); s1=mat2str(SNR'); s2=mat2str(Angs1'); figure k=sin([-90:0.5:90]'*pi/180); k1=[1:1:sir_data]; plot(k,pat,'.-g',k,pat2,'.-b'); legend('MVDR','MVDR-LMS'); xlabel('Angles'); ylabel('Power Gain in dB') title(sprintf('Array elements# =%g, SNR=%sdB from degree %s, experiment# =%g, snapshot=%g\nOTIR of MVDR=%gdB, OTIR of MVDR-LMS after %g iteration=%gdB',N,s1,s2,experimentno,snapshot,sir_end,sir_data,sir_int(end))) for u=1:L line(Angs1(u),[-70:0.1:0]) end line([-1:0.002:1],0); axis([-1 1 -70 10]); totaltime=now-starttime; tt=datestr(totaltime,13)