MATLAB Examples

Figure 52 (b). Displaced-beam (db) and Adjacent-beam (ab) pre-Doppler STAP performance with ambiguous clutter (β=3) and no jamming.

Contents

Coded by Ilias Konsoulas, 16 Sept. 2018. Code provided for educational purposes only. All rights reserved.

clc; clear; close all;

Radar System Operational Parameters

radar_oper_params;

Thermal Noise Power Computation

thermal_noise_power;

Thermal Noise Covariance Matrix

Rn = sigma2*eye(M*N);

Clutter Patch RCS computation

clutter_patch_rcs;

Calculate the Array Transmit and Element Receive Power Gains

Tx_Rx_power_gains;

Calculate the Clutter to Noise Ratio (CNR) for each azimuth angle

ksi = Pt*Gtgain.*Grec*lambda^2*sigma/((4*pi)^3*Pn*10^(Ls/10)*Rcik^4);   % Eq. (58)

Clutter Covariance Matrix Computations

beta = 3;   % beta parameter.
phia = 0;   % Velocity Misalignment Angle.

Rc = clutt_cov(ksi,beta);

Analytic Interference Covariance Matrix for Fully Adaptive STAP:

Ru = Rc + Rn;
InvRu = inv(Ru);

Target Space-Time Steering Vector

phit = 0; thetat = 0;      % Target azimuth and elevation angles in degrees.
fdt = 100;                 % Target Doppler Frequency.
fspt = d/lambda*cos(thetat*pi/180)*sin(phit*pi/180);
omegat = fdt/fr;
at = exp(-1i*2*pi*fspt*(0:N-1)).';     % Target Spatial Steering Vector.
fd = 0:.5:300;   Lfd = length(fd);
omegad = fd/fr;
SNRo = M*N;

