algrithem_database.m

基于Matlab开发的数字波束形成算法（含界面和6个小算法）

% hObject    handle to radiobutton6 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of radiobutton6

function PCI_Callback(hObject, eventdata, handles)
% hObject    handle to radiobutton6 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of radiobutton6

function PCA_MVB_Callback(hObject, eventdata, handles)
% hObject    handle to radiobutton6 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of radiobutton6

function SC_MVB_Callback(hObject, eventdata, handles)
% hObject    handle to radiobutton6 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of radiobutton6

function EC_Callback(hObject, eventdata, handles)
% hObject    handle to radiobutton6 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of radiobutton6

function reset_Callback(hObject, eventdata, handles)
% hObject    handle to reset (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

initialize_gui(gcbf, handles, true);     %改动处


function initialize_gui(fig_handle, handles, isreset)
% If the metricdata field is present and the reset flag is false, it means
% we are we are just re-initializing a GUI by calling it from the cmd line
% while it is up. So, bail out as we dont want to reset the data.
if isfield(handles, 'metricdata') && ~isreset
    return;
end

handles.metricdata.array_number = 16;
handles.metricdata.work_frequency  = 3e9;
handles.metricdata.array_distance  = 0.05;
handles.metricdata.work_wavelength  = 0.1;

handles.metricdata.SNR  = 0;
handles.metricdata.JNR1  = 30;
handles.metricdata.JNR2  = 30;
handles.metricdata.JNR3  = 30;

handles.metricdata.ang0  = 15;
handles.metricdata.ang1  = 25;
handles.metricdata.ang2  = 0;
handles.metricdata.ang3  = -15;

handles.metricdata.snap_number  = 100;
handles.metricdata.sample_frequency  = 18e9;

set(handles.array_number, 'String', handles.metricdata.array_number);
set(handles.work_frequency,  'String', handles.metricdata.work_frequency);
set(handles.array_distance, 'String', handles.metricdata.array_distance);
set(handles.work_wavelength,  'String', handles.metricdata.work_wavelength);

set(handles.SNR, 'String', handles.metricdata.SNR);
set(handles.JNR1,  'String', handles.metricdata.JNR1);
set(handles.JNR2, 'String', handles.metricdata.JNR2);
set(handles.JNR3,  'String', handles.metricdata.JNR3);

set(handles.ang0, 'String', handles.metricdata.ang0);
set(handles.ang1,  'String', handles.metricdata.ang1);
set(handles.ang2, 'String', handles.metricdata.ang2);
set(handles.ang3,  'String', handles.metricdata.ang3);

set(handles.snap_number, 'String', handles.metricdata.snap_number);
set(handles.sample_frequency,  'String', handles.metricdata.sample_frequency);

set(handles.AM_JAM1,'Value',1,'Enable','on');
set(handles.PM_JAM1,'Value',0,'Enable','on');
set(handles.AM_JAM2,'Value',1,'Enable','on');
set(handles.PM_JAM2,'Value',0,'Enable','on');
set(handles.AM_JAM3,'Value',1,'Enable','on');
set(handles.PM_JAM3,'Value',0,'Enable','on');

set(handles.MVDR,'Value',1,'Enable','on');
set(handles.GSC,'Value',0,'Enable','on');
set(handles.MWF,'Value',0,'Enable','on');
set(handles.LCEC,'Value',0,'Enable','on');
set(handles.PCI,'Value',0,'Enable','on');
set(handles.PCA_MVB,'Value',0,'Enable','on');
set(handles.SC_MVB,'Value',0,'Enable','on');
set(handles.EC,'Value',0,'Enable','on');
% Update handles structure
guidata(handles.figure1, handles);

function Calculate_Callback(hObject, eventdata, handles)
% hObject    handle to Calculate (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
M=handles.metricdata.array_number;
fc=handles.metricdata.work_frequency;
d=handles.metricdata.array_distance;
lamda=handles.metricdata.work_wavelength;

SNR=handles.metricdata.SNR;
JNR1=handles.metricdata.JNR1;
JNR2=handles.metricdata.JNR2;
JNR3=handles.metricdata.JNR3;

ang0=handles.metricdata.ang0;
ang1=handles.metricdata.ang1;
ang2=handles.metricdata.ang2;
ang3=handles.metricdata.ang3;

K=handles.metricdata.snap_number;
fs=handles.metricdata.sample_frequency;

a0=zeros(M,1);        % a0 期望信号方向矢量
a1=zeros(M,1);        % a1 干扰信号1的方向矢量
a2=zeros(M,1);        % a2 干扰信号2的方向矢量
a3=zeros(M,1);        % a3 干扰信号3的方向矢量
j=sqrt(-1);
for i=1:M
    a0(i)=exp(j*2*pi/lamda*d*(i-1)*sin(ang0*pi/180));
    a1(i)=exp(j*2*pi/lamda*d*(i-1)*sin(ang1*pi/180));
    a2(i)=exp(j*2*pi/lamda*d*(i-1)*sin(ang2*pi/180));
    a3(i)=exp(j*2*pi/lamda*d*(i-1)*sin(ang3*pi/180));
end
% A=[a0';a1';a2';a3']';
% A1=[a1';a2';a3']';
A=[a0,a1,a2,a3];
A1=[a1,a2,a3];

t=(1:1:K)/fs; 

s0=cos(2*pi*fc*t);        

AM_JAM1=get(handles.AM_JAM1,'value');
PM_JAM1=get(handles.PM_JAM1,'value');
AM_JAM2=get(handles.AM_JAM2,'value');
PM_JAM2=get(handles.PM_JAM2,'value');
AM_JAM3=get(handles.AM_JAM3,'value');
PM_JAM3=get(handles.PM_JAM3,'value');


if AM_JAM1==1
    s01=randn(1,K).*sin(2*pi*fc*t+12);
end

if PM_JAM1==1
    s01=sin(2*pi*fc*t+randn(1,K));
end

if AM_JAM2==1
    s02=randn(1,K).*sin(2*pi*fc*t+12);
end

if PM_JAM2==1
    s02=sin(2*pi*fc*t+randn(1,K));
end

if AM_JAM3==1
    s03=randn(1,K).*sin(2*pi*fc*t+12);
end

if PM_JAM3==1
    s03=sin(2*pi*fc*t+randn(1,K));
end


Ps0=s0*s0'/K;                       % 计算期望信号的功率
Ps1=s01*s01'/K;                     % 计算干扰信号1的功率
Ps2=s02*s02'/K;                     % 计算干扰信号2的功率
Ps3=s03*s03'/K;                     % 计算干扰信号2的功率

Ps=10.^(SNR/10);
Pj1=10.^(JNR1/10);
Pj2=10.^(JNR2/10);
Pj3=10.^(JNR3/10);

s=s0/Ps0^0.5*Ps^0.5;             % 具有输入信噪比的期望信号
s1=s01/Ps1^0.5*Pj1^0.5;          % 具有干噪比1的干扰信号1
s2=s02/Ps2^0.5*Pj2^0.5;          % 具有干噪比2的干扰信号2
s3=s03/Ps3^0.5*Pj3^0.5;          % 具有干噪比2的干扰信号2

S=[s;s1;s2;s3];
S1=[s1;s2;s3];
numJam=size(S1,1);
Noise=randn(M,K);
X=A*S+Noise;
X1=A1*S1+Noise;

R=X*X'/K;
Rjn=X1*X1'/K;


h0=a0/(a0'*a0)^0.5;
B0=zeros(M-1,M);
for i=1:M-1
    for k=i:i+1
        B0(i,k)=-1^(i+k)/h0(k);
    end
end
Rx0=B0*Rjn*B0';
rx0d0=B0*Rjn*h0;
Wx0=inv(Rx0)*rx0d0;

%%%%%%%%%%%%%%MWF_SIM%%%%%%%%%%%%%%%
%前向递推
r=M-1;               % r MWF的阶数，三种算法结果完全一样，说明MWF与MVDR的等效性
h=[];
% Pd=[];
delta=[];
L=[];
% D=[];
B=eye(M-1);
rxd=rx0d0;
Rxi=Rx0;
for i=1:r
    if i<r
        deltai=(rxd'*rxd)^0.5;
        delta=[delta,deltai];
        hi=rxd/deltai;
        h=[h;hi];
        if i==1
            L=[L;hi'];
        else
            L=[L;hi'*B];
        end
        Bi=zeros(M-i-1,M-i);
        for kk=1:M-i-1
            for k=kk:kk+1
                Bi(kk,k)=-1^(kk+k)/hi(k);
            end
        end
        B=Bi*B;
        rxd=Bi*Rxi*hi;
%         Pdi=hi'*Rxi*hi;
%         Pd=[Pd,Pdi];
        Rxi=Bi*Rxi*Bi';
    else
        L=[L;B];
    end
end

Rd=L*Rx0*L';
rdd0=L*rx0d0;
Wd=inv(Rd)*rdd0;

[V,D]=eig(Rjn);
%对特征值从小到大进行排列并对对应的特征向量进行调整
for rr=1:M
for i=1:M-1
    if D(i,i)>D(i+1,i+1)
        Dtemp=D(i+1,i+1);
        D(i+1,i+1)=D(i,i);
        D(i,i)=Dtemp;
        Vtemp=V(:,i+1);
        V(:,i+1)=V(:,i);
        V(:,i)=Vtemp;
    end
end
end
Us=V(:,(M-numJam+1:M));
Ds=D((M-numJam+1:M),(M-numJam+1:M));
T=[Us,a0];

Un=V(:,(1:M-numJam));
Dn=D((1:M-numJam),(1:M-numJam));

dd=diag(D);
[dd,IX]=sort(dd);
dd=flipud(dd);
IX=flipud(IX);
U=[];
for i=1:r
    U=[U,V(:,IX(i))];
end

MVDR=get(handles.MVDR,'value');
GSC=get(handles.GSC,'value');
MWF=get(handles.MWF,'value');
LCEC=get(handles.LCEC,'value');
PCI=get(handles.PCI,'value');
PCA_MVB=get(handles.PCA_MVB,'value');
SC_MVB=get(handles.SC_MVB,'value');
EC=get(handles.EC,'value');

if MVDR==1
    W1=inv(Rjn)*a0/(a0'*inv(Rjn)*a0);
end

if GSC==1
    W1=h0-B0'*Wx0;
end

if MWF==1
    W1=h0-B0'*L'*Wd;
end

if LCEC==1
    W1=((eye(size(Us,1))-Us*Us')*a0)/(a0'*(eye(size(Us,1))-Us*Us')*a0);   % LCEC_SIM
end

if PCI==1
    W1=((eye(size(Us,1))-Us*Us')*a0);                                     % PCI
end

if PCA_MVB==1
    W1=(T*inv(T'*Rjn*T)*T'*a0)/(a0'*T*inv(T'*Rjn*T)*T'*a0);                     % PCA_MVB
end

if SC_MVB==1
    W1=(Un*inv(Dn)*Un'*a0)/(a0'*Un*inv(Dn)*Un'*a0);                  %SC_MVB
end

if EC==1
    W1=(eye(M)-U*U')*a0;     %  EC
end
    

ang=-90:0.5:90;
a=zeros(M,1);
for k=1:length(ang)
    for i=1:M
        a(i)=exp(j*2*pi/lamda*d*(i-1)*sin(ang(k)*pi/180));
    end
    beam1(k)=W1'*a;            
end
beam1=abs(beam1);
beam1=beam1/max(beam1);
beam1=20*log10(beam1);

axes(handles.directional_diagram_axes)
plot(ang,beam1)
set(handles.directional_diagram_axes,'XMinorTick','on')
grid on