Figure 64. SINR loss for Element Space pre-Doppler STAP. Zero intrinsic cluller motion.

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;

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)
Ksic = sigma2*diag(ksi);

Clutter Covariance Matrix Computations

K2 = 3;             % Number of Pulses per sub-CPI.
M2 = M - K2 + 1;    % Number of sub-CPI's for K2=3.
beta = 1;   % beta parameter.
phia = 0;   % Velocity Misalignment Angle.
Ita = d/lambda*cos(theta);

% Calculate Spatial and Doppler Frequencies.
% Spatial frequency of the k-th clutter patch.
fsp = Ita*sin(phi*pi/180);
% Normalized Doppler Frequency:
omegac = beta*Ita*sin(phi*pi/180 + phia*pi/180);

% Clutter Steering Vector Pre-allocation for sub-CPI of size K:
a = zeros(N,Nc);
b = zeros(M,Nc);
b2 = zeros(K2,Nc);
Vc  = zeros(M*N,Nc);
Vc2 = zeros(K2*N,Nc);

for k=1:Nc
    a(:,k) = exp(1i*2*pi*fsp(k)*(0:N-1));       % Spatial Steering Vector.
    b(:,k) = exp(1i*2*pi*omegac(k)*(0:M-1));    % Time Steering Vector.
    b2(:,k) = exp(1i*2*pi*omegac(k)*(0:K2-1));  % Time Steering Vector for K = 3.
    Vc(:,k) = kron(b(:,k),a(:,k));              % Space-Time Steering Vector.
    Vc2(:,k) = kron(b2(:,k),a(:,k));            % Space-Time Steering Vector for K = 3.
end

