MATLAB Examples

Figure 60. Beamspace post-Doppler in a clutter-only environment. K = 16. Ks = Kt = 4.

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 = 1;   % 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.
bt = exp(-1i*2*pi*fdt*(0:M-1)).';  % Target Doppler Steering Vector.
ta = chebwin(N,30);
tb = chebwin(M,30);
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

Displaced Beam pre-Doppler Calculations

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

F40 = diag(chebwin(M1,40))*F0;

Beamformer G Matrix Construction for Displaced Filter Beamspace post-Doppler STAP

g0  = exp(-1i*2*pi*fspt*(0:N1-1)).';
g30 = chebwin(N1,30).*g0;            % 30-dB spatial taper.
G0  = toeplitz([g0;  zeros(Ks-1,1);],  [g0(1) zeros(1,Ks-1)]);  % N1 x Ks Beamformer matrix G.
G30 = toeplitz([g30; zeros(Ks-1,1);], [g30(1) zeros(1,Ks-1)]);  % N1 x Ks Beamformer matrix G.

SINR0df_mat  = zeros(M,Lfd);
SINR30df_mat = zeros(M,Lfd);
SINR0af_mat  = zeros(M,Lfd);
SINR30af_mat = zeros(M,Lfd);

Create Doppler Filterbank Matrix F for Adjacent Filter Beamspace post-Doppler STAP

U2 = zeros(M,M);
Pt = floor(Kt/2);
if  mod(Kt,2)   % If Kt is odd
    for m=1:M
        U2(:,m) = 1/sqrt(M)*exp(-1i*2*pi*omegadopplerbank(m)*(0:M-1));
    end
else            % if Kt is even:
    outomegadoppler = zeros(1,M);
    for m=1:M
        outomegadoppler(m) = (omegadopplerbank(m) + omegadopplerbank(m+1))/2;
        U2(:,m) = 1/sqrt(M)*exp(-1i*2*pi*outomegadoppler(m)*(0:M-1));
    end
end

td0af  = ones(M,1);              % Uniform Doppler Taper.
td40af = chebwin(M,40);          % 40-dB Chebyshev Doppler Taper.

% Create Doppler Filter Bank in Fab matrix for Adjacent Filter Beamspace post-Doppler method:
Faf0  = diag(td0af)*U2;
Faf40 = diag(td40af)*U2;

Fmaf0  = zeros(M,Kt,M);
Fmaf40 = zeros(M,Kt,M);

for m=1:M
    if mod(Kt,2)                  % if K is odd and >1.
        if (m-Pt>0)&&(m+Pt<=M)
            Fmaf0(:,:,m)  = Faf0(:,m-Pt:m+Pt);                       % Eq. 231.
            Fmaf40(:,:,m) = Faf40(:,m-Pt:m+Pt);
        elseif (m-Pt<=0)&&(m+Pt<=M)
            Fmaf0(:,:,m)  = [Faf0(:,M+(m-Pt):M)   Faf0(:,1:m+Pt)];   % Eq. 231.
            Fmaf40(:,:,m) = [Faf40(:,M+(m-Pt):M) Faf40(:,1:m+Pt)];
        elseif m+Pt>M
            Fmaf0(:,:,m)  = [Faf0(:,m-Pt:M)   Faf0(:,1:m+Pt-M)];     % Eq. 231.
            Fmaf40(:,:,m) = [Faf40(:,m-Pt:M) Faf40(:,1:m+Pt-M)];
        end

    else      % if K is even.

        if (m-Pt>0)&&(m+Pt<=M+1)
            Fmaf0(:,:,m)  = Faf0(:,m-Pt:m+Pt-1);                      % Eq. 231.
            Fmaf40(:,:,m) = Faf40(:,m-Pt:m+Pt-1);
        elseif (m-Pt<=0)&&(m+Pt<=M)
            Fmaf0(:,:,m)  = [Faf0(:,M+(m-Pt):M)   Faf0(:,1:m+Pt-1)];  % Eq. 231.
            Fmaf40(:,:,m) = [Faf40(:,M+(m-Pt):M) Faf40(:,1:m+Pt-1)];
        elseif m+Pt>M+1
            Fmaf0(:,:,m)  = [Faf0(:,m-Pt:M)     Faf0(:,1:m-M+Pt-1)];  % Eq. 231.
            Fmaf40(:,:,m) = [Faf40(:,m-Pt:M)   Faf40(:,1:m-M+Pt-1)];
        end
    end

