MATLAB Examples

# Figure 62. Beamspace post-Doppler in a clutter-plus-jamming environment. K = 9. Ks = Kt = 3.

## Two-Step-Nulling (TSN).

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

```clc; clear; close all; ```

## Radar System Operational Parameters:

```fo = 450*1e6; % Operating Frequency in Hz Pt = 200*1e3; % Peak Transmit Power 200 kW dc = 0.06; % Duty Factor or Duty Cycle 6% Gt = 22; % Transmit Gain in dB Gr = 10; % Column Receive Gain in dB B = 4*1e6; % Receiver Instantaneous Bandwidth in Hz Ls = 4; % System Losses in dB fr = 300; % PRF in Hz Tr = 1/fr; % PRI in sec. M = 18; % Number of Pulses per CPI. Tp = 200*1e-6; % Pulse Width in sec. %Antenna Array Parameter Value N = 18; % Number of Elements % Number of Elements z: 4 % Element Pattern: Cosine Gel = 4; % Element Gain in dB % Transmit Taper: Uniform be = -30; % Element Backlobe Level in db Nc = 360; % Number of clutter patches uniformly distributed in azimuth. c = 299792458; % Speed of Light in m/sec. lambda = c/fo; % Operating wavelength in meters. d = lambda/2; % Interelement Spacing % Azimuth angle in degrees: phi = -180:179; Lphi = length(phi); f = zeros(1,Lphi); AF = zeros(1,Lphi); % Array Factor pre-allocation. % Platform Parameters: beta = 1; % beta parameter > 1 signifies Doppler ambiguous clutter. ha = 9000; % in meters. % Noise Power Computations: k = 1.3806488*1e-23; % Boltzmann Constant in J/K. To = 270; % Standard room Temperature in K. F = 3; % Receiver Noise Figure in dB; Ts = To*10^(F/10); % Receiver Temperature in K. Nn = k*Ts; % Receiver Noise PSD in Watts/Hz. Pn = Nn*B; % Receiver Noise Power in Watts sigma2 = 1; % Mormalized Noise Power in Watts. % Clutter Patch Geometry computations: Rcik = 130000; % (clutter) range of interest in meters. dphi = 2*pi/Nc; % Azimuth angle increment in rad. dR = c/2/B; % Radar Range Resolution in meters. Re = 6370000; % Earth Radius in meters. ae = 4/3*Re; % Effective Earth Radius in meters. % psi = -asin((Rc*Rc - ha*(ha+ 2*ae))/(2*Rc*ae)); % Grazing angle at the clutter patch in rad % (spherical earth model). psi = asin(ha/Rcik); % Grazing angle at the clutter patch in rad (flat earth model). theta = psi; % Elevation (look-down angle). Flat earth assumption. gamma = 10^(-3/10); % Terrain-dependent reflectivity factor. bwdth = 2*asin(1.4*lambda/N/d/pi); % Antenna Azimuth 3-dB beamwidth in rad. phia = 0; % Velocity Misalignment angle in degrees. ```

## Clutter Covariance Matrix Computations:

```% Calculate the Voltage Element Pattern: for i =1:Lphi if abs(phi(i))<=90 f(i) = cos(phi(i)*pi/180); else f(i) = 10^(be/10)*cos(phi(i)*pi/180); end end % Calculate the Array Factor (AF) (Voltage): steering_angle = 0; % Angle of beam steering in degrees. for k=1:Lphi AF(k) = sum(exp(-1i*2*pi/lambda*d*(0:N-1)*(sin(phi(k)*pi/180) ... - sin(steering_angle*pi/180))).*cos(phi(k)*pi/180)); end % Calculate the Full Array Transmit Power Gain: Gtgain = 10^(Gt/10)*abs(AF).^2; % Calculate the Element Receive Power Gain: grgain = 10^(Gel/10)*abs(f).^2; % Clutter Patch RCS Calculation: PatchArea = Rcik*dphi*dR*sec(psi); sigma0 = gamma*sin(psi); sigma = sigma0*PatchArea; % Calculate the Clutter to Noise Ration (CNR) for each clutter patch: ksi = Pt*Gtgain.*grgain*10^(Gr/10)*lambda^2*sigma/((4*pi)^3*Pn*10^(Ls/10)*Rcik^4); Ksic = sigma2*diag(ksi); % Platform Velocity for beta parameter value: va = round(beta*d*fr/2); Ita = d/lambda*cos(theta); % Calculate Spatial and Doppler Frequencies. % Spatial frequency of the k-th clutter patch. fspc = 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: ac = zeros(N,Nc); bc = zeros(M,Nc); Vc = zeros(M*N,Nc); for k=1:Nc ac(:,k) = exp(-1i*2*pi*fspc(k)*(0:N-1)); % Clutter Spatial Steering Vector. bc(:,k) = exp(-1i*2*pi*omegac(k)*(0:M-1)); % Clutter Temporal Steering Vector. Vc(:,k) = kron(bc(:,k),ac(:,k)); % Clutter Space-Time Steering Vector. end Rc = real(Vc*Ksic*Vc'); Rn = sigma2*eye(M*N); ```

