untitled.m

这是一个宽带波束形成的matlab程序

clear all;
clc;
close all;
randn('state',0);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%参数设置%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Ns=512;%时域采样点数；
FF=180/180;%空域采样间隔(单位:度)
snapshots=512;
DB=-10;%信噪比；
B=10^(DB/20);
c=3*10^8;%光速；
%阵元设置
M=36;
%中心点的半径
r0=0;
r1=0.101;%第一圈的半径
r2=0.192;%第二圈的半径
r3=0.287;%第三圈的半径
r4=0.405;%第四圈的半径
r5=0.56;%第五圈的半径

         for q=1
            beta(q)=0;
            XX(q)=r0*cos(beta(q));
            YY(q)=r0*sin(beta(q));
         end
         for q=2:8
            beta(q)=2*pi*(q-2)/7;
            XX(q)=r1*cos(beta(q));
            YY(q)=r1*sin(beta(q));
        end
        for q=9:15
            beta(q)=2*pi*(q-9)/7;
            XX(q)=r2*cos(beta(q));
            YY(q)=r2*sin(beta(q));
        end
        for q=16:22
            beta(q)=2*pi*(q-16)/7;
            XX(q)=r3*cos(beta(q));
            YY(q)=r3*sin(beta(q));
        end
        for q=23:29
            beta(q)=2*pi*(q-23)/7;
            XX(q)=r4*cos(beta(q));
            YY(q)=r4*sin(beta(q));
        end
        for q=30:36
            beta(q)=2*pi*(q-30)/7;
            XX(q)=r5*cos(beta(q));
            YY(q)=r5*sin(beta(q));
        end
f0=2.0*10^9;                          %频率1.5G-2.4G；
lambda=c/f0;                          %波长；
frequency_band=9*10^8;              %带宽；
fm=frequency_band/2;                 %调制频率;
fs=5*10^9;                          %采样频率
Ts=1/fs;                             %采样周期；
t=Ts:Ts:Ns*Ts;                       %采样时间点
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%信号设置%%%%%%%%%%%%%%%%%%%%%%%
%信号频率
ff=[f0,f0,f0];
%信号个数
signalnum=length(ff);
%信号源DOA:[fai(0-90),sita(0-360)]
angle=[30 200
       10 100
       50 300];   
%生成阵列流行矩阵，并产生接收信号矩阵
for sig=1:signalnum
        S(sig,:)=B*(1+0.1*cos(2*pi*fm*t)).*exp(i*2*pi*ff(sig)*t)+(randn(1,Ns)+i*randn(1,Ns))/2^0.5;
end
for p=1:Ns
   % 生成阵列流行矩阵
    for k=1:M
        for sig=1:signalnum
            A0(k,sig)=exp(i*2*pi*ff(sig)*(cos(angle(sig,2)*pi/180)*sin(angle(sig,1)*pi/180)...
                     *XX(k)+sin(angle(sig,2)*pi/180)*sin(angle(sig,1)*pi/180)*YY(k))/c)*(cos(angle(sig,1)*pi/180))^0.5;
         end;
     end;
    %产生接收信号向量 
    x1=A0*S(:,p);
    Xs(:,p)=x1;                
end
%对接受信号进行混频处理，产生零中频信号
for k=1:M
    for p=1:Ns
        Xs(k,p)=Xs(k,p)*exp(-i*2*pi*f0*p*Ts);
    end
end
%对信号进行时域滤波
[b,a]=butter(5,(frequency_band)/(fs/2)); 
for k=1:M
        Ys(k,:)=filter(b,a,Xs(k,:));
