Figure 66.(a) SINR loss for example system, β=2.6, σ_υ = 0.2 m/sec, no misalignment.

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

K = 3;          % Number of Pulses per sub-CPI.
M1 = M - K + 1; % Number of sub-CPI's for K2=3.
beta = 2.6;      % beta parameter.
phia = 0;        % Velocity Misalignment Angle.
sigma_icm = 0.2; % Intrinsic Clutter Motion Standard Deviation in meters/sec.
kc = 2*sigma_icm/lambda;
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(K,Nc);
Rc = zeros(M*N,M*N);
Rcsub2 = zeros(K*N,K*N);

gammac =  exp(-(kc^2*Tr^2*(0:M-1).^2)/2);
gammac_sub = exp(-(kc^2*Tr^2*(0:K-1).^2)/2);
Gammac = toeplitz(gammac);
Gammac_sub = toeplitz(gammac_sub);

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:K-1));  % Time Steering Vector for K = 3
    Rc = Rc + ksi(k)*kron(Gammac.*(b(:,k)*b(:,k)'), a(:,k)*a(:,k)');
    Rcsub2 = Rcsub2 + ksi(k)*kron(Gammac_sub.*(b2(:,k)*b2(:,k)'), a(:,k)*a(:,k)');
end

Thermal Noise Covariance Matrices

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

Jammer Covariance Matrix Calculation

jamm_cov;

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

Spatial jammer-plus-noise Covariance Matrix

Phi_jn = Phi_j + sigma2*eye(N);

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.

1. Optimum Fully Adaptive STAP

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

InvRu = inv(Ru);
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. Element-Space pre-Doppler, sub-CPI with K=3.

bt2 = [1; -1; 1;];
tb2 = [1;  2; 1;];             % Binomial Doppler Taper for K = 3.

gt2 = kron(tb2.*bt2,ta.*at);
wsub2 = Rusub2\gt2;            % Weight Vector for K = 3;

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

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

for m=1:length(dopplerfilterbank2)
    um = exp(1i*2*pi*omegadopplerbank2(m)*(0:M1-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;
P = floor(K/2);

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 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 Two Step Nulling

Kt = 3;          % sub-CPI length (fixed).
Ks = 5;          % 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));

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.
    W30 = zeros(M*Ks(n),M1);
    GTSN = zeros(N,Ks(n));

    for p=0:Ks(n)-1
        Jp = [zeros(p,N1); eye(N1); zeros(N-N1-p,N1)]; % Beam Selector Matrix Jsm.
        GTSN(p+1:N1+p,p+1) = (Jp.'*Phi_jn*Jp)\g30;
    end

    utab30 = kron(tb.*bt,GTSN'*(ta.*at));

    % 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(n),Kt);
        aux2 = Wp30*Jp.';
        W30(:,p+1)  = aux2(:);
    end

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

end

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

LSINRpreTSN   = SINR30/SNRo;

5. Beamspace post-Doppler: Displaced Filters with Two Step Nulling

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

% 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;
wTSN = zeros(M*N,length(Ks));

for n=1:length(Ks)

    N1 = N - Ks(n) +1;             % Number of possible sub-appertures.

    % Beamformer G Matrix Construction for Displaced Filter Beamspace post-Doppler STAP
    % with Two Step Nulling (TSN):
    taper_level = 30;                     % Tapering level in db.
    g0 = exp(-1i*2*pi*fspt*(0:N1-1)).';
    gtap = chebwin(N1,taper_level).*g0;   % spatial taper application.

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

    SINRTSNdf_cube = zeros(M1,Lfd,length(Ks));

    % 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(:,n) = TmTSN*wmTSN;

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

    end

    SINR30(n,:) = max(abs(SINRTSNdf_cube(:,:,n)));

end

LSINRpostTSN = SINR30/SNRo;

Plot the SINR Losses

figure('NumberTitle', 'off','Name', ...
       'Figure 66(a) SINR loss for example system, β=2.6, no misalinment, σ_υ = 0.2 m/sec');
plot(fd,10*log10(LSINRopt),'LineWidth',1.5); hold on;
plot(fd,10*log10(LSINRsub2),'r','LineWidth',1.5);
plot(fd,10*log10(LSINR30),'g','LineWidth',1.5);
plot(fd,10*log10(LSINRpreTSN(1,:)),'k','LineWidth',1.5);
plot(fd,10*log10(LSINRpostTSN(1,:)),'m','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 (15)','Beam-Post TSN (15)','Location','South')
grid on;