end

Create Beamformer Matrix G for Adjacent Filter Beamspace post-Doppler STAP

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

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

Create Beam Selector Matrix Jsm for Adjacent Filter Beamspace post-Doppler STAP

Ps = floor(Ks/2);
if mod(Ks,2)  % if Ks is odd
    Jsm = [zeros(N/2-Ps,Ks); eye(Ks); zeros(N/2-Ps-1,Ks)]; % Beam Selector Matrix Jsm.
    G0af = Godd;                                           % N x Ks Beamformer matrix.
    G30af = diag(chebwin(N,30))*G0af;
else
    Jsm = [zeros(N/2-Ps,Ks); eye(Ks); zeros(N/2-Ps,Ks)];   % Beam Selector Matrix Jsm.
    G0af = Geven;                                          % N x Ks Beamformer matrix.
    G30af = diag(chebwin(N,30))*G0af;
end

G0maf = G0af*Jsm;
G30maf = G30af*Jsm;

A. SINR Computations for Displaced Filter Beamspace post-Doppler STAP

for m=1:M

    Fm0  = toeplitz([F0(:,m);  zeros(Kt-1,1)],[F0(1,m)  zeros(1,Kt-1)]);
    Fm30 = toeplitz([F40(:,m); zeros(Kt-1,1)],[F40(1,m) zeros(1,Kt-1)]);

    Tm0  = kron(Fm0,G0);
    Tm30 = kron(Fm30,G30);

    Rum0  = Tm0'*Ru*Tm0;
    Rum30 = Tm30'*Ru*Tm30;

    bdfb = exp(-1i*2*pi*omegadopplerbank(m)*(0:M-1)).';
    gt = kron(bdfb,at);

    utm0  = Tm0'*gt;
    utm30 = Tm30'*gt;

    wm0  = Rum0\utm0;
    wm30 = Rum30\utm30;

    w0  = Tm0*wm0;
    w30 = Tm30*wm30;

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

SINR0df  = max(abs(SINR0df_mat));
SINR30df = max(abs(SINR30df_mat));

LSINR0df  = SINR0df/SNRo;
LSINR30df = SINR30df/SNRo;

B. SINR Computations for Adjacent Filter Beamspace post-Doppler STAP

for m=1:M

    F0maf  = Fmaf0(:,:,m);
    F30maf = Fmaf40(:,:,m);

    Tm0af = kron(F0maf,G0maf);
    Tm30af = kron(F30maf,G30maf);

    Rum0af  = Tm0af'*Ru*Tm0af;
    Rum30af = Tm30af'*Ru*Tm30af;

    bdfb = exp(-1i*2*pi*omegadopplerbank(m)*(0:M-1)).';
    gt = kron(bdfb,at);

    utm0af  = Tm0af'*gt;
    utm30af = Tm30af'*gt;

    wm0af  = Rum0af\utm0af;
    wm30af = Rum30af\utm30af;

    w0af  = Tm0af*wm0af;
    w30af = Tm30af*wm30af;

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

end

SINR0af  = max(abs(SINR0af_mat));
SINR30af = max(abs(SINR30af_mat));

LSINR0af  = SINR0af/SNRo;
LSINR30af = SINR30af/SNRo;

Plot the SINR Losses

figure('NumberTitle', 'off','Name', ...
       'Figure 60. Beamspace post-Doppler in a clutter only environment. K=16, Ksm = Ktm = 4.');

%  Displaced Filter Beamspace post-Doppler STAP Results:
plot(fd,10*log10(LSINRopt),'k','LineWidth',1.5)
hold on;
plot(fd,10*log10(LSINR0df),'LineWidth',1.5)
plot(fd,10*log10(LSINR30df),'g','LineWidth',1.5)
grid on;
ylabel('SINR Loss (dB)');
xlabel('Target Doppler Frequency (Hz)');
ylim([-30 1]);
xlim([-5 305]);

%  Adjacent Filter Beamspace post-Doppler STAP Results:
plot(fd,10*log10(LSINR0af),'r','LineWidth',1.5)
hold on;
plot(fd,10*log10(LSINR30af),'m','LineWidth',1.5)
grid on;
ylabel('SINR Loss (dB)');
xlabel('Target Doppler Frequency (Hz)');
ylim([-30 1]);
xlim([-5 305]);
legend('Optimum','Displaced-Uniform', 'Displaced-Tapered', 'Adjacent-Uniform', 'Adjacent-Tapered','Location','South');