用相干信号子空间法写一段基于主瓣干扰抑制的阻塞矩阵预处理的宽带波束形成算法的MATLAB代码

% 宽带波束形成算法示例代码
% 使用相干信号子空间法及阻塞矩阵预处理抑制主瓣干扰

clc; clear all; close all;

%% 参数设置
K = 4; % 子带数
N = 8; % 天线数
M = 16; % 信号数
theta = [-30 0 30]; % 信号入射角度
f = [0.1 0.3 0.5]; % 子带频率
SNR = 20; % 信噪比

%% 信号发射矩阵
S = exp(1j*2*pi*(0:M-1)'*sind(theta)); % 信号发射矩阵

%% 接收信号模拟
noise = sqrt(0.5/SNR)*(randn(N,K,M)+1j*randn(N,K,M)); % 噪声信号
A = exp(-1j*2*pi*f'*(0:N-1)); % 天线阵列响应矩阵
X = zeros(N,K,M); % 接收信号矩阵
for k = 1:K
X(:,k,:) = A*S(:,k,:) + noise(:,k,:); % 接收信号
end

%% 相干信号子空间法波束形成
Rxx = zeros(N,N,K); % 信号协方差矩阵
for k = 1:K
Rxx(:,:,k) = X(:,k,:)*X(:,k,:)'/M;
end
Rx = sum(Rxx,3)/K; % 平均信号协方差矩阵
[V,D] = eig(Rx); % 特征分解
[lambda,idx] = sort(diag(D),'descend'); % 特征值排序
V = V(:,idx); % 特征向量排序
V1 = V(:,1:M-K); % 干扰信号子空间
V2 = V(:,M-K+1:end); % 信号子空间
W = V2*V2'/V1/V1'; % 波束权向量
P = W*diag(diag(W)'>0); % 阻塞矩阵预处理
y = zeros(K,M); % 波束输出
for k = 1:K
y(k,:) = P*X(:,k,:)*W; % 波束输出
end

%% 结果显示
figure;
subplot(121);
plot(-90:90, 20*log10(abs(W'*A*S(:,1))), '-r', ...
-90:90, 20*log10(abs(W'*A*S(:,2))), '-g', ...
-90:90, 20*log10(abs(W'*A*S(:,3))), '-b');
legend('S1', 'S2', 'S3');
xlabel('Angle (degree)'); ylabel('Amplitude (dB)');
title('Antenna array pattern without blocking matrix');
grid on;
subplot(122);
plot(-90:90, 20*log10(abs(P*W'*A*S(:,1))), '-r', ...
-90:90, 20*log10(abs(P*W'*A*S(:,2))), '-g', ...
-90:90, 20*log10(abs(P*W'*A*S(:,3))), '-b');
legend('S1', 'S2', 'S3');
xlabel('Angle (degree)'); ylabel('Amplitude (dB)');
title('Antenna array pattern with blocking matrix');
grid on;