MATLAB Examples

Example 11.5.3. - Diagonal Loading of the Sample Correlation Matrix.

In this example, we explore the use of diagonal loading of the sample correlation matrix to control the sidelobe levels of the SMI adaptive beamformer using the same set of parameters as in Example 11.5.2. The beampatterns for the SMI adaptive beamformer and the diagonally loaded SMI adaptive beamformer are shown in Figure 11.24 along with the beampattern for the optimum beamformer for $-10< \phi < 90$. The diagonal loading level was set to 5 dB above the thermal noise power, that is, $\sigma_{I}^2 = 10^0.5$, and the sample support was K = 100.

Copyright 2016-2026, Ilias S. Konsoulas.

Contents

Workspace Initialization.

clc; clear; close all;

Signal Definitions.

M  = 40;       % Number of Array Elements.
K1 = 100;      % Number of Signal Snapshots Used for Test Case #1.
K2 = 500;      % Number of Signal Snapshots Used for Test Case #2.
K3 = 1000;     % Number of Signal Snapshots Used for Test Case #3.
K4 = 5000;     % Number of Signal Snapshots Used for Test Case #4.

lambda = 1;                                % Incoming Signal Wavelength in (m).
d      = lambda/2;                         % Interelement Distance in (m).
SNR    = 20;                               % target SNR in dBs
phi_s  = 0;                                % target azimuth angle in degrees
phi_i  = 10;                               % interference angle in degrees
INR    = 70;                               % Interference #1 INR in dBs
sigma_load = sqrt(10^(5/10));              % Diagonal Loading level of 5dB.

u_s = (d/lambda)*sin(phi_s*pi/180);        % Target Normalized Spatial Frequency.
u_i = (d/lambda)*sin(phi_i*pi/180);        % Interferer Normalized Spatial Frequency.

