Projection Approximation Subspace Tracking PAST 算法

PAST 算法

建立模型

  在线估计一个映射矩阵 $\mathbf{W}\in {\mathbb{C}}^{M\times K}$，使得对于观测序列 $\mathbf{x}(i)\in {\mathbb{C}}^M$ 和对应的低维投影向量 $\mathbf{y}(i)\in {\mathbb{C}}^K$ 在 $t$ 时刻的最小二乘损失函数最小：
$$J_t(\mathbf{W})=\sum _{i=1}^t{\beta }^{t-i}{\Vert \mathbf{x}(i)-\mathbf{Wy}(i)\Vert }^2$$
其中 $\mathbf{y}(i)={\mathbf{W}}^H(t-1)\mathbf{x}(i)$ 以及 $0<\beta <1$。

闭式解与递推

  把 $J_t(\mathbf{W})$ 对 $\mathbf{W}$ 求导并令导数为零：
$$\frac{\partial J_t}{\partial {\mathbf{W}}^{\ast }}=-\sum _{i=1}^t{\beta }^{t-i}\left[\mathbf{x}(i)-\mathbf{Wy}(i)\right]{\mathbf{y}}^H(i)=0_{M\times K}$$
得到正规方程：
$$\sum _{i=1}^t{\beta }^{t-i}\mathbf{x}(i){\mathbf{y}}^H(i)=\mathbf{W}(t)\sum _{i=1}^t{\beta }^{t-i}\mathbf{y}(i){\mathbf{y}}^H(i)$$
即：
$$\mathbf{W}(t)={\mathbf{C}}_{xy}(t){\mathbf{C}}_{yy}^{-1}(t)$$
其中 ${\mathbf{C}}_{xy}(t)={\sum }_{i=1}^t{\beta }^{t-i}\mathbf{x}(i){\mathbf{y}}^H(i)$ 以及 ${\mathbf{C}}_{yy}(t)={\sum }_{i=1}^t{\beta }^{t-i}\mathbf{y}(i){\mathbf{y}}^H(i)$。
  由定义可得：
$$\begin{array}{rl}{\mathbf{C}}_{xy}(t)& =\beta {\mathbf{C}}_{xy}(t-1)+\mathbf{x}{\mathbf{y}}^H\\ {\mathbf{C}}_{yy}(t)& =\beta {\mathbf{C}}_{yy}(t-1)+\mathbf{y}{\mathbf{y}}^H\end{array}$$
其中 $\mathbf{x}=\mathbf{x}(t)$ 以及 $\mathbf{y}=y(t)$。
  接下来令 $\mathbf{P}(t)={\mathbf{C}}_{yy}^{-1}(t-1)$，对 $\mathbf{P}(t)$ 应用 Woodbury 定理：
$$\begin{array}{rl}\mathbf{P}(t)& ={\left(\beta {\mathbf{C}}_{yy}(t-1)+\mathbf{y}{\mathbf{y}}^H\right)}^{-1}\\ & ={\beta }^{-1}\left[\mathbf{P}(t-1)-\frac{\mathbf{P}(t-1)\mathbf{y}{\mathbf{y}}^H\mathbf{P}(t-1)}{\beta +{\mathbf{y}}^H\mathbf{P}(t-1)\mathbf{y}}\right]\end{array}$$
由此可得：
$$\mathbf{W}(t)={\mathbf{C}}_{xy}(t)\mathbf{P}(t)=\left[\beta {\mathbf{C}}_{xy}(t-1)+\mathbf{x}{\mathbf{y}}^H\right]\mathbf{P}(t)$$
其中：
$$\begin{array}{rl}\beta {\mathbf{C}}_{xy}(t-1)\mathbf{P}(t)& ={\mathbf{C}}_{xy}(t-1)\left[\mathbf{P}(t-1)-\frac{\mathbf{P}(t-1)\mathbf{y}{\mathbf{y}}^H\mathbf{P}(t-1)}{\beta +{\mathbf{y}}^H\mathbf{P}(t-1)\mathbf{y}}\right]\\ & =\mathbf{W}(t-1)\left[{\mathbf{I}}_K-\frac{\mathbf{y}{\mathbf{y}}^H\mathbf{P}(t-1)}{\beta +{\mathbf{y}}^H\mathbf{P}(t-1)\mathbf{y}}\right]\end{array}$$
以及：
$$\begin{array}{rl}\mathbf{P}(t)\mathbf{y}& ={\beta }^{-1}\left[\mathbf{P}(t-1)\mathbf{y}-\frac{\mathbf{P}(t-1)\mathbf{y}{\mathbf{y}}^H\mathbf{P}(t-1)\mathbf{y}}{\beta +{\mathbf{y}}^H\mathbf{P}(t-1)\mathbf{y}}\right]\\ & ={\beta }^{-1}\frac{\beta \mathbf{P}(t-1)\mathbf{y}+\mathbf{P}(t-1)\mathbf{y}{\mathbf{y}}^H\mathbf{P}(t-1)\mathbf{y}-\mathbf{P}(t-1)\mathbf{y}{\mathbf{y}}^H\mathbf{P}(t-1)\mathbf{y}}{\beta +{\mathbf{y}}^H\mathbf{P}(t-1)\mathbf{y}}\\ & =\frac{\mathbf{P}(t-1)\mathbf{y}}{\beta +{\mathbf{y}}^H\mathbf{P}(t-1)\mathbf{y}}\end{array}$$
所以有：
$$\begin{array}{rl}\mathbf{W}(t)& =\mathbf{W}(t-1)-\frac{\mathbf{W}(t-1)\mathbf{y}{\mathbf{y}}^H\mathbf{P}(t-1)}{\beta +{\mathbf{y}}^H\mathbf{P}(t-1)\mathbf{y}}+\frac{\mathbf{x}{\mathbf{y}}^H\mathbf{P}(t-1)}{\beta +{\mathbf{y}}^H\mathbf{P}(t-1)\mathbf{y}}\\ & =\mathbf{W}(t-1)+\frac{\left[\mathbf{x}-\mathbf{W}(t-1)\mathbf{y}\right]{\mathbf{y}}^H\mathbf{P}(t-1)}{\beta +{\mathbf{y}}^H\mathbf{P}(t-1)\mathbf{y}}\end{array}$$
  最后整理得：