LSINR Computation for Optimum Fully Adaptive Case

  LSINRopt = zeros(1,Lfd);
    for n=1:Lfd
         bt = exp(-1i*2*pi*omegad(n)*(0:M-1)).'; % Target Doppler Steering Vector
         vt = kron(bt,at);
         w = InvRu*vt;
         LSINRopt(n) = real(w'*vt)/SNRo;
    end

A. Displaced Beam pre-Doppler Calculations

Kt = 2;            % sub-CPI length (fixed).
Ks = 2:6;          % various sub-apperture lengths or number of beams used.
M1 = M - Kt +1;    % Number of possible sub-CPI's.

Create Doppler Filterbank Matrix F:

dopplerfilterbank = linspace(0,300,M1+1);
omegadopplerbank = dopplerfilterbank/fr;
F0 = zeros(M1,M1);
for m1=1:M1
    F0(:,m1) = 1/sqrt(M1)*exp(-1i*2*pi*omegadopplerbank(m1)*(0:M1-1));
end

F30 = diag(chebwin(M1,30))*F0;

Beamformer G Matrix Construction for Displaced Beam pre-Doppler

tb = [1;   1;];                      % Temporal Binomial Taper.
bt = [1;  -1;];                      % Doppler Steering Vector.
ta = ones(N,1);
% at = ones(N,1);                    % Steering Vector Pointing Broadside (phi = 0);
wmdb0   = zeros(M*N,M1,length(Ks));
wmdb30 = zeros(M*N,M1,length(Ks));

for n=1:length(Ks)
    N1 = N - Ks(n) + 1;              % Number of possible sub-appertures.
    g0 = exp(-1i*2*pi*fspt*(0:N1-1)).';
    g30 = chebwin(N1,30).*g0;        % 30-dB spatial taper.
    G0   = toeplitz([g0; zeros(Ks(n)-1,1);], [g0(1) zeros(1,Ks(n)-1)]);   % N1 x Ks Beamformer matrix G.
    G30 = toeplitz([g30; zeros(Ks(n)-1,1);], [g30(1) zeros(1,Ks(n)-1)]);  % N1 x Ks Beamformer matrix G.
    utab0 = kron(tb.*bt,G0'*(ta.*at));                                    % Desired response. Eq. 236.
    utab30 = kron(tb.*bt,G30'*(ta.*at));                                  % Desired response. Eq. 236.
    W0 = zeros(M*Ks(n),M1);
    W30 = zeros(M*Ks(n),M1);

W matrix construction.

    for p=0:M1-1
        Jp = [zeros(p,Kt); eye(Kt); zeros(M-Kt-p,Kt)];   % Selector Matrix Jp.
        Rupdb0  = kron(Jp,G0)'*Ru*kron(Jp,G0);
        Rupdb30 = kron(Jp,G30)'*Ru*kron(Jp,G30);
        wp0 = Rupdb0\utab0;
        wp30 = Rupdb30\utab30;
        Wp0 = reshape(wp0,Ks(n),Kt);
        Wp30 = reshape(wp30,Ks(n),Kt);
        aux1 = Wp0*Jp.';
        aux2 = Wp30*Jp.';
        W0(:,p+1)  = aux1(:);
        W30(:,p+1)  = aux2(:);
    end

    for m=1:M1
        wmdb0(:,m,n)  = kron(eye(M),G0)*W0*F0(:,m);
        wmdb30(:,m,n) = kron(eye(M),G30)*W30*F30(:,m);
    end
end

LSINR Computations for Displaced Beams

SINR0_cube   = zeros(M1,Lfd,length(Ks));
SINR30_cube = zeros(M1,Lfd,length(Ks));
SINR0 = zeros(length(Ks),Lfd);
SINR30 = zeros(length(Ks),Lfd);
for n=1:length(Ks)
    for m=1:M1
        for n1=1:Lfd
            bt1 = exp(-1i*2*pi*omegad(n1)*(0:M-1)).'; % Dummy Target Doppler Steering Vector
            vt = kron(bt1,at);
            SINR0_cube(m,n1,n)  = abs(wmdb0(:,m,n)'*vt)^2/real(wmdb0(:,m,n)'*Ru*wmdb0(:,m,n));
            SINR30_cube(m,n1,n) = abs(wmdb30(:,m,n)'*vt)^2/real(wmdb30(:,m,n)'*Ru*wmdb30(:,m,n));
        end
    end
    SINR0(n,:)  = max(abs(SINR0_cube(:,:,n)));
    SINR30(n,:) = max(abs(SINR30_cube(:,:,n)));
end

LSINR0  = SINR0/SNRo;
LSINR30 = SINR30/SNRo;

B. Adjacent Beams pre-Doppler Calculations

selbeamdist = 6.35;             % Selected Beam Angular Distance in degrees.
beamangles1 = [-9:-1  0:8]*selbeamdist; thetabeam = theta;

% beamangles1 = linspace(-4.5,4,N); thetabeam = theta;      % Inc = .5.
% beamangles1 = linspace(-9,8,N);   thetabeam = theta;      % Inc = 1.
% beamangles1 = linspace(-18,16,N); thetabeam = theta;      % Inc = 2.
% beamangles1 = linspace(-27,24,N); thetabeam = theta;      % Inc = 3.
% beamangles1 = linspace(-36,32,N); thetabeam = theta;      % Inc = 4.
% beamangles1 = linspace(-45,40,N); thetabeam = theta;      % Inc = 5.
% beamangles1 = linspace(-54,48,N); thetabeam = theta;      % Inc = 6.
% beamangles1 = linspace(-63,56,N); thetabeam = theta;      % Inc = 7.
% beamangles1 = linspace(-72,64,N); thetabeam = theta;      % Inc = 8.
% beamangles1 = linspace(-81,72,N); thetabeam = theta;      % Inc = 9.
% beamangles1 = linspace(-90,80,N); thetabeam = theta;      % Inc = 10.
beamangleinc = beamangles1(2) - beamangles1(1);
beamangles2 = beamangles1 + beamangleinc/2;

beamfreqs1 = d/lambda*cos(thetabeam*pi/180)*sin(beamangles1*pi/180);
beamfreqs2 = d/lambda*cos(thetabeam*pi/180)*sin(beamangles2*pi/180);

Godd  = zeros(N,N);
Geven = zeros(N,N);

for n=1:N
    Godd(:,n)  = 1/sqrt(N)*exp(-1i*2*pi*beamfreqs1(n)*(0:N-1));
    Geven(:,n) = 1/sqrt(N)*exp(-1i*2*pi*beamfreqs2(n)*(0:N-1));
end

wmab0  = zeros(M*N,M1,length(Ks));
wmab30 = zeros(M*N,M1,length(Ks));

for n=1:length(Ks)
    N1 = N - Ks(n) + 1;                                              % Number of possible sub-appertures.
    P = floor(Ks(n)/2);
    if mod(Ks(n),2)  % if Ks is odd
        Js = [zeros(N/2-P,Ks(n)); eye(Ks(n)); zeros(N/2-P-1,Ks(n))]; % Beam Selector Matrix Js.
        Gab0 = Godd*Js;                                              % N x Ks Beamformer matrix.
        Gab30 = diag(chebwin(N,30))*Gab0;
    else
        Js = [zeros(N/2-P,Ks(n)); eye(Ks(n)); zeros(N/2-P,Ks(n))];   % Beam Selector Matrix Js.
        Gab0 = Geven*Js;
        Gab30 = diag(chebwin(N,30))*Gab0;% N x Ks Beamformer matrix.
    end

W matrix construction

    ut2ab0   = kron(tb.*bt,Gab0'*(ta.*at));                          % Desired response. Eq. 236.
    ut2ab30 = kron(tb.*bt,Gab30'*(ta.*at));                          % Desired response. Eq. 236.
    Wab0 = zeros(M*Ks(n),M1);
    Wab30 = zeros(M*Ks(n),M1);

    for p=0:M1-1
        Jp = [zeros(p,Kt); eye(Kt); zeros(M-Kt-p,Kt)];               % Selector Matrix Jp.
        Rup0   = kron(Jp,Gab0)'*Ru*kron(Jp,Gab0);
        Rup30 = kron(Jp,Gab30)'*Ru*kron(Jp,Gab30);
        wpab0 = Rup0\ut2ab0;
        wpab30 = Rup30\ut2ab30;
        Wpab0 = reshape(wpab0,Ks(n),Kt);
        Wpab30 = reshape(wpab30,Ks(n),Kt);
        aux1 = Wpab0*Jp.';
        aux2 = Wpab30*Jp.';
        Wab0(:,p+1)  = aux1(:);
        Wab30(:,p+1) = aux2(:);
    end

    for m=1:M1
        wmab0(:,m,n)  = kron(eye(M),Gab0)*Wab0*F0(:,m);
        wmab30(:,m,n) = kron(eye(M),Gab30)*Wab30*F30(:,m);
    end
end

LSINR Computations for Adjacent Beams

SINR0ab_mat  = zeros(M1,Lfd,length(Ks));
SINR30ab_mat = zeros(M1,Lfd,length(Ks));
SINR0ab = zeros(length(Ks),Lfd);
SINR30ab = zeros(length(Ks),Lfd);
for n=1:length(Ks)
    for m=1:M1
        for n1=1:Lfd
            bt = exp(-1i*2*pi*omegad(n1)*(0:M-1)).'; % Dummy Target Doppler Steering Vector
            vt = kron(bt,at);
            SINR0ab_mat(m,n1,n)  = abs(wmab0(:,m,n)'*vt)^2/real(wmab0(:,m,n)'*Ru*wmab0(:,m,n));
            SINR30ab_mat(m,n1,n) = abs(wmab30(:,m,n)'*vt)^2/real(wmab30(:,m,n)'*Ru*wmab30(:,m,n));
        end
    end
    SINR0ab(n,:)  = max(abs(SINR0ab_mat(:,:,n)));
    SINR30ab(n,:) = max(abs(SINR30ab_mat(:,:,n)));
end

LSINR0ab  = SINR0ab/SNRo;
LSINR30ab = SINR30ab/SNRo;

Plot the SINR Losses

figure('NumberTitle', 'off','Name', ...
      ['Figure 52. (b) Displaced-beam pre-Doppler STAP performance with doppler-ambiguous clutter (β = ',...
                                                               num2str(beta),') and no jamming.'],...
       'Position',[1 1 1100 850]);

%  Displaced Beams Results
subplot(2,2,1);
%    plot(fd,10*log10(LSINRopt),'LineWidth',1.5)
%    hold on;
plot(fd,10*log10(LSINR0),'LineWidth',1.5)
grid on;
ylabel('SINR Loss (dB)');
xlabel('Target Doppler Frequency (Hz)');
ylim([-35 0]);
xlim([-5 305]);
legend('2 Beams', '3 Beams', '4 Beams', '5 Beams', '6 Beams','Location','South');
title('Displaced Beams, Untapered');

subplot(2,2,3);
%   plot(fd,10*log10(LSINRopt),'LineWidth',1.5)
%   hold on;
plot(fd,10*log10(LSINR30),'LineWidth',1.5)
grid on;
ylabel('SINR Loss (dB)');
xlabel('Target Doppler Frequency (Hz)');
ylim([-35 0]);
xlim([-5 305]);
legend('2 Beams', '3 Beams', '4 Beams', '5 Beams', '6 Beams','Location','South');
title('Displaced Beams, 30 dB Chebyshev Taper');

% Adjacent Beams Results
subplot(2,2,2);
%   plot(fd,10*log10(LSINRopt),'LineWidth',1.5)
%   hold on;
plot(fd,10*log10(LSINR0ab),'LineWidth',1.5)
grid on;
ylabel('SINR Loss (dB)');
xlabel('Target Doppler Frequency (Hz)');
ylim([-35 0]);
xlim([-5 305]);
legend('2 Beams', '3 Beams', '4 Beams', '5 Beams', '6 Beams','Location','South');
title('Adjacent Beams, Untapered');

subplot(2,2,4);
%   plot(fd,10*log10(LSINRopt),'LineWidth',1.5)
%   hold on;
plot(fd,10*log10(LSINR30ab),'LineWidth',1.5)
grid on;
ylabel('SINR Loss (dB)');
xlabel('Target Doppler Frequency (Hz)');
ylim([-35 0]);
xlim([-5 305]);
legend('2 Beams', '3 Beams', '4 Beams', '5 Beams', '6 Beams','Location','South');
title('Adjacent Beams, 30 dB Chebyshev Taper');