v_s = exp(-1i*2*pi*u_s*(0:M-1).')/sqrt(M); % Target Steering Vector.
v_i = exp(-1i*2*pi*u_i*(0:M-1).')/sqrt(M); % Interferer Steering Vector.

Calculation of the true Interference-plus-Noise Autocorrelation matrix:

R_ipn = 10^(INR/10)*(v_i*v_i') + eye(M);

% InvR = inv(R_ipn); However, the following lines are faster:
InvRv_s = R_ipn\v_s;

Optimum Beamformer computed for a phi_s = 0 deg.

Instead of: copt = InvR*v_s/(v_s'*InvR*v_s); More efficiently we use:

copt  = InvRv_s/(v_s'*InvRv_s);

Sample Matrix Inversion Monte Carlo Runs.

Run the experiment for 100 times (Monte-Carlo runs) and average over the results to get the correct curves of Fig. 11.24.

csmi1 = zeros(M,1); csmi2 = zeros(M,1);
csmi3 = zeros(M,1); csmi4 = zeros(M,1);
csmi1_loaded = zeros(M,1); csmi2_loaded = zeros(M,1);
csmi3_loaded = zeros(M,1); csmi4_loaded = zeros(M,1);

MC_Runs = 100; NSamples = 1e3;
angle = -10:100/NSamples:90-60/NSamples;

SMI_Beam_Pat1        = zeros(MC_Runs,NSamples);
SMI_Beam_Pat2        = zeros(MC_Runs,NSamples);
SMI_Beam_Pat3        = zeros(MC_Runs,NSamples);
SMI_Beam_Pat4        = zeros(MC_Runs,NSamples);
SMI_Beam_Pat1_loaded = zeros(MC_Runs,NSamples);
SMI_Beam_Pat2_loaded = zeros(MC_Runs,NSamples);
SMI_Beam_Pat3_loaded = zeros(MC_Runs,NSamples);
SMI_Beam_Pat4_loaded = zeros(MC_Runs,NSamples);

for k=1:MC_Runs

    % Uncorrelated unit power thermal noise samples drawn from a complex Gaussian distribution
    w1 = (randn(M,K1) + 1i*randn(M,K1))/sqrt(2);
    w2 = (randn(M,K2) + 1i*randn(M,K2))/sqrt(2);
    w3 = (randn(M,K3) + 1i*randn(M,K3))/sqrt(2);
    w4 = (randn(M,K4) + 1i*randn(M,K4))/sqrt(2);

    % The interference (jammer) vector are generated by:
    i_x1=(10^(INR/20))*v_i*(randn(1,K1) + 1i*randn(1,K1))/sqrt(2);
    i_x2=(10^(INR/20))*v_i*(randn(1,K2) + 1i*randn(1,K2))/sqrt(2);
    i_x3=(10^(INR/20))*v_i*(randn(1,K3) + 1i*randn(1,K3))/sqrt(2);
    i_x4=(10^(INR/20))*v_i*(randn(1,K4) + 1i*randn(1,K4))/sqrt(2);

    % The interference-plus-noise signal x_i+n is generated by:
    iplusn1 = i_x1 + w1;
    iplusn2 = i_x2 + w2;
    iplusn3 = i_x3 + w3;
    iplusn4 = i_x4 + w4;

    % Compute the sample interference + noise correlation matrices:
    R1 = 1/K1*(iplusn1*iplusn1');
    R2 = 1/K2*(iplusn2*iplusn2');
    R3 = 1/K3*(iplusn3*iplusn3');
    R4 = 1/K4*(iplusn4*iplusn4');

    % Apply diagonal loading to the sample interference + noise correlation matrices:
    R1_loaded = R1 + sigma_load^2*eye(M,M);
    R2_loaded = R2 + sigma_load^2*eye(M,M);
    R3_loaded = R3 + sigma_load^2*eye(M,M);
    R4_loaded = R4 + sigma_load^2*eye(M,M);

    InvR1v_s = R1\v_s;
    InvR2v_s = R2\v_s;
    InvR3v_s = R3\v_s;
    InvR4v_s = R4\v_s;

    csmi1 = InvR1v_s/(v_s'*InvR1v_s);  % SMI Beamformer computed for a phi_s = 0 deg and K1 = 100 snapshots.
    csmi2 = InvR2v_s/(v_s'*InvR2v_s);  % SMI Beamformer computed for a phi_s = 0 deg and K2 = 500 snapshots.
    csmi3 = InvR3v_s/(v_s'*InvR3v_s);  % SMI Beamformer computed for a phi_s = 0 deg and K3 = 1000 snapshots.
    csmi4 = InvR4v_s/(v_s'*InvR4v_s);  % SMI Beamformer computed for a phi_s = 0 deg and K4 = 5000 snapshots.

    InvR1_loadedv_s = R1_loaded\v_s;
    InvR2_loadedv_s = R2_loaded\v_s;
    InvR3_loadedv_s = R3_loaded\v_s;
    InvR4_loadedv_s = R4_loaded\v_s;

    csmi1_loaded = InvR1_loadedv_s/(v_s'*InvR1_loadedv_s);
    csmi2_loaded = InvR2_loadedv_s/(v_s'*InvR2_loadedv_s);
    csmi3_loaded = InvR3_loadedv_s/(v_s'*InvR3_loadedv_s);
    csmi4_loaded = InvR4_loadedv_s/(v_s'*InvR4_loadedv_s);

    for l=1:NSamples
         u = (d/lambda)*sin(angle(l)*pi/180);
         v = exp(-1i*2*pi*u*(0:M-1).')/sqrt(M); % Azimuth Scanning Steering Vector.
         SMI_Beam_Pat1(k,l) = csmi1'*v;
         SMI_Beam_Pat2(k,l) = csmi2'*v;
         SMI_Beam_Pat3(k,l) = csmi3'*v;
         SMI_Beam_Pat4(k,l) = csmi4'*v;
         SMI_Beam_Pat1_loaded(k,l) = csmi1_loaded'*v;
         SMI_Beam_Pat2_loaded(k,l) = csmi2_loaded'*v;
         SMI_Beam_Pat3_loaded(k,l) = csmi3_loaded'*v;
         SMI_Beam_Pat4_loaded(k,l) = csmi4_loaded'*v;
    end

    SMI_Beam_Pat1(k,:) = 10*log10(abs(SMI_Beam_Pat1(k,:)).^2);
    SMI_Beam_Pat2(k,:) = 10*log10(abs(SMI_Beam_Pat2(k,:)).^2);
    SMI_Beam_Pat3(k,:) = 10*log10(abs(SMI_Beam_Pat3(k,:)).^2);
    SMI_Beam_Pat4(k,:) = 10*log10(abs(SMI_Beam_Pat4(k,:)).^2);
    SMI_Beam_Pat1_loaded(k,:) = 10*log10(abs(SMI_Beam_Pat1_loaded(k,:)).^2);
    SMI_Beam_Pat2_loaded(k,:) = 10*log10(abs(SMI_Beam_Pat2_loaded(k,:)).^2);
    SMI_Beam_Pat3_loaded(k,:) = 10*log10(abs(SMI_Beam_Pat3_loaded(k,:)).^2);
    SMI_Beam_Pat4_loaded(k,:) = 10*log10(abs(SMI_Beam_Pat4_loaded(k,:)).^2);

end

% Average the caclulated beampatterns:
Ave_SMI_Beam_Pat1 = 1/MC_Runs*sum(SMI_Beam_Pat1);
Ave_SMI_Beam_Pat2 = 1/MC_Runs*sum(SMI_Beam_Pat2);
Ave_SMI_Beam_Pat3 = 1/MC_Runs*sum(SMI_Beam_Pat3);
Ave_SMI_Beam_Pat4 = 1/MC_Runs*sum(SMI_Beam_Pat4);

Ave_SMI_Beam_Pat1_loaded = 1/MC_Runs*sum(SMI_Beam_Pat1_loaded);
Ave_SMI_Beam_Pat2_loaded = 1/MC_Runs*sum(SMI_Beam_Pat2_loaded);
Ave_SMI_Beam_Pat3_loaded = 1/MC_Runs*sum(SMI_Beam_Pat3_loaded);
Ave_SMI_Beam_Pat4_loaded = 1/MC_Runs*sum(SMI_Beam_Pat4_loaded);

Compute the Beampattern of the Optimal Beamformer:

Opt_Beam_Pat  = zeros(NSamples,1);
for l=1:NSamples
    u = (d/lambda)*sin(angle(l)*pi/180);
    v = exp(-1i*2*pi*u*(0:M-1).')/sqrt(M);
    Opt_Beam_Pat(l)  = copt'*v;
end

Plot the Average Beampatterns and compare with the Optimum Beampattern.

figure('NumberTitle', 'off','Name','Figure 11.24','Position',[150 100 900 600]);
subplot(2,2,1);
plot(angle,10*log10(abs(Opt_Beam_Pat).^2),'r',angle,Ave_SMI_Beam_Pat1,'k','LineWidth',1.5);
hold on;
plot(angle,Ave_SMI_Beam_Pat1_loaded,'LineWidth',1.5);
title('K = 100');
legend('Optimum BF BP','Avereged SMI BP','Averaged SMI Loaded BP','Location','NorthEast');
xlim([-10 90]);
ylim([-50 5]);
grid on;
xlabel('Angle (deg)');
ylabel('Power (dB)');

subplot(2,2,2);
plot(angle,10*log10(abs(Opt_Beam_Pat).^2),'r',angle,Ave_SMI_Beam_Pat2,'k','LineWidth',1.5)
hold on;
plot(angle,Ave_SMI_Beam_Pat2_loaded,'LineWidth',1.5);
title('K = 500');
legend('Optimum BF BP','Avereged SMI BP','Averaged SMI Loaded BP','Location','NorthEast');
xlim([-10 90]);
ylim([-50 5]);
grid on;
xlabel('Angle (deg)');
ylabel('Power (dB)');

subplot(2,2,3);
plot(angle,10*log10(abs(Opt_Beam_Pat).^2),'r',angle,Ave_SMI_Beam_Pat3,'k','LineWidth',1.5)
hold on;
plot(angle,Ave_SMI_Beam_Pat3_loaded,'LineWidth',1.5);
title('K = 1000');
legend('Optimum BF BP','Avereged SMI BP','Averaged SMI Loaded BP','Location','NorthEast');
xlim([-10 90]);
ylim([-50 5]);
grid on;
xlabel('Angle (deg)');
ylabel('Power (dB)');

subplot(2,2,4);
plot(angle,10*log10(abs(Opt_Beam_Pat).^2),'r',angle,Ave_SMI_Beam_Pat4,'k','LineWidth',1.5);
hold on;
plot(angle,Ave_SMI_Beam_Pat4_loaded,'LineWidth',1.5);
title('K = 5000');
legend('Optimum BF BP','Avereged SMI BP','Averaged SMI Loaded BP','Location','NorthEast');
xlim([-10 90]);
ylim([-50 5]);
grid on;
xlabel('Angle (deg)');
ylabel('Power (dB)');

tightfig;

Compute and Plot the Correlation Matrices eigenvalues:

eig_true    = sort(real(eig(R_ipn)),'descend');
eig_sample1 = sort(real(eig(R1)),   'descend');
eig_sample2 = sort(real(eig(R2)),   'descend');
eig_sample3 = sort(real(eig(R3)),   'descend');
eig_sample4 = sort(real(eig(R4)),   'descend');

eig_sample1_loaded = sort(real(eig(R1_loaded)),'descend');
eig_sample2_loaded = sort(real(eig(R2_loaded)),'descend');
eig_sample3_loaded = sort(real(eig(R3_loaded)),'descend');
eig_sample4_loaded = sort(real(eig(R4_loaded)),'descend');

Plot the Thermal Noise Eigenvalues of all three kinds of correlation matrices:

figure('NumberTitle', 'off','Name','Figure 11.25 - Noise Eigenvalues','Position',[150 100 900 600]);
subplot(2,2,1);
plot(2:M,10*log10(eig_true(2:end)),'b.');
hold on;
plot(2:M,10*log10(eig_sample1(2:end)),'r.');
plot(2:M,10*log10(eig_sample1_loaded(2:end)),'g.');
line([0 40],[5  5],'LineStyle','--','LineWidth',2,'Color','m');
title('Correlation Matrix Noise Eigenvalues');
legend('True','Sample for K=100','Sample Loaded','Location','SouthWest');
ylim([-10 8]);
xlabel('Eigenvalue Number');
ylabel('Eigenvalue Power (dB)','Rotation',90);
grid on;

subplot(2,2,2);
plot(2:M,10*log10(eig_true(2:end)),'b.');
hold on;
plot(2:M,10*log10(eig_sample2(2:end)),'r.');
plot(2:M,10*log10(eig_sample2_loaded(2:end)),'g.');
line([0 40],[5  5],'LineStyle','--','LineWidth',2,'Color','m');
title('Correlation Matrix Noise Eigenvalues');
legend('True','Sample for K=500','Sample Loaded','Location','SouthWest');
ylim([-10 8]);
xlabel('Eigenvalue Number');
ylabel('Eigenvalue Power (dB)','Rotation',90);
grid on;

subplot(2,2,3);
plot(2:M,10*log10(eig_true(2:end)),'b.');
hold on;
plot(2:M,10*log10(eig_sample3(2:end)),'r.');
plot(2:M,10*log10(eig_sample3_loaded(2:end)),'g.');
line([0 40],[5  5],'LineStyle','--','LineWidth',2,'Color','m');
title('Correlation Matrix Noise Eigenvalues');
legend('True','Sample for K=1000','Sample Loaded','Location','SouthWest');
ylim([-10 8]);
xlabel('Eigenvalue Number');
ylabel('Eigenvalue Power (dB)','Rotation',90);
grid on;

subplot(2,2,4);
plot(2:M,10*log10(eig_true(2:end)),'b.');
hold on;
plot(10*log10(eig_sample4(2:end)),'r.');
plot(10*log10(eig_sample4_loaded(2:end)),'g.');
line([0 40],[5  5],'LineStyle','--','LineWidth',2,'Color','m');
title('Correlation Matrix Noise Eigenvalues');
legend('True','Sample for K=5000','Sample Loaded','Location','SouthWest');
ylim([-10 8]);
xlabel('Eigenvalue Number');
ylabel('Eigenvalue Power (dB)','Rotation',90);
grid on;

tightfig;