$$\begin{array}{rl}\mathbf{P}(t)& ={\beta }^{-1}\left[\mathbf{P}(t-1)-\frac{\mathbf{P}(t-1)\mathbf{y}{\mathbf{y}}^H\mathbf{P}(t-1)}{\beta +{\mathbf{y}}^H\mathbf{P}(t-1)\mathbf{y}}\right]\\ \mathbf{W}(t)& =\mathbf{W}(t-1)+\frac{\left[\mathbf{x}-\mathbf{W}(t-1)\mathbf{y}\right]{\mathbf{y}}^H\mathbf{P}(t-1)}{\beta +{\mathbf{y}}^H\mathbf{P}(t-1)\mathbf{y}}\end{array}$$

算法步骤

$$\begin{array}{rl}\mathbf{y}(t)& ={\mathbf{W}}^H(t-1)\mathbf{x}(t)\\ \mathbf{h}(t)& =\mathbf{P}(t-1)\mathbf{y}(t)\\ \mathbf{g}(t)& =\mathbf{h}(t)/\left[\beta +{\mathbf{y}}^H(t)\mathbf{h}(t)\right]\\ \mathbf{P}(t)& =\frac{1}{\beta }\left[\mathbf{P}(t-1)-\mathbf{g}(t){\mathbf{h}}^H(t)\right]\\ \mathbf{e}(t)& =\mathbf{x}(t)-\mathbf{W}(t-1)\mathbf{y}(t)\\ \mathbf{W}(t)& =\mathbf{W}(t-1)+\mathbf{e}(t){\mathbf{g}}^H(t)\end{array}$$

O-PAST 算法

建立模型

  OPAST 算法在每次迭代之后添加一步正交化操作：
$$\mathbf{W}(t):=\mathbf{W}(t){\left[{\mathbf{W}}^H(t)\mathbf{W}(t)\right]}^{-\frac{1}{2}}$$
正交化保证 ${\mathbf{W}}^H(t)\mathbf{W}(t)={\mathbf{I}}_K$。

闭式解与递推

  易得：
$$\begin{array}{rl}{\mathbf{W}}^H(t)\mathbf{W}(t)& =\left[{\mathbf{W}}^H(t-1)+\mathbf{g}(t){\mathbf{e}}^H(t)\right]\left[\mathbf{W}(t-1)+\mathbf{e}(t){\mathbf{g}}^H(t)\right]\\ & ={\mathbf{I}}_K+\Vert \mathbf{e}(t){\Vert }^2\mathbf{g}(t){\mathbf{g}}^H(t)={\mathbf{I}}_K+\mathbf{z}(t){\mathbf{z}}^H(t)\end{array}$$
其中 $\mathbf{z}(t)=\Vert \mathbf{e}(t)\Vert \mathbf{g}(t)$。上式的化简利用了 ${\mathbf{W}}^H(t-1)\mathbf{W}(t-1)={\mathbf{I}}_K$ 以及：
$$\begin{array}{rl}{\mathbf{W}}^H(t-1)\mathbf{e}(t)& ={\mathbf{W}}^H(t-1)\mathbf{x}(t)-{\mathbf{W}}^H(t-1)\mathbf{W}(t-1)\mathbf{y}(t)\\ & ={\mathbf{W}}^H(t-1)\mathbf{x}(t)-\mathbf{y}(t)=0_{K\times 1}\end{array}$$
  进一步可得到：
