Matrix pencil矩阵铅笔算法(原始论文记录与复现)（二）

《Estimating two-dimensional frequencies by matrix enhancement and matrix pencil》¹
这篇上一部分见

---

上一部分

Matrix pencil矩阵铅笔算法(原始论文记录与复现)(一)

本文的补充

MEMP的pairing部分

从 { y i ; i = 1 , ⋯   , I } , { z i ; i = 1 , ⋯   , I } \left\{ y_i;i=1,\cdots ,I \right\} ,\left\{ z_i;i=1,\cdots ,I \right\} {yi​;i=1,⋯,I},{zi​;i=1,⋯,I}中选择对应的 { ( y i , z i ) ; i = 1 , ⋯   , I } . \left\{ \left( y_i,z_i \right) ;i=1,\cdots ,I \right\}. {(yi​,zi​);i=1,⋯,I}. 去
$$Maximize⟶J_n\left(i,j\right)=\sum _{i=1}^I||{\mathbf{u}}_i^H{\mathbf{e}}_L\left(y_i,z_j\right)|{|}^2\left(1\right)$$
具体做法如下：
$$\begin{array}{r}1:\text{遍历}y\text{令}i=1\\ 2:\text{计算}{J}_n\left(i,j\right)forj=1,2,\cdots ,I\\ 3:\text{找到}2\text{中使得}{J}_n\left(i,j\right)\text{最大的}j,\text{得到}index\left(i,j\left(i\right)\right)\text{于是得到}pair\left(y_i,z_{j\left(i\right)}\right)\\ 4:i⟵i+1\\ 5:\text{计算}{J}_n\left(i,j\right)forj=1,2,\cdots ,Ibutj\ne j\left(k\right)wherek=1,2,\cdots ,i-1\\ 6:\text{找到}5\text{中使得}{J}_n\left(i,j\right)\text{最大的}j,\text{得到}index\left(i,j\left(i\right)\right)\text{于是得到}pair\left(y_i,z_{j\left(i\right)}\right)\\ 7:\text{回到}4\text{直至}i=I-1\end{array}$$

MEMP算法完整步骤