end
Xs=Ys;
R=Xs*Xs';
%对角加载
L=60;
XX11=conj(Xs');
%%%%%%%%%%%%%%%%%%%%%%频域聚焦过程%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
XX1=fft(XX11,251);
XX2=conj(XX1');
figure 
plot([-125:125]/250*(fs/2),abs(fftshift(fft(XX11(:,1),251))));
%中心频率为2.0GHz
f0=2.0*10^9;
%%分频过程
ffd=10*10^6/250;  
fff(126)=f0;
for ffk=1:125
    fff(ffk)=f0-ffd*(125-ffk+1);
end
for ffk=127:251
    fff(ffk)=f0+ffd*(ffk-125-1);
end
%确定参考频率
f0=2.0*10^9;
%%分频产生阵列流行矩阵过程
for fffk=1:125                    
    for k=1:M
        for sig=1:signalnum
            A(k,sig)=exp(i*2*pi*fff(fffk)*(cos(angle(sig,2)*pi/180)*sin(angle(sig,1)*pi/180)...
                     *XX(k)+sin(angle(sig,2)*pi/180)*sin(angle(sig,1)*pi/180)*YY(k))/c)*(cos(angle(sig,1)*pi/180))^0.5;;
        end
    end
AA=A*A0';
[U S V]=svd(AA);
%求聚焦矩阵
T1=V*U';
%生成第40个频点的聚焦矩阵
if fffk==40
    TT=T1;
end
%生成第80个频点的聚焦矩阵
if fffk==80
    TT1=T1;
end
XX2(:,fffk)=T1*XX2(:,fffk);
end
for fffk=127:251                   
    for k=1:M
        for sig=1:signalnum
            A(k,sig)=exp(i*2*pi*fff(fffk)*(cos(angle(sig,2)*pi/180)*sin(angle(sig,1)*pi/180)...
                     *XX(k)+sin(angle(sig,2)*pi/180)*sin(angle(sig,1)*pi/180)*YY(k))/c)*(cos(angle(sig,1)*pi/180))^0.5;;
        end
    end
AA=A*A0';
[U S V]=svd(AA);
T2=V*U';
%生成第130个频点的聚焦矩阵
if fffk==130
    TT2=T2;
end
%生成第200个频点的聚焦矩阵
if fffk==200
    TT3=T2;
end
XX2(:,fffk)=T2*XX2(:,fffk);
end
%聚焦过程
R0=zeros(36,36);
for kkkk=1:251
XX0=XX2(:,kkkk);
%求相关矩阵
R0=conj(XX0')*(conj(XX0'))'+R0;
end
R0=R0/251;
f1=[1;0;0];
%对Rx进行对角加载
R0=R0+10^(L/20)*eye(M);
Wopt0=inv(R0)*(A0)*inv((A0'*inv(R0)*A0))*f1;
%不同频点处对应的权值
W1=T1'*Wopt0;
W2=T2'*Wopt0;
W3=TT'*Wopt0;
W4=TT1'*Wopt0;
W5=TT2'*Wopt0;
W6=TT3'*Wopt0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%波束输出%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
lambda=c/fff(251);
sig=1;
da=1;
p=1;
for theta=0:da:90
      for fai=0:da:360
        for q=1
            beta(q)=0;
            A11(q,p)=exp(j*(2*pi./lambda*r0.*sin(theta*pi/180).*cos(fai...
                   *pi/180-beta(q))));
        end
        for q=2:8
            beta(q)=2*pi*(q-2)/7;
            A11(q,p)=exp(j*(2*pi./lambda*r1.*sin(theta*pi/180).*cos(fai...
                   *pi/180-beta(q))));
        end
        for q=9:15
            beta(q)=2*pi*(q-9)/7;
            A11(q,p)=exp(j*(2*pi./lambda*r2.*sin(theta*pi/180).*cos(fai...
                   *pi/180-beta(q))));
        end
        for q=16:22
            beta(q)=2*pi*(q-16)/7;
            A11(q,p)=exp(j*(2*pi./lambda*r3.*sin(theta*pi/180).*cos(fai...
                   *pi/180-beta(q))));
        end
        for q=23:29
            beta(q)=2*pi*(q-23)/7;
            A11(q,p)=exp(j*(2*pi./lambda*r4.*sin(theta*pi/180).*cos(fai...
                   *pi/180-beta(q))));
        end
        for q=30:36
            beta(q)=2*pi*(q-30)/7;
            A11(q,p)=exp(j*(2*pi./lambda*r5.*sin(theta*pi/180).*cos(fai...
                   *pi/180-beta(q))));
        end
B11=A11*((cos(angle(sig,1)*pi/180))^0.5);
F(fai+1,theta+1)=abs(W2'*B11); 
    end
end
F1=F(201,:);
Fmax=max(max(F));
F1max=max(max(F1));
%%%%%%%%%%%%%%%%%%%%%%%%%%%方位面切面波束图输出%%%%%%%%%%%%%%%%%%%%%
lambda=c/fff(126);
sig=1;
da=1;
p=1;
for theta=0:da:90
      for fai=0:da:360
        for q=1
            beta(q)=0;
            A11(q,p)=exp(j*(2*pi./lambda*r0.*sin(theta*pi/180).*cos(fai...
                   *pi/180-beta(q))));
        end
        for q=2:8
            beta(q)=2*pi*(q-2)/7;
            A11(q,p)=exp(j*(2*pi./lambda*r1.*sin(theta*pi/180).*cos(fai...
                   *pi/180-beta(q))));
        end
        for q=9:15
            beta(q)=2*pi*(q-9)/7;
            A11(q,p)=exp(j*(2*pi./lambda*r2.*sin(theta*pi/180).*cos(fai...
                   *pi/180-beta(q))));
        end
        for q=16:22
            beta(q)=2*pi*(q-16)/7;
            A11(q,p)=exp(j*(2*pi./lambda*r3.*sin(theta*pi/180).*cos(fai...
                   *pi/180-beta(q))));
        end
        for q=23:29
            beta(q)=2*pi*(q-23)/7;
            A11(q,p)=exp(j*(2*pi./lambda*r4.*sin(theta*pi/180).*cos(fai...
                   *pi/180-beta(q))));
        end
        for q=30:36
            beta(q)=2*pi*(q-30)/7;
            A11(q,p)=exp(j*(2*pi./lambda*r5.*sin(theta*pi/180).*cos(fai...
                   *pi/180-beta(q))));
        end
B11=A11*((cos(angle(sig,1)*pi/180))^0.5);
F(fai+1,theta+1)=abs(Wopt0'*B11); 
    end
end
F2=F(201,:);
%第40个频点方位切面图输出
lambda=c/fff(40);
sig=1;
da=1;
p=1;
for theta=0:da:90
      for fai=0:da:360
        for q=1
            beta(q)=0;
            A11(q,p)=exp(j*(2*pi./lambda*r0.*sin(theta*pi/180).*cos(fai...
                   *pi/180-beta(q))));
        end
        for q=2:8
            beta(q)=2*pi*(q-2)/7;
            A11(q,p)=exp(j*(2*pi./lambda*r1.*sin(theta*pi/180).*cos(fai...
                   *pi/180-beta(q))));
        end
        for q=9:15
            beta(q)=2*pi*(q-9)/7;
            A11(q,p)=exp(j*(2*pi./lambda*r2.*sin(theta*pi/180).*cos(fai...
                   *pi/180-beta(q))));
        end
        for q=16:22
            beta(q)=2*pi*(q-16)/7;
            A11(q,p)=exp(j*(2*pi./lambda*r3.*sin(theta*pi/180).*cos(fai...
                   *pi/180-beta(q))));
        end
        for q=23:29
            beta(q)=2*pi*(q-23)/7;
            A11(q,p)=exp(j*(2*pi./lambda*r4.*sin(theta*pi/180).*cos(fai...
                   *pi/180-beta(q))));
        end
        for q=30:36
            beta(q)=2*pi*(q-30)/7;
            A11(q,p)=exp(j*(2*pi./lambda*r5.*sin(theta*pi/180).*cos(fai...
                   *pi/180-beta(q))));
        end
B11=A11*((cos(angle(sig,1)*pi/180))^0.5);
F(fai+1,theta+1)=abs(W3'*B11); 
    end
end
F3=F(201,:);
F3max=max(max(F3));
F2max=max(max(F2));
%第80个频点方位切面图输出
lambda=c/fff(80);
sig=1;
da=1;
p=1;
for theta=0:da:90
      for fai=0:da:360
        for q=1
            beta(q)=0;
            A11(q,p)=exp(j*(2*pi./lambda*r0.*sin(theta*pi/180).*cos(fai...
                   *pi/180-beta(q))));
        end
        for q=2:8
            beta(q)=2*pi*(q-2)/7;
            A11(q,p)=exp(j*(2*pi./lambda*r1.*sin(theta*pi/180).*cos(fai...
                   *pi/180-beta(q))));
        end
        for q=9:15
            beta(q)=2*pi*(q-9)/7;
            A11(q,p)=exp(j*(2*pi./lambda*r2.*sin(theta*pi/180).*cos(fai...
                   *pi/180-beta(q))));
        end
        for q=16:22
            beta(q)=2*pi*(q-16)/7;
            A11(q,p)=exp(j*(2*pi./lambda*r3.*sin(theta*pi/180).*cos(fai...
                   *pi/180-beta(q))));
        end
        for q=23:29
            beta(q)=2*pi*(q-23)/7;
            A11(q,p)=exp(j*(2*pi./lambda*r4.*sin(theta*pi/180).*cos(fai...
                   *pi/180-beta(q))));
        end
        for q=30:36
            beta(q)=2*pi*(q-30)/7;
            A11(q,p)=exp(j*(2*pi./lambda*r5.*sin(theta*pi/180).*cos(fai...
                   *pi/180-beta(q))));
        end
B11=A11*((cos(angle(sig,1)*pi/180))^0.5);
F(fai+1,theta+1)=abs(W4'*B11); 
    end
end
F4=F(201,:);
F4max=max(max(F4));
%第130个频点方位切面图输出
lambda=c/fff(130);
sig=1;
da=1;
p=1;
for theta=0:da:90
      for fai=0:da:360
        for q=1
            beta(q)=0;
            A11(q,p)=exp(j*(2*pi./lambda*r0.*sin(theta*pi/180).*cos(fai...
                   *pi/180-beta(q))));
        end
        for q=2:8
            beta(q)=2*pi*(q-2)/7;
            A11(q,p)=exp(j*(2*pi./lambda*r1.*sin(theta*pi/180).*cos(fai...
                   *pi/180-beta(q))));
        end
        for q=9:15
            beta(q)=2*pi*(q-9)/7;
            A11(q,p)=exp(j*(2*pi./lambda*r2.*sin(theta*pi/180).*cos(fai...
                   *pi/180-beta(q))));
        end
        for q=16:22
            beta(q)=2*pi*(q-16)/7;
            A11(q,p)=exp(j*(2*pi./lambda*r3.*sin(theta*pi/180).*cos(fai...
                   *pi/180-beta(q))));
        end
        for q=23:29
            beta(q)=2*pi*(q-23)/7;
            A11(q,p)=exp(j*(2*pi./lambda*r4.*sin(theta*pi/180).*cos(fai...
                   *pi/180-beta(q))));
        end
        for q=30:36
            beta(q)=2*pi*(q-30)/7;
            A11(q,p)=exp(j*(2*pi./lambda*r5.*sin(theta*pi/180).*cos(fai...
                   *pi/180-beta(q))));
        end
B11=A11*((cos(angle(sig,1)*pi/180))^0.5);
F(fai+1,theta+1)=abs(W5'*B11); 
    end
end
F5=F(201,:);
F5max=max(max(F5));
%第200个频点方位切面图输出
lambda=c/fff(200);
sig=1;
da=1;
p=1;
for theta=0:da:90
      for fai=0:da:360
        for q=1
            beta(q)=0;
            A11(q,p)=exp(j*(2*pi./lambda*r0.*sin(theta*pi/180).*cos(fai...
                   *pi/180-beta(q))));
        end
        for q=2:8
            beta(q)=2*pi*(q-2)/7;
            A11(q,p)=exp(j*(2*pi./lambda*r1.*sin(theta*pi/180).*cos(fai...
                   *pi/180-beta(q))));
        end
        for q=9:15
            beta(q)=2*pi*(q-9)/7;
            A11(q,p)=exp(j*(2*pi./lambda*r2.*sin(theta*pi/180).*cos(fai...
                   *pi/180-beta(q))));
        end
        for q=16:22
            beta(q)=2*pi*(q-16)/7;
            A11(q,p)=exp(j*(2*pi./lambda*r3.*sin(theta*pi/180).*cos(fai...
                   *pi/180-beta(q))));
        end
        for q=23:29
            beta(q)=2*pi*(q-23)/7;
            A11(q,p)=exp(j*(2*pi./lambda*r4.*sin(theta*pi/180).*cos(fai...
                   *pi/180-beta(q))));
        end
        for q=30:36
            beta(q)=2*pi*(q-30)/7;
            A11(q,p)=exp(j*(2*pi./lambda*r5.*sin(theta*pi/180).*cos(fai...
                   *pi/180-beta(q))));
        end
B11=A11*((cos(angle(sig,1)*pi/180))^0.5);
F(fai+1,theta+1)=abs(W6'*B11); 
    end
end
F6=F(201,:);
F6max=max(max(F6));
figure;
xx=0:da:360;
yy=0:da:90;
[Xx,Yy]=meshgrid(xx,yy);
mesh(Xx,Yy,20*log10(F'/Fmax));grid on;
xlabel('（度）');
ylabel('（度）');
zlabel('增益（dB）');
axis([0 360 0 90 -60 0]);
figure;
xx=0:da:360;
yy=0:da:90;
[Xx,Yy]=meshgrid(xx,yy);
mesh(Xx,Yy,20*log10(F'/Fmax));grid on;
view(0,0);
xlabel('方位角（度）');
zlabel('dB');
axis([0 360 0 90 -60 0]);
figure;
xx1=0:da:360;
yy1=0:da:90;
[X11,Y11]=meshgrid(xx1,yy1);
mesh(X11,Y11,20*log10(F'/Fmax));grid on;
view(90,0);
ylabel('俯仰角（度）');
zlabel('dB');
axis([0 360 0 90 -60 0]);
figure;
x=0:da:90;
plot(x,20*log10(F1'/F1max),x,20*log10(F2'/F2max),x,20*log10(F3'/F3max),x,20*log10(F4'/F4max),x,20*log10(F5'/F5max),x,20*log10(F6'/F6max));grid on;
xlabel('俯仰角（度）');
ylabel('dB');
axis([0 90 -60 0]);