Traditional Comm笔记【9】：阵列信号处理及MATLAB实现（第2版）阅读随笔(四)

Chapter6 阵列信号处理常用MATLAB代码

(1). LCMV波束形成MATLAB程序：

“线性约束最小方差准则（LCMV）”在对有用信号形式和信号来向完全未知情况下，在某种约束条件下使阵列输出的方差最小。

"LCMV"方法的代价函数可以表示为 $J(\omega )={\omega }^{H}R\omega$，约束条件为 ${\omega }^{H}a(\theta )=f$。取 $f=1$，得到最佳解为 $\omega =R^{-1}c[c^H{R}^{-1}c{]}^{-1}$。

clc; 
close all
clear all; 

M=18;                                     % 天线数
L=100;                                    % 快拍数
thetas=10;                                % 信号入射角度 
thetai=[-30 30];                          % 干扰入射角度 
n=[0:M-1]';                               % 构造一个一维列矩阵 

vs=exp(-1i*pi*n*sin(thetas/180*pi));      % 信号方向矢量 
vi=exp(-1i*pi*n*sin(thetai/180*pi));      % 干扰方向矢量 
f=16000;                                  % 载波频率
t=[0:1:L-1]/200;                          % 构造时间变量
snr=10;                                   % 信噪比 
inr=10;                                   % 干噪比
xs=sqrt(10^(snr/10))*vs*exp(1i*2*pi*f*t); % 构造有用信号
xi=sqrt(10^(inr/10)/2)*vi*[randn(length(thetai),L)+1i*randn(length(thetai),L)];
% 构造干扰信号
noise=[randn(M,L)+1i*randn(M,L)]/sqrt(2); % 噪声
X=xi+noise;                               % 含噪信号
R=X*X'/L;                                 % 构造协方差矩阵
wop1=inv(R)*vs/(vs'*inv(R)*vs);           % 波束形成
sita=48*[-1:0.001:1];                     % 扫描方向范围
v=exp(-1i*pi*n*sin(sita/180*pi));         % 扫描方向向量 
B=abs(wop1'*v);                           % 求不同角度的增益
plot(sita,20*log10(B/max(B)),'k'); 
title('波束图');xlabel('角度/degree');ylabel('波束图/dB'); 
grid on 
axis([-48 48 -50 0]); 
hold off

figure1 LCMV波束形成算法仿真结果

(2). LMS自适应波束形成MATLAB程序：

“最小均方算法（LMS）”采用迭代模式，在每个迭代步骤 $n$ 时刻的权向量加上一个校正量后，即组成 $n+1$ 时刻的权向量，用它逼近最佳权向量。

clear all
close all
clc
M=16;                  % 天线数
K=2;                   % 信源数
theta=[0 30];          % 信号入射角度
d=0.3;                 % 天线间距
N=500;                 % 采样数
Meann=0;   varn=1;     % 噪声均值、方差
SNR=20;
INR=20;
pp=zeros(100,N);pp1=zeros(100,N);
rvar1=sqrt(varn)*10^(SNR/20);        % 信号功率
rvar2=sqrt(varn)*10^(INR/20);        % 干扰功率
for q=1:100
    s=[rvar1*exp(1i*2*pi*(50*0.001*[0:N-1]));rvar2*exp(1i*2*pi*(100*0.001*[0:N-1]+rand))];  % 生成源信号
    A=exp(-1i*2*pi*d*[0:M-1].'*sin(theta*pi/180));   % 方向向量
    e=sqrt(varn/2)*(randn(M,N)+1i*randn(M,N));       % 噪声
    Y=A*s+e;                                         % 接收信号

    % LMS算法
    L=200;
    de=s(1,:);
    mu=0.0005;
    w=zeros(M,1);
    for k=1:N
        y(k)=w'*Y(:,k);                              % 预测下一个采样和误差
    e(k)=de(k)-y(k);                                 % 误差
    w=w+mu*Y(:,k)*conj(e(k));                        % 调整权向量
    end
end

% 波束形成
beam=zeros(1,L);
for i=1:L
    a=exp(-1i*2*pi*d*[0:M-1].'*sin(-pi/2+pi*(i-1)/L));
    beam(i)=20*log10(abs(w'*a));
end

% 作图
figure
angle=-90:180/200:(90-180/200);
plot(angle,beam);
grid on
xlabel('方向角/degree');
ylabel('幅度响应/dB');
figure
for k=1:N
    en(k)=(abs(e(k))).^2;
end
semilogy(en);hold on
xlabel('迭代次数');
ylabel('MSE');

运行结果：

figure2 LMS自适应波束形成算法仿真结果 figure3 LMS自适应波束形成算法误差图

(3). root-MUSIC算法MATLAB程序：

求根"MUSIC"，即"root-MUSIC"算法是"MUSIC"算法的一种多项式求根形式，其基本思想是"Pisarenko"分解。相比于"MUSIC"算法，"root-MUSIC"算法无须谱峰搜索，降低了复杂度。

clear all
close all
derad=pi/180;                           % 角度->弧度
radeg=180/pi;                           % 弧度->角度
twpi=2*pi;
kelm=8;                                 % 阵元数
dd=0.5;                                 % 阵元间距
d=0:dd:(kelm-1)*dd;
iwave=3;                                % 信源数
theta=[10 20 30];                       % 波达方向
snr=10;                                 % 信噪比
n=200;                                  % 采样数
A=exp(-1i*twpi*d.'*sin(theta*derad));   % 方向向量
S=randn(iwave,n);                       % 信源信号
X0=A*S;                                 % 接收信号
X=awgn(X0,snr,'measured');              % 添加噪声
Rxx=X*X';                               % 计算协方差矩阵
InvS=inv(Rxx);
[EVx,Dx]=eig(Rxx);                      % 特征值分解
EVAx=diag(Dx)';
[EVAx,Ix]=sort(EVAx);                   % 特征值从小到大排序
EVAx=fliplr(EVAx);                      % 左右翻转，从大到小排序
EVx=fliplr(EVx(:,Ix));
  
% root-MUSIC
Unx=EVx(:,iwave+1:kelm);                % 噪声子空间
syms z;
pz=z.^([0:kelm-1]');
pz1=(z^(-1)).^([0:kelm-1]);
fz=z.^(kelm-1)*pz1*Unx*Unx'*pz;         % 构造多项式
a=sym2poly(fz);                         % 符号多项式->数值多项式
zx=roots(a);                            % 求根
rx=zx';
[as,ad]=(sort(abs((abs(rx)-1))));
DOAest=asin(sort(-angle(rx(ad([1,3,5])))/pi))*180/pi

运行结果：

DOAest=
   9.9925    20.1515    29.9966

(4). 谱峰搜索传播算子算法MATLAB程序：

谱峰搜索传播算子算法利用方向向量和 $Q$ 矩阵的正交性构造谱函数，通过一维谱峰搜索，得到波达方向估计。相比于“MUSIC”算法，它无须对接数据的协方差矩阵进行特征值分解，降低了复杂度。算法的主要步骤如下：

1. 根据接收信号构造协方差矩阵 $R$。
2. 对协方差矩阵分块，计算传播算子 $P$。
3. 构造 $Q$ 矩阵。
4. 使 $\theta$ 变化，按照 $P_{PM}(\theta )=1/[a^H(\theta )Q{Q}^{H}a(\theta )]$ 计算谱函数，通过寻找谱峰来得到波达方向的估计值。

clear all
close all
derad=pi/180;
radeg=180/pi;
twpi=2*pi;
kelm=16;                                % 阵元数
dd=0.5;                                 % 阵元间距
d=0:dd:(kelm-1)*dd;
iwave=3;                                % 信源数
theta=[10 20 30];                       % DOA
pw=[1 0.8 0.7]';                        % 信号功率
nv=ones(1,kelm);                        % 归一化噪声方差
snr=20;                                 % 信噪比
snr0=10^(snr/10);               
n=200;                                  % 样本数量

A=exp(-1i*twpi*d.'*sin(theta*derad));   % 方向向量

K=length(d);
cr=zeros(K,K);
L=length(theta);
data=randn(L,n);
data=sign(data);
twpi=2.0*pi;
derad=pi/180.0;
s=diag(pw)*data;
A1=exp(-1i*twpi*d.'*sin([0:0.2:90]*derad));
received_signal=A*s;
cx=received_signal+diag(sqrt(nv/snr0/2))*(randn(K,n)+1i*randn(K,n));
Rxx=cx*cx'/n;

% PM算法
G=Rxx(:,1:iwave);
H=Rxx(:,iwave+1:end);
P=inv(G'*G)*G'*H;                       % 传播算子矩阵
Q=[P',-diag(ones(1,kelm-iwave))];       % Q矩阵

for iang=1:361
    angle(iang)=(iang-181)/2;
    phim=derad*angle(iang);
    a=exp(-1i*twpi*d*sin(phim)).';
    SP(iang)=1/(a'*Q'*Q*a);
end
SP=abs(SP);
SPmax=max(SP);
SP=10*log10(SP/SPmax);

% 作图
figure
h=plot(angle,SP,'-k');
set(h,'Linewidth',2)
xlabel('angle(degree)')
ylabel('magnitude(dB)')
axis([0 60 -60 0])
set(gca,'XTick',[0:10:60])
grid on
hold on
legend('Propagator Method')

运行结果：

figure4 谱峰搜索传播算子算法仿真运行结果

Reference

[1] 阵列信号处理及MATLAB实现（第2版）$\text{c◯}$ 张小飞 李建峰 徐大专 等 著.