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.

Contents

Workspace Initialization.

clc; clear; close all;

Signal Definitions.

M   = 40;        % Number of Array Elements.
K1 = 100;      % Number of Signal Snapshots Used.
K2 = 500;      % Number of Signal Snapshots Used.
K3 = 1000;    % Number of Signal Snapshots Used.
K4 = 5000;    % Number of Signal Snapshots Used.

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^0.5); % Diagonal Loading level.

u_s = (d/lambda)*sin(phi_s*pi/180);  % Normalized Spatial Frequency.
u_i  = (d/lambda)*sin(phi_i*pi/180);   % Normalized Spatial Frequency for Source of Interference.

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);   % Interference Steering Vector.

% Calculation of the true i+n 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);   % We better use:
copt  = InvRv_s/(v_s'*InvRv_s);

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

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

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

for k=1:Experiment_Repetitions

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

Ave_SMI_Beam_Pat1 = 1/Experiment_Repetitions*sum(SMI_Beam_Pat1);
Ave_SMI_Beam_Pat2 = 1/Experiment_Repetitions*sum(SMI_Beam_Pat2);
Ave_SMI_Beam_Pat3 = 1/Experiment_Repetitions*sum(SMI_Beam_Pat3);
Ave_SMI_Beam_Pat4 = 1/Experiment_Repetitions*sum(SMI_Beam_Pat4);

Ave_SMI_Beam_Pat1_loaded = 1/Experiment_Repetitions*sum(SMI_Beam_Pat1_loaded);
Ave_SMI_Beam_Pat2_loaded = 1/Experiment_Repetitions*sum(SMI_Beam_Pat2_loaded);
Ave_SMI_Beam_Pat3_loaded = 1/Experiment_Repetitions*sum(SMI_Beam_Pat3_loaded);
Ave_SMI_Beam_Pat4_loaded = 1/Experiment_Repetitions*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');
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');
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');
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');
xlim([-10 90]);
ylim([-50 5]);
grid on;
xlabel('Angle (deg)');
ylabel('Power (dB)');

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 eigenvalues of all three kinds of correlation matrices:

figure('NumberTitle', 'off','Name','Figure 11.25','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 Eigenvalues');
legend('True','Sample for K=100','Sample Loaded','Location','Best');
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 Eigenvalues');
legend('True','Sample for K=500','Sample Loaded','Location','Best');
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 Eigenvalues');
legend('True','Sample for K=1000','Sample Loaded','Location','Best');
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 Eigenvalues');
legend('True','Sample for K=5000','Sample Loaded','Location','Best');
xlabel('eigenvalue Number');
ylabel('Eigenvalue Power (dB)','Rotation',90);
grid on;