%用均匀圆阵MVDR方法进行宽带波束形成

%% 基本参数

f0 = 4*10^6;       %信号中心频率

bandwide = 0.5*10^6;  %带宽0.05MHz

interval = 100;     %每间隔interval(Hz)形成点频信号，合成宽带

NarrowNumber = bandwide/interval;  %所需点频个数

NarrowF = linspace(f0-bandwide/2, f0+bandwide/2, NarrowNumber); %各点频的频率值

R = 50;             %阵元半径

c = 3*10^8;         %传播速度：光速

bosu=340;           % 波速

degrad=pi/180;

p=1;

seta1=150*degrad;                           %信号到来方位角

fai1=45*degrad;                            %信号到来仰角

kk = 64;                                   %进行MVDR所需采样数

Nfft = 256;                                %进行FFT子带化点数，即进行Nfft点快速傅立叶变换

Nlearn = kk * Nfft;                        %时域数据快拍数

M=40;                                      %均匀圆阵总数M

ElementFai=(0:(M-1))*9*degrad+4.5*degrad;  %各阵元位置向量图

MNumber = 10;                               %选出的阵元个数

SelectElement=[ElementFai(1) ElementFai(5) ElementFai(9) ElementFai(13) ElementFai(17) ElementFai(21) ElementFai(25) ElementFai(29) ElementFai(33) ElementFai(37)];

A=zeros(MNumber,p);

k=[0:MNumber-1]';

for t=1:p

A(:,t)=exp(-j*2*pi*R*f0/bosu*(cos(seta1*degrad)*cos(2*pi*k/M)+cos(seta1*degrad)*sin(2*pi*k/M)));%圆阵的阵列流型

end

s=zeros(p,Nlearn);

Am=[3,5];        %幅度

phi=[20,30];     %初相

for ii=1:p

s(ii,:)=s(ii,:)+Am(ii)*cos(2*pi*f0*t+phi(ii)*degrad);%宽带信号

end

F=linspace(f0-bandwide/2, f0+bandwide/2, Nfft);

snr=10;

n=exp(j*2*pi*randn(MNumber,Nlearn))/snr;%噪声

X=A*s + n;               % 阵列输出

%% 开始处理

%对时域数据分段作FFT变换

for mn=1:MNumber

for m=1:kk

xxF(mn,m,:) = fft( X(mn,((m-1)*Nfft+1):((m-1)*Nfft+Nfft)), Nfft);

end

end

temp1=cos(seta1-SelectElement);

F=linspace(f0-bandwide/2, f0+bandwide/2, Nfft);

%求窄带信号的自适应最佳权

for nfft=1:Nfft

h=zeros(MNumber);

for k=1:kk

h=h+xxF(:,k,nfft)*xxF(:,k,nfft)';

end

Rxxf(:,:,nfft)=h./kk;     %得到某窄频的自相关矩阵

windows=ones(MNumber,1);

a1=exp(-j*2*pi*R*F(nfft)/c*sin(fai1)*temp1);

s=(a1.*windows.').';

Wopt(:,nfft)=(inv(Rxxf(:,:,nfft))*s)/(s'*inv(Rxxf(:,:,nfft))*s);         %MVDR算法的程序

end

%% 扫描矢量

NcurveAzi=180;                     %空间扫描角度采样数

seta=linspace(0,2*pi,NcurveAzi);   %阵元空间扫描指向角度

for ne=1:Nfft

for na=1:NcurveAzi

temp=cos(seta(na)-SelectElement);

a=exp(-j*2*pi*R*F(ne)/c*sin(fai1)*temp);%扫描时方向矢量

f=(1+temp)/2;        %扫描时的阵元方向图

BF(na,ne)=abs(sum(a.*Wopt(:,nfft)'));

end

%%仿真结果

plot(seta*180/pi,20*log10(BF(:,1)/max(BF(:,1))),'b')

grid

xlabel('\phi方位角(度)')

ylabel('波束输出/分贝')

title('均匀圆阵MVDR算法宽带波束形成图')

end

hold off