MATLAB Examples

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

## Contents

clc; clear; close all;

## 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 = 2;   % beta parameter.
phia = 0;   % Velocity Misalignment Angle.

Rc = clutt_cov(ksi,beta);

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);
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 Beams 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);
F0 = zeros(M1,M1);
for m1=1:M1
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);
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 Beam pre-Doppler Calculations

selbeamdist = 6.35;                % Selected Beam 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. (a) Displaced-beam pre-Doppler STAP performance with doppler-ambiguous clutter (β = ',num2str(beta),') and no jamming.'],...
'Position',[1 1 1000 800]);

% 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');

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');