$$\begin{array}{rl}{\left[{\mathbf{W}}^H(t)\mathbf{W}(t)\right]}^{-\frac{1}{2}}& ={\mathbf{I}}_K+\frac{1}{\Vert \mathbf{z}(t){\Vert }^2}\left(\frac{1}{\sqrt{1+\Vert \mathbf{z}(t){\Vert }^2}}\right)\mathbf{z}{\mathbf{z}}^H\\ & ={\mathbf{I}}_K+\tau (t)\mathbf{g}(t){\mathbf{g}}^H(t)\end{array}$$
其中：
$$\tau (t)=\frac{1}{\Vert \mathbf{g}(t){\Vert }^2}\left(\frac{1}{\sqrt{1+\Vert \mathbf{e}(t){\Vert }^2\Vert \mathbf{g}(t){\Vert }^2}}-1\right)$$
由此可得：
$$\begin{array}{rl}& \mathbf{W}(t):=\mathbf{W}(t){\left[{\mathbf{W}}^H(t)\mathbf{W}(t)\right]}^{-\frac{1}{2}}\\ =& \left[\mathbf{W}(t-1)+\mathbf{e}(t){\mathbf{g}}^H(t)\right]\left[{\mathbf{I}}_K+\tau (t)\mathbf{g}(t){\mathbf{g}}^H(t)\right]\\ =& \mathbf{W}(t-1)+\tau (t)\mathbf{W}(t-1)\mathbf{g}(t){\mathbf{g}}^H(t)+\left[1+\tau (t)\Vert \mathbf{g}(t){\Vert }^2\right]\mathbf{e}(t){\mathbf{g}}^H(t)\\ =& \mathbf{W}(t-1)+\tau (t)\mathbf{W}(t-1)\mathbf{g}(t){\mathbf{g}}^H(t)+\frac{\mathbf{e}(t){\mathbf{g}}^H(t)}{\sqrt{1+\Vert \mathbf{e}(t){\Vert }^2\Vert \mathbf{g}(t){\Vert }^2}}\end{array}$$

代码实现

  PAST 与 OPAST 算法可作为 DoA 估计的预处理步骤，其估计得到的矩阵 $\mathbf{W}\in {\mathbb{C}}^{M\times K}$ 可直接用于后续的 ESPRIT 算法进行 DoA 估计。

clc; clear; close all;M = 16;
theta = [0,5,10];
K = length(theta);
L = 10240;
SNRdB = -10;S = (randn(K,L)+1j*randn(K,L))/sqrt(2);
N = 1*sqrt(10^(-SNRdB/10)/2) * (randn(M,L) + 1j*randn(M,L));A = exp(-1j*pi*(0:M-1)'*sind(theta));
X = A*S+N;beta = 1;
DOA_past = PAST(X,K,beta); disp(DOA_past);
DOA_opast = OPAST(X,K,beta); disp(DOA_opast);function [DOA] = PAST(X, K, beta)[M,L] = size(X);% W = randn(M,K) + 1j*randn(M,K);  % random initialW = eye(M);W = W(:,1:K);P = eye(K);for t = 1:Lx = X(:,t);y = W' * x;                   % K x 1h = P * y;                    % K x 1g = h / (beta + y' * h);      % K x 1P = (P - g * h') / beta;      % K x Ke = x - W * y;                % M x 1W = W + e * g';               % M x KendUs = W;U1 = Us(1:M-1,:);U2 = Us(2:M,:);Psi = pinv(U1)*U2;[~,Phi] = eig(Psi,'vector');DOA = sort(asind(-angle(Phi)/pi));
endfunction [DOA] = OPAST(X, K, beta)[M,L] = size(X);% W = randn(M,K) + 1j*randn(M,K);  % random initialW = eye(M);W = W(:,1:K);P = eye(K);for t = 1:Lx = X(:,t);y = W' * x;                   % K x 1h = P * y;                    % K x 1g = h / (beta + y' * h);      % K x 1P = (P - g * h') / beta;      % K x Ke = x - W * y;                % M x 1W = W + e * g';               % M x Ktau = 1/(g'*g)*(1/sqrt(1+(e'*e)*(g'*g))-1);W = W*(eye(K)+tau*g*g');endUs = W;U1 = Us(1:M-1,:);U2 = Us(2:M,:);Psi = pinv(U1)*U2;[~,Phi] = eig(Psi,'vector');DOA = sort(asind(-angle(Phi)/pi));
end