前言
记录阵列信号处理的学习过程。
一、自适应波束形成的最佳权矢量
一般来说,并不希望直接求解方程,其理由如下:①由于移动用户环境是时变的,所以权向量的解必须能及时更新;②由于估计最佳解需要的数据是含噪声的,所以希望使用一种更新技术,可利用已求出的权向量求平滑最佳响应的估计,以减小噪声的影响。因此,希望使用自适应算法周期更新权向量。
实际上就是自适应滤波的工程,使得输出不断向期望信号收敛。这里以MMSE(最小均方误差)为例,说明如何把他变成一种自适应算法。
1.MMSE方法
MMSE准则就是使误差
y
(
k
)
−
d
(
k
)
y(k)-d(k)
y(k)−d(k)的均方值最小化,即代价函数取
J
(
w
q
)
=
E
[
∣
w
q
H
x
(
k
)
−
d
q
(
k
)
∣
2
]
J(w_q)=E[|w^H_qx(k) - d_q(k)|^2]
J(wq)=E[∣wqHx(k)−dq(k)∣2]实际上就是数学中的问题,给了y,要去求使得y取得最小最小值的x,在一维中我们直接求导数就可以,在二维和高维中,我们利用雅可比行列式,实际上是和一元差不多。在一元函数中,我们寻找最优点往往就是先给定初始
x
x
x,然后以一个步长往梯度的反方向步进,在这里也是一样。
我们将上式子展开
J
(
w
q
)
=
w
H
E
[
x
(
k
)
x
H
(
k
)
]
w
q
−
E
[
d
q
(
k
)
x
H
(
k
)
]
w
q
−
w
q
H
E
[
x
(
k
)
d
q
∗
(
k
)
]
+
E
[
d
q
(
k
)
d
q
∗
(
k
)
]
J(w_q)=w^HE[x(k)x^H(k)]w_q-E[d_q(k)x^H(k)]w_q-w_q^HE[x(k)d^*_q(k)]+E[d_q(k)d^*_q(k)]
J(wq)=wHE[x(k)xH(k)]wq−E[dq(k)xH(k)]wq−wqHE[x(k)dq∗(k)]+E[dq(k)dq∗(k)]对上式子对
w
q
w_q
wq求导数(矩阵求导,下次抽空了写一下),求得
∂
∂
w
q
J
(
w
q
)
=
2
E
[
x
(
k
)
x
H
(
k
)
]
w
q
−
2
E
[
x
(
k
)
d
q
∗
(
k
)
]
=
2
R
x
w
q
−
r
x
d
\frac{\partial}{\partial w_q}J(w_q)=2E[x(k)x^H(k)]w_q-2E[x(k)d^*_q(k)]=2R_xw_q-r_{xd}
∂wq∂J(wq)=2E[x(k)xH(k)]wq−2E[x(k)dq∗(k)]=2Rxwq−rxd
只要求出了梯度,我们就可以迭代,一步一步沿着梯度的方向去步进。
2.LMS算法
MMSE方法可以用LMS算法实现。
考虑随机梯度算法,其更新权矢量的一般公式为
w
q
(
k
+
1
)
=
w
q
(
k
)
−
1
2
u
∇
w_q(k+1) = w_q(k) - \frac{1}{2}u\nabla
wq(k+1)=wq(k)−21u∇
式中
∇
=
∂
∂
w
q
(
k
)
J
(
w
q
(
k
)
)
\nabla = \frac{\partial}{\partial w_q(k)}J(w_q(k))
∇=∂wq(k)∂J(wq(k)),
u
u
u是收敛因子,实际上就是步长,这个步长,选长了,容易收敛不了,就一直在最优点旁边左右移动,选短了,需要的迭代次数增加。
上面已经求出
∇
=
2
(
E
[
x
(
k
)
x
H
(
k
)
]
w
q
(
k
)
−
E
[
x
(
k
)
d
q
∗
(
k
)
]
)
\nabla = 2(E[x(k)x^H(k)]w_q(k) - E[x(k)d^*_q(k)])
∇=2(E[x(k)xH(k)]wq(k)−E[x(k)dq∗(k)])
但是
E
[
]
E[]
E[]这个东西是个什么呢?在实际运算中,我们如何求
E
E
E呢?实际上
E
E
E是一个总体的特征,当样本足够多时,我们可以去估计它,但是我们现在要求实时更新,也就是采一次样,那么就利用LMS算法去更新一次,现在有
E
E
E是不可以的,LMS算法的基本思路是把数学期望用各自的瞬时值代替,即得到
k
k
k时刻的梯度估计值如下
∇
~
(
k
)
=
2
x
(
k
)
[
x
H
(
k
)
w
q
(
k
)
−
d
q
∗
(
k
)
]
=
−
2
x
(
k
)
f
(
k
)
\tilde{\nabla}(k) = 2x(k)[x^H(k)w_q(k)-d^*_q(k)]=-2x(k)f(k)
∇~(k)=2x(k)[xH(k)wq(k)−dq∗(k)]=−2x(k)f(k)实际上
E
[
∇
~
]
=
∇
E[\tilde{\nabla} ] = \nabla
E[∇~]=∇,即
∇
~
\tilde{\nabla}
∇~是
∇
\nabla
∇的一个无偏估计,在某些情况下,我们是可以替代的。
我们将梯度值代回权矢量的更新公式,得到了LMS算法为
w
(
k
+
1
)
=
w
(
k
+
1
)
+
u
x
(
k
)
f
(
k
)
w(k+1) = w(k+1)+ux(k)f(k)
w(k+1)=w(k+1)+ux(k)f(k)
算法 | 初始化 | 更新公式 |
---|---|---|
LMS | w 0 = 0 w_0=0 w0=0 | y ( k ) = w H ( k ) x ( k ) f ( k ) = d ( k ) − y ( k ) w ( k + 1 ) = w ( k ) + u x ( k ) f ∗ ( k ) y(k)=w^H(k)x(k) \\ f(k)=d(k)-y(k) \\ w(k+1) = w(k) + ux(k)f^*(k) y(k)=wH(k)x(k)f(k)=d(k)−y(k)w(k+1)=w(k)+ux(k)f∗(k) |
二、MATLAB实现
1.代码逻辑
实际上LMS实现比较简单。
1.
k
k
k时刻计算
y
(
k
)
=
w
H
(
k
)
x
(
k
)
y(k)=w^H(k)x(k)
y(k)=wH(k)x(k)
2.计算
f
(
k
)
=
d
(
k
)
−
y
(
k
)
f(k)=d(k)-y(k)
f(k)=d(k)−y(k)
3.更新权矢量
w
(
k
+
1
)
=
w
(
k
)
+
u
x
(
k
)
f
∗
(
k
)
w(k+1) = w(k) + ux(k)f^*(k)
w(k+1)=w(k)+ux(k)f∗(k)
4.进行第
k
+
1
k+1
k+1次迭代
2.代码实现
clear;
clc;
M = 16;
thetas = [0 30 60];
lambda = 10;
d = lambda /2;
N = 1000;
n = 0:N-1;
f0 = 2000;
s = [1*sin(2*pi*f0 *n/(8*f0));...
2*sin(2*pi*2*f0 *n/(8*f0));...
3*sin(2*pi*3*f0 *n/(8*f0))
];
% 生成方向矢量
A = exp(-1i * 2 * pi * d * (0:M-1)' * sind(thetas) / lambda);
St = A*s + randn(M,N);
% LMS 算法 开始 进行自适应滤波
di = s(1,:); % 第一行为期望信号
u = 0.0005;
w = zeros(M,1); % 初始化权重向量
for k = 1:N
y(k) = w'*St(:,k);
e(k) = di(k) - y(k);
w = w + u * St(:,k) * conj(e(k));
end
scan_theta = [-89:90];
beam = zeros(1,length(scan_theta));
for i = 1 :length(scan_theta)
% 构造该方向的方向向量
v = exp(-1i * 2 * pi *d* (0:M-1)'.*sind(scan_theta(i))/lambda);
beam(i) = abs(w'*v);
end
figure;
plot(scan_theta,20*log10(beam/max(beam)))
title('方向图')
figure;
for k = 1:N
en(k) = (abs(e(k)))^2;
end
semilogy(en); hold on;
xlabel('迭代次数')
ylabel('MSE')
title('MSE')
3.实验结果
当以0度入射的信号为期望信号时
当以30度信号为期望信号时。
可见权值都收敛到了正确的方向图方向。
主要参考:张小飞.阵列信号处理及MATLAB实现[M].电子工业出版社