s t e p 1 : 从带噪声的 x ′ ( m , n ) 构造 X e ′ ( 31 ) , K 和 L 满足 ( 42 ) , 最好满足 ( 41 ) s t e p 2 : 对 X e ′ 做 S V D 得到 U s ′ ，得到 U 1 ′ 、 U 2 ′ 、 U 1 P ′ 、 U 2 P ′ s t e p 3 : 计算 U 2 ′ − λ U 1 ′ 和 U 2 P ′ − λ U 1 P ′ 的 G E    ( Q Z    a l g o r i t h m    t o    c o m p u t e    t h e    G E s    o f    U 1 ′ H U 2 ′ − λ U 1 ′ H U 1 ′ a n d U 1 P ′ H U 2 P ′ − λ U 1 P ′ H U 1 P ′ ) s t e p 4 : p a i r i n g    { ( y i , z i ) ; i = 1 , ⋯   , I }    t o    m a x m i z e ( 44 ) s t e p 5 : C a l c u l a t e    { ( f 1 i , f 2 i ) ; i = 1 , ⋯   , I }    f r o m    e s t i m a t e d    { ( y i , z i ) ; i = 1 , ⋯   , I }    w i t h    exp ⁡ e x p r e s s i o n    f o r m . step1:\text{从带噪声的}x^{\prime}\left( m,n \right) \text{构造}\boldsymbol{X}_{e}^{\prime}\left( 31 \right) ,K\text{和}L\text{满足}\left( 42 \right) ,\text{最好满足}\left( 41 \right) \\ step2:\text{对}\boldsymbol{X}_{e}^{\prime}\text{做}SVD\text{得到}\boldsymbol{U}_{s}^{\prime}\text{，得到}\boldsymbol{U}_{1}^{\prime}\text{、}\boldsymbol{U}_{2}^{\prime}\text{、}\boldsymbol{U}_{1P}^{\prime}\text{、}\boldsymbol{U}_{2P}^{\prime} \\ step3:\text{计算}\boldsymbol{U}_{2}^{\prime}-\lambda \boldsymbol{U}_{1}^{\prime}\text{和}\boldsymbol{U}_{2P}^{\prime}-\lambda \boldsymbol{U}_{1P}^{\prime}\text{的}GE \\ \,\, \left( QZ\,\,algorithm\,\,to\,\,compute\,\,the\,\,GEs\,\,of\,\,{\boldsymbol{U}_{1}^{\prime}}^H\boldsymbol{U}_{2}^{\prime}-\lambda {\boldsymbol{U}_{1}^{\prime}}^H\boldsymbol{U}_{1}^{\prime}and{\boldsymbol{U}_{1P}^{\prime}}^H\boldsymbol{U}_{2P}^{\prime}-\lambda {\boldsymbol{U}_{1P}^{\prime}}^H\boldsymbol{U}_{1P}^{\prime} \right) \\ step4:pairing\,\,\left\{ \left( y_i,z_i \right) ;i=1,\cdots ,I \right\} \,\,to\,\,maxmize\left( 44 \right) \\ step5:Calculate\,\,\left\{ \left( {f_1}_i,{f_2}_i \right) ;i=1,\cdots ,I \right\} \,\,from\,\,estimated\,\,\left\{ \left( y_i,z_i \right) ;i=1,\cdots ,I \right\} \,\,with\,\,\exp expression\,\,form. step1:从带噪声的x′(m,n)构造Xe′​(31),K和L满足(42),最好满足(41)step2:对Xe′​做SVD得到Us′​，得到U1′​、U2′​、U1P′​、U2P′​step3:计算U2′​−λU1′​和U2P′​−λU1P′​的GE(QZalgorithmtocomputetheGEsofU1′​HU2′​−λU1′​HU1′​andU1P′​HU2P′​−λU1P′​HU1P′​)step4:pairing{(yi​,zi​);i=1,⋯,I}tomaxmize(44)step5:Calculate{(f1​i​,f2​i​);i=1,⋯,I}fromestimated{(yi​,zi​);i=1,⋯,I}withexpexpressionform.
此外，${\mathbf{X}}_e$有一种更加增强后的矩阵形式——
$$\begin{gathered} {\mathbf{X}}_{ee}=\left[{\mathbf{X}}_e,{\mathbf{P}}_e{\mathbf{X}}_e^{\ast }\right] \\ {\mathbf{P}}_e=\left[\begin{array}{cccc}& & & 1\\ & & 1& \\ & \cdots & & \\ 1& & & \end{array}\right] \\ 因为range\left({\mathbf{X}}_{ee}\right)=range\left({\mathbf{X}}_e\right)=range\left({\mathbf{E}}_L\right) \end{gathered}$$

实验仿真

我们考虑这样的信号模型
$$\begin{gathered} x^{\prime }\left(m;n\right)=\sum _{i=1}^3\exp \left(j2\pi f_{1i}m+j2\pi f_{2i}n\right)+w\left(m;n\right) \\ 0⩽m⩽19,0⩽n⩽19. \\ \left(f_{11},f_{21}\right)=\left(0.26,0.24\right) \\ \left(f_{12},f_{22}\right)=\left(0.24,0.24\right) \\ \left(f_{13},f_{23}\right)=\left(0.24,0.26\right) \end{gathered}$$

每次独立仿真200次测试，
代码

