yuanzhenolen1.m

圆震波束形成程序

%  圆阵 olen 优化
 clc;
clear;
N=8;
f=2000;
fs=20*f;
c=1500;
lmda=c/f;
r=0.5*lmda;%-------
T=0.2;
snr=20;
sita0=40*pi/180;
t=1/fs:1/fs:T;
x=sin(2*pi*f*t);
 ps=sum(abs(x).^2)/length(t);
 pn=ps*10^(-snr/10);
x1=[zeros(1,500),x,zeros(1,500)];
n=1:N;
  tau0(n)=r*cos(sita0-2*pi*(n-1)/N)/c;
a=exp(j*2*pi*f*tau0);  

xj=a'*x+sqrt(pn)*randn(N,length(t));
sita=(-180:180)*pi/180;


R=zeros(N,N);
sitam=[-180:1:180]*pi/180;kshi=zeros(length(sitam));

for m=1:length(sitam)
    a1(m,:)=exp(-j*2*pi*f*r*cos(sitam(m)-2*pi*(n-1)/N)/c);
    R=R+kshi(m)*(a1(m,:)'*a1(m,:));
end
mu=2;
fai=pn*(eye(N)+R);
 wn=mu*inv(fai)*a.';    
for m=1:length(sita)
    a2(m,:)=exp(j*2*pi*f*r*cos(sita(m)-2*pi*(n-1)/N)/c);
end
for m=1:length(sita)
    p(m)=abs(wn.'*a2(m,:)');
end
p=p/max(p);
pp=p;
% toc
% figure
% plot(sita*180/pi,abs(pp)/max(abs(pp)),'o-')
%   grid on
%   figure
%   plot(sita*180/pi,20*log10(pp))
%  grid on

%----------------------
% 主瓣宽度  迭代次数
k=1;
d1=1/(10^(1.5));  %D/20% 旁瓣高度 
R1=zeros(N,N);
 for m=1:length(sitam)
        if sitam(m)>=-25*pi/180&sitam(m)<=110*pi/180
           kshi1(m)=0;
        else
            gama(m)=kshi(m)+k*(p(m)-d1);
            kshi1(m)=max(0,gama(m));
        end
 end
        for m=1:length(sitam)
            a1(m,:)=exp(j*2*pi*f*r*cos(sitam(m)-2*pi*(n-1)/N)/c);
            R1=R1+kshi1(m)*(a1(m,:)'*a1(m,:));
        end

    fai=pn*(eye(N)+R1);
    wn1=mu*inv(fai)*a.';
        for m=1:length(sita)
         p(m)=abs(wn1.'*a2(m,:)');
        end
       p=p/max(p); 

ii=1;
while any(abs(real(wn1-wn)./real(wn))>0.001)|any(abs(imag(wn1-wn)./imag(wn))>0.001)
    ii;
    ii=ii+1;
    wn=wn1;
    R1=zeros(N,N);
        for m=1:length(sitam)
 
        if sitam(m)>=-25*pi/180&sitam(m)<=110*pi/180
            kshi1(m)=0;
        else
            gama(m)=kshi(m)+k*(p(m)-d1);
            kshi1(m)=max(0,gama(m));
        end
        end
        for m=1:length(sitam)
            m;
            a1(m,:)=exp(-j*2*pi*f*r*cos(sitam(m)-2*pi*(n-1)/N)/c);
            R1=R1+kshi1(m)*(a1(m,:)'*a1(m,:));
        end
    mu=1;
    fai=pn*(eye(N)+R1);
    wn1=mu*inv(fai)*a.';
        for m=1:length(sita)
            p(m)=abs(wn1.'*a2(m,:)');
        end
       p=p/max(p); 
       kshi=kshi1;    
end
figure

plot(sita*180/pi,20*log10(pp)),grid on,ylabel('归一化输出功率/db')

figure,plot(sita*180/pi,20*log10(p)),grid on,ylabel('归一化输出功率/db')