## Jammer Covariance Matrix Calculation.

``` % Spatial frequency of the j-th jammer. J = 2; % Number of Jammers. thetaj = 0; phij = [-40 25]; % Jammer elevation and azimuth angles in degrees. R_j = [370 370]*1e3; Sj = 1e-3; % Jammer EIRP in Watts/Hz. fspj = d/lambda*cos(thetaj*pi/180)*sin(phij*pi/180); Lr = 2.0; % System Losses on Receive in dB. Aj = zeros(N,J); for j=1:J Aj(:,j) = exp(1i*2*pi*fspj(j)*(0:N-1)); % j-th Jammer Spatial Steering Vector. end indices = zeros(1,J); for j=1:J indices(j) = find(phi == phij(j)); end grgn = grgain(indices); Jo = (Sj*grgn*lambda^2)./((4*pi)^2.*10^(Lr/10).*R_j.^2); ksi_j = Jo./Nn; Ksi_j = sigma2*diag(ksi_j); % Spatial Jamming Covariance Matrix Phi_j = Aj*Ksi_j*Aj'; Rj = kron(eye(M),Phi_j); ```

## Analytic Interference Covariance Matrix for Fully Adaptive STAP:

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

## Spatial jammer-plus-noise Covariance Matrix:

``` Phi_jn = Phi_j + sigma2*eye(N); ```

## 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 = 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; ```

## with and without Two Step Nulling:

```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); for p=0:Ks-1 Jp = [zeros(p,N1); eye(N1); zeros(N-N1-p,N1)]; % Beam Selector Matrix Jsm. GTSN(p+1:N1+p,p+1) = inv(Jp.'*Phi_jn*Jp)*gtap; end GNoTSN = toeplitz([gtap; zeros(Ks-1,1);], [gtap(1) zeros(1,Ks-1)]); % N1 x Ks Beamformer matrix G. SINRNoTSNdf_mat = zeros(M,Lfd); SINRTSNdf_mat = zeros(M,Lfd); % Solve a separate adaptive problem in each Doppler bin m: for m=1:M ```

## SINR Computations for Displaced Filter Beamspace post-Doppler STAP:

``` Fm0 = toeplitz([F0(:,m); zeros(Kt-1,1)],[F0(1,m) zeros(1,Kt-1)]); Fm30 = toeplitz([F30(:,m); zeros(Kt-1,1)],[F30(1,m) zeros(1,Kt-1)]); TmNoTSN = kron(Fm0,GNoTSN); TmTSN = kron(Fm30,GTSN); Rum0 = TmNoTSN'*Ru*TmNoTSN; Rum30 = TmTSN'*Ru*TmTSN; bdfb = exp(-1i*2*pi*omegadopplerbank(m)*(0:M-1)).'; gt = kron(bdfb,at); utmNoTSM = TmNoTSN'*gt; utmTSM = TmTSN'*gt; wmNoTSN = Rum0\utmNoTSM; wmTSN = Rum30\utmTSM; wNoTSN = TmNoTSN*wmNoTSN; 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); SINRNoTSNdf_mat(m,n) = abs(wNoTSN'*vt)^2/real(wNoTSN'*Ru*wNoTSN); SINRTSNdf_mat(m,n) = abs(wTSN'*vt)^2/real(wTSN'*Ru*wTSN); end ```
``` end SINRNoTSNdf = max(abs(SINRNoTSNdf_mat)); SINRTSNdf = max(abs(SINRTSNdf_mat)); LSINRNoTSNdf = SINRNoTSNdf/SNRo; LSINRTSNdf = SINRTSNdf/SNRo; ```

## Plot the SINR Losses:

``` figure('NumberTitle', 'off','Name', ... 'Figure 62. Beamspace post-Doppler in a clutter-plus-jamming environment. K=9, Ksm = Ktm = 3.'); % Displaced Filter Beamspace post-Doppler STAP Results: plot(fd,10*log10(LSINRopt),'k','LineWidth',1.5) hold on; plot(fd,10*log10(LSINRNoTSNdf),'LineWidth',1.5) plot(fd,10*log10(LSINRTSNdf),'g','LineWidth',1.5) grid on; ylabel('SINR Loss (dB)'); xlabel('Target Doppler Frequency (Hz)'); title('Ks = Kt = 3'); ylim([-30 1]); xlim([-5 305]); legend('Optimum','No TSN','With TSN','Location','South'); ```