I = 3;
M= 20;
N = 20;
f1=[0.26,0.24,0.24];
f2=[0.24,0.24,0.26];
% ploting
f = figure;
xlabel('f1');
ylabel('f2');
title('SNR=20db,K=6,L=6');
xlim([0.2,0.3]);
ylim([0.2,0.3]);
line([f1(1),f1(1)],[0,1],'linestyle','--');
hold on;
line([f1(2),f1(2)],[0,1],'linestyle','--');
hold on;
line([f1(3),f1(3)],[0,1],'linestyle','--');
hold on;
line([0,0.4],[f2(1),f2(1)],'linestyle','--');
hold on;
line([0,0.4],[f2(2),f2(2)],'linestyle','--');
hold on;
line([0,0.4],[f2(3),f2(3)],'linestyle','--');
hold on;
scatter(f1(1),f2(1),'filled','r');
scatter(f1(2),f2(2),'filled','r');
scatter(f1(3),f2(3),'filled','r');
x = zeros(M,N);
SNR = 20; %db
noise_var = 1/(10^((SNR/10)));
for idx_m = 1:M
    for idx_n = 1:N
        for idx_i = 1:I
            x(idx_m, idx_n) = x(idx_m, idx_n)+ exp(1j*2*pi*f1(idx_i)*(idx_m-1)+1j*2*pi*f2(idx_i)*(idx_n-1));
            %x(idx_m, idx_n) = x(idx_m, idx_n)+ exp(1j*2*pi*f1(idx_i)*(idx_m-1)+1j*2*pi*f2(idx_i)*(idx_n-1)) 
        end
    end
