lcmv matlab,matlab宽带时域LCMV算法

这是我通过看文献自己编的一个Forst结构的宽带LCMV算法，可是运行出来没有主波束，不知道哪里出了问题，请大神指点，谢谢！！

clc;clear all;close all;

ima=sqrt(-1);c=1500;fs=50000;

fl=1000; fu=5000;

k=21;      %子带数量

Nf = linspace(fl,fu, k); %各点频的频率值1*k

L=16;              %FIR滤波器长度

%%阵列参数设定%%

N=16; %阵元数

theta=-90:1:90;  %搜索范围确定

f=fl:200:fu;

thetas=10;                 %来波方向

thetai=50;                  %干扰方向

Ns=5000;  %采样点数

d=1/2*c/fu;  %阵元间距与波长的比值

n=[0:N-1]';

%  带通滤波器

Wn1=2*[1000 5000]/fs;

bb=fir1(128,Wn1,'bandpass');

as=1*randn(1,Ns);%目标信号

NN=500*randn(1,Ns);

ai=filter(bb,1,NN);%干扰信号

%ai=1*randn(1,Ns);

noise=0.1*(randn(N,Ns)+ima*randn(N,Ns));%环境噪声

for i=1:k

for l=1:L

ss=exp(-ima*2*pi*Nf(i)*(d/c*sin(thetas*pi/180)*n+(l-1)*1/fs));

si=exp(-ima*2*pi*Nf(i)*(d/c*sin(thetai*pi/180)*n+(l-1)*1/fs));

a1(:,l)=ss;

a2(:,l)=si;

%NN=500000*randn(1,Ns);

%ai=filter(bb,1,NN);

%as=5*randn(1,Ns);

s=ss*as;

in=si*ai;%N*Ns矩阵

%noise=1*(randn(N,Ns)+ima*randn(N,Ns));

x((1+16*(l-1):16+16*(l-1)),:)=s+in+noise;%干扰加噪声N*Ns

end

vs=reshape(a1,N*L,1);%目标信号的导向矢量

vi=reshape(a2,N*L,1);%干扰的导向矢量

C=[vs vi];%约束矩阵

F=[exp(-ima*pi*Nf(i)*(L-1)*1/fs);10^-6];%阵列响应矢量

%F=[1;0];

Rx=1/Ns*x*x';%协方差矩阵

Wopt=Rx\C/(C'/Rx*C)*F;%权向量

%%%%%波束形成%%%%%%%%%%

for m=1:length(theta);

for mm=1:L

a=exp(-ima*2*pi*Nf(i)*(d/c*sin(theta(m)*pi/180)*n+(mm-1)*1/fs));

a0(:,mm)=a;

end

aa=reshape(a0,N*L,1);

p(i,m)=abs(Wopt'*aa);%17*91

end

end

% 归一化

figure

for k1=1:k

ppp(k1)=max(p(k1,:));

pp=max(p(k1,:));%第k1行的最大元素的值

out1=p(k1,:)./pp;

plot(theta,20.*log10(out1));

xlabel('方位角/度')

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

title('零陷')

hold on

end

grid on

P=20*log10(p./max(ppp));

figure

[X,Y]=meshgrid(theta,f);

mesh(X,Y,P);

ylabel('频率/Hz'),xlabel('角度/°'),zlabel('归一化方向图增益/db')

%figure

%pcolor(X,Y,P);shading interp;

%ylabel('频率/Hz'),xlabel('角度/°')

%hold on;

grid on;