Rc = real(Vc*Ksic*Vc');
Rcsub2 = real(Vc2*Ksic*Vc2');

Thermal Noise Covariance Matrices

Rn = sigma2*eye(M*N);
Rnsub2 = sigma2*eye(K2*N);

Jammer Covariance Matrix Calculation.

Spatial frequency of the j-th jammer.

jamm_cov;

Rjsub2 = kron(eye(K2),Phi_j);

Beamformer G Matrix Construction for Displaced Filter Beamspace post-Doppler STAP with Two Step Nulling (TSN):

Ks = 3;
N1 = N - Ks + 1;         % Number of possible sub-appertures.
taper_level = 30;        % Tapering level in db.
phit = 0; thetat = 0;    % Target azimuth and elevation angles in degrees.
fspt = d/lambda*cos(thetat*pi/180)*sin(phit*pi/180);
g0 = exp(-1i*2*pi*fspt*(0:N1-1)).';
gtap = chebwin(N1,taper_level).*g0;           % spatial taper application.

% Spatial jammer-plus-noise Covariance Matrix:
Phi_jn = Phi_j + sigma2*eye(N);

% Calculation of Beamformer Matrix for Two-Step Nulling Approach:
GTSN = zeros(N,Ks);
for p=0:Ks-1
    Jp = [zeros(p,N1); eye(N1); zeros(N-N1-p,N1)]; % Beam Selector Matrix Jp.
    GTSN(p+1:N1+p,p+1) = (Jp.'*Phi_jn*Jp)\gtap;
end

Analytic Interference Covariance Matrix for Fully Adaptive STAP

Ru = Rc + Rj + Rn;

Analytic Interference Covariance Matrix for second sub-CPI (K=3)

Rusub2 = Rcsub2 + Rjsub2 + Rnsub2;

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);
omegact = fdt/fr;
at = exp(1i*2*pi*fspt*(0:N-1)).';   % Target Spatial Steering Vector.
ta = chebwin(N,30);                 % 30 dB Chebychev Spatial Tapper.
bt2 = [1; -1; 1];
tb2 = [1;  2; 1];                   % Binomial Doppler Taper for K = 3.
gt2 = kron(tb2.*bt2,ta.*at);

Tapered Fully Adaptive STAP Solution

wsub2 = Rusub2\gt2;   % Weight Vector for K = 3;

W2 = zeros(M*N,M2);
% W2 matrix construction. Equation (180).
% Assume that weight vectors for each sub-CPI are equal,.
for k=1:M2
    W2((k-1)*N+1:(k+2)*N,k) = wsub2;
end

td2 =  chebwin(M2,30);             % 30 dB Chebychev Doppler Taper.

1. SINR Computations for Optimum Fully Adaptive STAP:

phit = 0;                          % Target Azimuth Angle.
% fspt = d/lambda*cos(theta*pi/180)*sin(phi*pi/180);
fspt = 0;
at = exp(1i*2*pi*fspt*(0:N-1)).';  % Target Spatial Steering Vector.

fd = 0:.5:300;   Lfd = length(fd);
omegad = fd/fr;
dopplerfilterbank2 = linspace(0,300,M2+1);
omegadopplerbank2 = dopplerfilterbank2/fr;
LSINRopt = zeros(1,Lfd);
SINRsub2_mat = zeros(length(dopplerfilterbank2),Lfd);

InvRu = inv(Ru);
InvRusub2 = inv(Rusub2);
SNRo = M*N;

% LSINR Computation for Optimum Fully Adaptive Case:
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; %#ok<MINV>
    LSINRopt(n) = real(w'*vt)/SNRo;
end

2. SINR Computations for Element-Space pre-Doppler sub-CPI with K=3

for m=1:length(dopplerfilterbank2)
    um = exp(1i*2*pi*omegadopplerbank2(m)*(0:M2-1)).';
    fm2 = td2.*um;
    wm2 = W2*fm2;              % m-th Doppler bin composite weight vector.

    for n=1:Lfd
        bt = exp(1i*2*pi*omegad(n)*(0:M-1)).';  % Dummy Target Doppler Steering Vector
        vt = kron(bt,at);
        SINRsub2_mat(m,n) = abs(wm2'*vt)^2/real(wm2'*Ru*wm2);
    end
end

SINRsub2 = max(abs(SINRsub2_mat));
LSINRsub2 = SINRsub2/SNRo;

3. Element-space post-Doppler: PRI-Staggered Post-Doppler method

% Doppler Filter Bank Creation
dopplerfilterbank = linspace(0,300,M+1);
omegadopplerbank = dopplerfilterbank/fr;
fd = 0:.5:300;   Lfd = length(fd);
omegad = fd/fr;
SNRo = M*N;

K = 3;
P = floor(K/2);
M1= M - K +1;
U1 = zeros(M1,M);
for m=1:M
    U1(:,m) =  1/sqrt(M)*exp(-1i*2*pi*omegadopplerbank(m)*(0:M1-1));
end

td30 = chebwin(M1,30);        % 30-dB Chebyshev Doppler Taper.
F30 = diag(td30)*U1;

% Create Doppler Filter Bank in Fm Matrix for PRI-Staggered Post-Doppler method:
Fm30 = zeros(M,K,M);

for m=1:M
    Fm30(:,:,m) = toeplitz([F30(:,m); zeros(K-1,1)],[F30(1,m) zeros(1,K-1)]);
end

SINR30_mat   = zeros(M,Lfd);
SINRab30_mat = zeros(M,Lfd);

% PRI-Staggered Post-Doppler SINR Computations:
for m=1:M

    f30m = Fm30(:,:,m);
    R30um = kron(f30m,eye(N))'*Ru*kron(f30m,eye(N));
    bdfb = exp(-1i*2*pi*omegadopplerbank(m)*(0:M-1)).';
    gt = kron(bdfb,ta.*at);         % Desired Vector common for both methods.
    gt30m =  kron(f30m,eye(N))'*gt;
    w30m = R30um\gt30m;             % Calculate K*N X 1 Adaptive Weight for m-th Doppler Bin.
    w30 = kron(f30m,eye(N))*w30m;

    for n=1:Lfd
        bt = exp(1i*2*pi*omegad(n)*(0:M-1)).'; % Dummy Target Doppler Steering Vector
        vt = kron(bt,at);
        SINR30_mat(m,n) = abs(w30'*vt)^2/real(w30'*Ru*w30);
    end
end

% PRI-Staggered SINR Loss Computations:
SINR30 = max(abs(SINR30_mat));
LSINR30 = SINR30/SNRo;

4. Beamspace pre-Doppler. Displaced Beams with TSN Calculations

Kt = 3;            % sub-CPI length (fixed).
Ks = 3;            % 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;   2;  1];          % Temporal Binomial Taper.
bt = [1;  -1;  1];          % Doppler Steering Vector.
ta = ones(N,1);
wmdb30 = zeros(M*N,M1,length(Ks));

utab30 = kron(tb.*bt,GTSN'*(ta.*at));    % Desired response. Eq. 236.
W30 = zeros(M*Ks,M1);

% W matrix construction.
for p=0:M1-1
    Jp = [zeros(p,Kt); eye(Kt); zeros(M-Kt-p,Kt)]; % Selector Matrix Jp.
    Rupdb30 = kron(Jp,GTSN)'*Ru*kron(Jp,GTSN);
    wp30 = Rupdb30\utab30;
    Wp30 = reshape(wp30,Ks,Kt);
    aux2 = Wp30*Jp.';
    W30(:,p+1)  = aux2(:);
end

for m=1:M1
    wmdb30(:,m) = kron(eye(M),GTSN)*W30*F30(:,m);
end

% LSINR Computations for Displaced Beams.
SINR30_cube = zeros(M1,Lfd);
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);
        SINR30_cube(m,n1) = abs(wmdb30(:,m)'*vt)^2/real(wmdb30(:,m)'*Ru*wmdb30(:,m));
    end
end
SINR30 = max(abs(SINR30_cube));
LSINRpreTSN30 = SINR30/SNRo;

5. Beamspace post-Doppler: Displaced Filters with TSN Calculations

Kt = 3;              % sub-CPI length (fixed).
Ks = 3;              % sub-apperture length or number of beams used.
M1 = M - Kt +1;      % Number of possible sub-CPI's.
N1 = N - Ks +1;      % Number of possible sub-appertures.

% Create Doppler Filterbank Matrix F for Displaced Filter Beamspace post-Doppler STAP:
dopplerfilterbank = linspace(0,300,M+1);
omegadopplerbank = dopplerfilterbank/fr;
F0 = zeros(M1,M);
for m=1:M
    F0(:,m) = 1/sqrt(M)*exp(-1i*2*pi*omegadopplerbank(m)*(0:M1-1));
end

F30 = diag(chebwin(M1,30))*F0;
SINRTSNdf_mat = zeros(M,Lfd);

SINR Computations for Displaced Filter Beamspace post-Doppler STAP

% Solve a separate adaptive problem in each Doppler bin m:

for m=1:M
    Fm30 = toeplitz([F30(:,m); zeros(Kt-1,1)],[F30(1,m) zeros(1,Kt-1)]);
    TmTSN = kron(Fm30,GTSN);
    Rum30 = TmTSN'*Ru*TmTSN;
    bdfb = exp(-1i*2*pi*omegadopplerbank(m)*(0:M-1)).';
    gt = kron(bdfb,at);
    utmTSM = TmTSN'*gt;
    wmTSN = Rum30\utmTSM;
    wTSN = TmTSN*wmTSN;

    for n=1:Lfd
        bt = exp(1i*2*pi*omegad(n)*(0:M-1)).'; % Dummy Target Doppler Steering Vector
        vt = kron(bt,at);
        SINRTSNdf_mat(m,n) = abs(wTSN'*vt)^2/real(wTSN'*Ru*wTSN);
    end
end

SINRTSNdf = max(abs(SINRTSNdf_mat));
LSINRTSNdf = SINRTSNdf/SNRo;

Plot the SINR Losses

figure('NumberTitle', 'off','Name', ...
    'Figure 64(a) SINR loss for example system, β=1, no misalinment, no intrinsic clutter motion');
plot(fd,10*log10(LSINRopt),'k','LineWidth',1.5);   hold on;
plot(fd,10*log10(LSINRsub2),'m','LineWidth',1.5)
plot(fd,10*log10(LSINR30),'color',[0 0.7 0],'LineWidth',1.5)
plot(fd,10*log10(LSINRpreTSN30),'LineWidth',1.5)
plot(fd,10*log10(LSINRTSNdf),'r','LineWidth',1.5)

ylabel('SINR Loss (dB)');
xlabel('Target Doppler Frequency (Hz)');
ylim([-30 1]);
xlim([-5 305]);
legend('Optimum Fully Adaptive','Elem-Pre (54)','Elem-Post (54)','Beam-Pre TSN (9)','Beam-Post TSN (9)', 'Location','South')
grid on;