end
sim_times = 200;
for idx_simluate = 1:sim_times
    noise = sqrt(noise_var/2) * (randn(M,N)+1j*randn(M,N));
    x_withnoise = x + noise;
   
    [u,s,v] = svd(x_withnoise);
    K = 6;
    L = 6;
    % get X_e_bar(X_en) with size of(KL,(M-K+1)(N-L+1))
    X_en =zeros(K*L,(M-K+1)*(N-L+1));

    %preprocess X_0 to X_M-1 (M blcok)
    X_m = zeros(M,L,N-L+1);
    for idx_m = 1:M
        for idx_L = 1:L
            X_m(idx_m,idx_L,:) = x_withnoise(idx_m,idx_L:N-L+idx_L);
        end
    end
    for idx_row = 1:K
        for idx_col = 1:M-K+1
            X_en((idx_row-1) * L + 1 : idx_row * L,(idx_col - 1) * (N-L+1)+1: idx_col * (N-L+1)) = X_m(idx_row + idx_col -1, :,:);
        end
    end
    P = zeros(K*L,K*L);
    for idx_p = 1:K*L
        P(idx_p,K*L-idx_p+1) = 1;
    end
    X_ee = cat(2,X_en,P*conj(X_en));
    [U,S,V]= svd(X_ee);
    Us = U(:,1:I); %(KL,I)
    Un = U(:,I+1:K*L);%(KL,min-I)
    U1 = Us(1:(K-1)*L,:);
    U2 = Us(L+1:K*L,:);
    lambda_y = eig(U1'*U2,U1'*U1); % yi extraction
    y= mod(angle(lambda_y)/2/pi+1,1); % % f1i extraction

    P =zeros(K*L,K*L);
    for i = 1:K*L
        P(i, ceil(i/K)+mod(i-1,K)*L ) = 1;
    end
    Usp = P*Us;%(KL,I)
    U1p = Usp(1:(L-1)*K,:);
    U2p = Usp(K+1:K*L,:);
    lambda_z = eig(U1p'*U2p,U1p'*U1p);
    z= mod(angle(lambda_z)/2/pi+1,1);

    %pairing
    YL = zeros(K,I);
    ZL = zeros(K,I);
    for i = 1:K
        YL(i,:) = lambda_y.^(i-1);
    end
    for i = 1:K
        ZL(i,:) = lambda_z.^(i-1);
    end
    temp= kron(YL(:,1),ZL(:,1)); % (L*L,1)
    % pairing
    z_flag = zeros(1,I); % indicate select_z is selected or not.
    point_yz =zeros(1,I); % point_yz[y] = z
    for idx_y = 1:I
        yi = YL(:,idx_y);
        base_value = -Inf;
        J = zeros(1,I);
        for select_z = 1:I
            if z_flag(select_z)==1
                continue;
            else
                zi = ZL(:,select_z);
                eL = kron(yi,zi); %(KL,1);
                for un_idx = 1:I
                    J(select_z) = J(select_z) + (abs(Us(:,un_idx)'*eL))^2;
                end
            end

        end
        [m,i] = max(J);
        point_yz(idx_y) = i;
        z_flag(point_yz(idx_y)) = 1;
    end
    scatter(y(1),z(point_yz(1)),'k');
    scatter(y(2),z(point_yz(2)),'k');
    scatter(y(3),z(point_yz(3)),'k');  
end
grid on
ax = gca; % 获取当前坐标轴对象
ax.XGrid = 'on'; % 打开X轴网格
ax.YGrid = 'on'; % 打开Y轴网格
%ax.GridLineStyle = '-'; % 设置网格线样式
ax.XMinorGrid = 'on'; % 打开X轴次要网格
ax.YMinorGrid = 'on'; % 打开Y轴次要网格
ax.MinorGridLineStyle = ':'; % 设置次要网格线样式
%ax.GridAlpha = 0.5; % 设置网格线透明度

SNR=20db:

10db:

---

其他问题

估计$\mathbf{A}$

首先回顾一下矩阵广义逆、违逆的相关知识。
直接甩结论！！！
违逆是广义逆的特殊形式
·左逆和右逆都违逆的特殊形式
*当矩阵列满秩时，左逆存在且为违逆。
*当矩阵行满秩时，右逆存在且为违逆。
广义逆（Generalized Inverse）：${\mathbf{A}}^{-}$是$\mathbf{A}$的广义逆如果$${\mathbf{AA}}^{-}\mathbf{A}=\mathbf{A}$$
广义逆不唯一。
Moore-Penrose逆：$\mathbf{X}$是$\mathbf{A}$的Moore-Penrose逆，如果满足一下Moore-Penrose条件：
1: $\mathbf{AXA}=\mathbf{A}$
2: $\mathbf{XAX}=\mathbf{X}$
3: ${\left(\mathbf{AX}\right)}^H=\mathbf{AX}$
4: ${\left(\mathbf{XA}\right)}^H=\mathbf{XA}$
其中${\left(\cdot \right)}^H$表示Hermitian共轭转置，${\left(\cdot \right)}^H={\stackrel{-}{\left(\cdot \right)}}^T$,可以证明Moore-Penrose逆存在且唯一。第一条就是广义逆的定义，因为Moore-Penrose逆是特殊的广义逆。
求解方法：
$$\begin{gathered} \text{设}r=rank\left(A_{m\times n}\right) \\ \text{将}A\text{分解为}{A}_{m\times n}=P_{m\times m}{\left[\begin{array}{cc}{I}_r& O\\ O& O\end{array}\right]}_{m\times n}{Q}_{n\times n},\text{其中}P,Q\text{可逆},I_r\text{表示}r\text{阶单位阵},A\text{可逆时},O\text{为}0\text{阶阵} \\ \text{令}{R}_{m\times r}=P_{m\times m}{\left[\begin{array}{c}{I}_r\\ O\end{array}\right]}_{m\times r},S={\left[I_{r}O\right]}_{r\times n}{Q}_{n\times n} \\ \text{则}{X}_{n\times m}={S^H}_{n\times r}{{\left({R^H}_{r\times m}{A}_{m\times n}{S^H}_{n\times r}\right)}^{-1}}_{r\times r}{R^H}_{r\times m} \end{gathered}$$
违逆（Pseudo-inverse）：
$$\begin{gathered} X\text{是}A\text{的违逆，如果} \\ X=V\Xi U^H \\ \text{其中酉矩阵}U,V\text{来自对}A\text{的}SVD\text{分解}A=U\Sigma V^H \\ \text{设}\Sigma =\left[\begin{array}{cccc}{\sigma }_1& & & \\ & \ddots & & \\ & & {\sigma }_r& \\ & & & O\end{array}\right],\text{则}\Xi =\left[\begin{array}{cccc}\frac{1}{{\sigma }_1}& & & \\ & \ddots & & \\ & & \frac{1}{{\sigma }_r}& \\ & & & O\end{array}\right]\left(A\text{可逆时，}O\text{为}0\text{阶阵}\right) \end{gathered}$$
上述违逆内容来自于Gil的18.06课程，Gil介绍的这种求违逆的方法就是求Moore-Penrose逆的方法之一。
将上式子带入4个条件都是成立的。
$$\begin{gathered} AXA=A:\text{左}=U\Sigma V^{H}V\Xi U^{H}U\Sigma V^H=U\Sigma V^H=A \\ XAX=X:\text{左}=V\Xi U^{H}U\Sigma V^{H}V\Xi U^H=V\Xi U^H=X \\ {\left(AX\right)}^H=AX:\text{左}={\left(U\Sigma V^{H}V\Xi U^H\right)}^H=U\Sigma V^{H}V\Xi U^H \\ {\left(XA\right)}^H=XA:\text{左}={\left(V\Xi U^{H}U\Sigma V^H\right)}^H=V\Xi U^{H}U\Sigma V^H \end{gathered}$$
同时因为Moore-Penrose逆对任意矩阵都存在且唯一，而违逆对任意矩阵都存在(因为SVD对任意矩阵都存在)因此违逆是且只能是Moore-Penrose逆。
$$\begin{gathered} \text{设}r=rank\left(A_{m\times n}\right) \\ \text{将}A\text{分解为}{A}_{m\times n}=P_{m\times m}{\left[\begin{array}{cc}{I}_r& O\\ O& O\end{array}\right]}_{m\times n}{Q}_{n\times n},\text{其中}P,Q\text{可逆},I_r\text{表示}r\text{阶单位阵},A\text{可逆时},O\text{为}0\text{阶阵} \\ \text{令}{R}_{m\times r}=P_{m\times m}{\left[\begin{array}{c}{I}_r\\ O\end{array}\right]}_{m\times r},S={\left[I_{r}O\right]}_{r\times n}{Q}_{n\times n} \\ \text{则}{X}_{n\times m}={S^H}_{n\times r}{{\left({R^H}_{r\times m}{A}_{m\times n}{S^H}_{n\times r}\right)}^{-1}}_{r\times r}{R^H}_{r\times m} \\ \text{将}A\text{做}SVD\text{，有}A=U\Sigma V^H \\ \Sigma =\left[\begin{array}{cccc}{\sigma }_1& & & \\ & \ddots & & \\ & & {\sigma }_r& \\ & & & O\end{array}\right],\text{则}\Xi =\left[\begin{array}{cccc}\frac{1}{{\sigma }_1}& & & \\ & \ddots & & \\ & & \frac{1}{{\sigma }_r}& \\ & & & O\end{array}\right]\left(A\text{可逆时，}O\text{为}0\text{阶阵}\right) \\ \text{现在证明一下}{X}_{n\times m}={S^H}_{n\times r}{{\left({R^H}_{r\times m}{A}_{m\times n}{S^H}_{n\times r}\right)}^{-1}}_{r\times r}{R^H}_{r\times m}=V\Xi U^H\text{这个等式成立}\left(Moore-Penrose\text{逆}=\text{违逆}\right) \\ prove:\text{我们的目标是证明}{X}_{n\times m}=S^H{\left(R^{H}A{S}^H\right)}^{-1}{R}^H=V\Xi U^H\text{即可} \\ \text{左边}=S^H{\left(R^{H}U\Sigma V^H{S}^H\right)}^{-1}{R}^H=S^H{\left(S^H\right)}^{-1}{\left(V^H\right)}^{-1}\Xi U^{-1}{\left(R^H\right)}^{-1}{R}^H={\left(V^H\right)}^{-1}\Xi U^{-1}=V\Xi U^H\left(\text{最后一步用到酉矩阵的性质}{V}^H=V^{-1},U^H=U^H\right) \end{gathered}$$
违逆和Moore-Penrose逆用$A^{†}$表示
左逆与右逆：
$$\begin{gathered} if{A}_{m\times n}\text{列满秩}\left(r=n\right),\text{左逆}{A}_{left}^{-1}\text{存在且}{A_{left}^{-1}}_{n\times m}={{\left({A^H}_{n\times m}{A}_{m\times n}\right)}^{-1}}_{n\times n}{A^H}_{n\times m}\text{此时有}{A}_{left}^{-1}A={\left(A^{H}A\right)}^{-1}{A}^{H}A=I \\ if{A}_{m\times n}\text{行满秩}\left(r=m\right),\text{右逆}{A}_{right}^{-1}\text{存在且}={A_{right}^{-1}}_{n\times m}={A^H}_{n\times m}{{\left(A_{m\times n}{A^H}_{n\times m}\right)}^{-1}}_{m\times m}\text{此时有}A{A}_{right}^{-1}=A{A}^H{\left(A{A}^H\right)}^{-1}=I \end{gathered}$$

---

上一篇文章中有一旦y和z被估计出来，下面两个矩阵可以得到
$$\begin{gathered} {\mathbf{Z}}_L={\left[\begin{array}{cccc}{z}_1^0& z_2^0& \cdots & z_I^0\\ z_1& z_2& \cdots & z_I\\ ⋮& ⋮& \cdots & ⋮\\ z_1^{L-1}& z_2^{L-1}& \cdots & z_I^{L-1}\end{array}\right]}_{L\times I} \\ {\mathbf{Z}}_R={\left[\begin{array}{cccc}{z}_1^0& z_1& \cdots & z_1^{N-L}\\ z_2^0& z_2& \cdots & z_2^{N-L}\\ ⋮& ⋮& \cdots & ⋮\\ z_I^0& z_I& \cdots & z_I^{N-L}\end{array}\right]}_{I\times N-L+1} \end{gathered}$$

$$\begin{gathered} {\mathbf{X}}_e={\mathbf{E}}_L{\mathbf{AE}}_R \\ {\mathbf{E}}_L={\left[\begin{array}{l}{\mathbf{Z}}_L\\ {\mathbf{Z}}_L{\mathbf{Y}}_d\\ \cdots \\ {\mathbf{Z}}_L{\mathbf{Y}}_d^{K-1}\end{array}\right]}_{KL\times I} \\ {\mathbf{E}}_R={\left[{\mathbf{Z}}_R,{\mathbf{Y}}_d{\mathbf{Z}}_R,\cdots ,{\mathbf{Y}}_d^{M-K}{\mathbf{Z}}_R\right]}_{I\times \left(N-L+1\right)\left(M-K+1\right)} \end{gathered}$$
从${\mathbf{X}}_e$式子可以得到
$$\mathbf{A}={\mathbf{E}}_{\mathbf{L}}^{+}{\mathbf{X}}_{\mathbf{e}}{\mathbf{E}}_{\mathbf{R}}^{+}$$
其中${\mathbf{E}}_{\mathbf{L}}^{+}$为左逆，由于${\mathbf{E}}_{\mathbf{L}}$列满秩，左逆存在且为${\mathbf{E}}_{\mathbf{L}}^{+}={\left({\mathbf{E}}_L^H{\mathbf{E}}_L\right)}^{-1}{\mathbf{E}}_L^H$
其中${\mathbf{E}}_{\mathbf{R}}^{+}$为右逆，由于${\mathbf{E}}_{\mathbf{R}}$行满秩，右逆存在且为${\mathbf{E}}_R^{+}={\mathbf{E}}_R^H{\left({\mathbf{E}}_R{\mathbf{E}}_R^H\right)}^{-1}$

估计A的仿真代码

I = 3;
M= 20;
N = 20;
% f1=[0.26,0.24,0.24];
% f2=[0.24,0.24,0.26];
f1=[0.1,0.2,0.3];
y_before = exp(1j*2*pi*f1);
f2=[0.4,0.5,0.6];
z_before = exp(1j*2*pi*f2);
a = [1+1j,3+4j,6+2j];
%a = [1,1,1];
% ploting
f = figure;
xlabel('f1');
ylabel('f2');
title('SNR=10db,K=3,L=3');
xlim([0.2,0.3]);
ylim([0.2,0.3]);
line([f1(1),f1(1)],[0,1],'linestyle','--');
hold on;
line([f1(2),f1(2)],[0,1],'linestyle','--');
hold on;
line([f1(3),f1(3)],[0,1],'linestyle','--');
hold on;
line([0,0.4],[f2(1),f2(1)],'linestyle','--');
hold on;
line([0,0.4],[f2(2),f2(2)],'linestyle','--');
hold on;
line([0,0.4],[f2(3),f2(3)],'linestyle','--');
hold on;
scatter(f1(1),f2(1),'filled','r');
scatter(f1(2),f2(2),'filled','r');
scatter(f1(3),f2(3),'filled','r');


x = zeros(M,N);
SNR = 30; %db
noise_var = 1/(10^((SNR/10)));
for idx_m = 1:M
    for idx_n = 1:N
        for idx_i = 1:I
            x(idx_m, idx_n) = x(idx_m, idx_n)+ a(idx_i)*exp(1j*2*pi*f1(idx_i)*(idx_m-1)+1j*2*pi*f2(idx_i)*(idx_n-1));
            %x(idx_m, idx_n) = x(idx_m, idx_n)+ exp(1j*2*pi*f1(idx_i)*(idx_m-1)+1j*2*pi*f2(idx_i)*(idx_n-1)) 
        end
    end
end
sim_times = 1;
for idx_simluate = 1:sim_times
    noise = sqrt(noise_var/2) * (randn(M,N)+1j*randn(M,N));
    x_withnoise = x;
    [u,s,v] = svd(x_withnoise);
    K = 3;
    L = 3;
    % get X_e_bar(X_en) with size of(KL,(M-K+1)(N-L+1))
    X_en =zeros(K*L,(M-K+1)*(N-L+1));

    %preprocess X_0 to X_M-1 (M blcok)
    X_m = zeros(M,L,N-L+1);
    for idx_m = 1:M
        for idx_L = 1:L
            X_m(idx_m,idx_L,:) = x_withnoise(idx_m,idx_L:N-L+idx_L);
        end
    end
    for idx_row = 1:K
        for idx_col = 1:M-K+1
            X_en((idx_row-1) * L + 1 : idx_row * L,(idx_col - 1) * (N-L+1)+1: idx_col * (N-L+1)) = X_m(idx_row + idx_col -1, :,:);
        end
    end
    % check if satisfies X_en == E_L *A*E_R
    % min(X_en-E_L_before*diag(a)*E_R_before)
    Z_L_before = zeros(L,I);
    Z_R_before = zeros(I,N-L+1);
    Y_d_before = diag(y_before);
    for idx_L = 1:L
        Z_L_before(idx_L,:) = z_before.^(idx_L -1);
    end
    for idx_R = 1:N-L+1
        Z_R_before(:,idx_R) = (z_before.').^(idx_R-1);
    end
    E_L_before = zeros(K*L,I);
    E_R_before = zeros(I,(N-L+1)*(M-K+1));
    for idx_k = 1:K
        E_L_before(1+(idx_k-1)*L:idx_k*L,:) = Z_L_before*(Y_d_before^(idx_k-1));
    end
    for idx_k = 1:M-K+1
        E_R_before(:,1+(idx_k-1)*(N-L+1):idx_k*(N-L+1)) = (Y_d_before^(idx_k-1)) * Z_R_before;
    end   
    % shuffling matrix generation
    P = zeros(K*L,K*L);
    for idx_p = 1:K*L
        P(idx_p,K*L-idx_p+1) = 1;
    end
    X_ee = cat(2,X_en,P*conj(X_en));
    [U,S,V]= svd(X_ee);
    Us = U(:,1:I); %(KL,I)
    Un = U(:,I+1:K*L);%(KL,min-I)
    U1 = Us(1:(K-1)*L,:);
    U2 = Us(L+1:K*L,:);
    lambda_y = eig(U1'*U2,U1'*U1); % yi extraction
    y= mod(angle(lambda_y)/2/pi+1,1); % % f1i extraction
    P =zeros(K*L,K*L);
    for i = 1:K*L
        P(i, ceil(i/K)+mod(i-1,K)*L ) = 1;
    end
    Usp = P*Us;%(KL,I)
    U1p = Usp(1:(L-1)*K,:);
    U2p = Usp(K+1:K*L,:);
    lambda_z = eig(U1p'*U2p,U1p'*U1p);
    z= mod(angle(lambda_z)/2/pi+1,1);
    %pairing
    YL = zeros(K,I);
    ZL = zeros(K,I);
    for i = 1:K
        YL(i,:) = lambda_y.^(i-1);
    end
    for i = 1:K
        ZL(i,:) = lambda_z.^(i-1);
    end
    temp= kron(YL(:,1),ZL(:,1)); % (L*L,1)
    % pairing
    z_flag = zeros(1,I); % indicate select_z is selected or not.
    point_yz =zeros(1,I); % point_yz[y] = z
    for idx_y = 1:I
        yi = YL(:,idx_y);
        base_value = -Inf;
        J = zeros(1,I);
        for select_z = 1:I
            if z_flag(select_z)==1
                continue;
            else
                zi = ZL(:,select_z);
                eL = kron(yi,zi); %(KL,1);
                for un_idx = 1:I
                    J(select_z) = J(select_z) + (abs(Us(:,un_idx)'*eL))^2;
                end
            end

        end
        [m,i] = max(J);
        point_yz(idx_y) = i;
        z_flag(point_yz(idx_y)) = 1;
    end
    scatter(y(1),z(point_yz(1)),'k');
    scatter(y(2),z(point_yz(2)),'k');
    scatter(y(3),z(point_yz(3)),'k');
    Y_d_after=diag([exp(1j*2*pi*y(1)),exp(1j*2*pi*y(2)),exp(1j*2*pi*y(3))]);% I x I
    zd=exp(1j*2*pi*[z(point_yz(1)),z(point_yz(2)),z(point_yz(3))]);
    Z_L_after = zeros(L,I);
    Z_R_after = zeros(I,N-L+1);
    for idx_L = 1:L
        Z_L_after(idx_L,:) = zd.^(idx_L -1);
    end
    for idx_R = 1:N-L+1
        Z_R_after(:,idx_R) = (zd.').^(idx_R-1);
    end
    E_L_after = zeros(K*L,I);
    E_R_after = zeros(I,(N-L+1)*(M-K+1));
    for idx_k = 1:K
        E_L_after(1+(idx_k-1)*L:idx_k*L,:) = Z_L_after*(Y_d_after^(idx_k-1));
    end
    for idx_k = 1:M-K+1
        E_R_after(:,1+(idx_k-1)*(N-L+1):idx_k*(N-L+1)) = (Y_d_after^(idx_k-1)) * Z_R_after;
    end   
    E_L_Moore = (ctranspose(E_L_after)*E_L_after)^(-1)*ctranspose(E_L_after);
    E_R_Moore = (ctranspose(E_R_after))*(E_R_after*ctranspose(E_R_after))^(-1);
    A_estimate = E_L_Moore * X_en * E_R_Moore 
end

这里我设定的
$$\{\begin{array}{c}{f}_1=[0.1,0.2,0.3]\\ f_2=[0.4,0.5,0.6]\\ a=[1+1j,3+4j,6+2j]\end{array}⟶\{\begin{array}{c}y=[0.3,0.2,0.1]\\ z=[0.6,0.5,0.4]\\ a=\mathrm{diag}\left({\mathbf{E}}_L^{+}{X}_e{\mathbf{E}}_R^{+}\right)\end{array}$$

个人总结

matrix pencil矩阵铅笔算法是一种角度估计的高精度算法，在2D频率估计的过程中，使用一系列矩阵分解和变换操作使得需要估计的频率参数可以被映射为二维采样矩阵的范德蒙德域的基向量空间(不知道术语是否贴切)上，需要估计的参数就是各基向量的广义特征值！

---

1. Y. Hua, “Estimating two-dimensional frequencies by matrix enhancement and matrix pencil,” in IEEE Transactions on Signal Processing, vol. 40, no. 9, pp. 2267-2280, Sept. 1992, doi: 10.1109/78.157226. ↩︎