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$有一种更加增强后的矩阵形式——
X e e = [ X e , P e X e ∗ ] P e = [ 1 1 ⋯ 1 ] 因为 r a n g e ( X e e ) = r a n g e ( X e ) = r a n g e ( E L ) \boldsymbol{X}_{ee}=\left[ \boldsymbol{X}_e,\boldsymbol{P}_e\boldsymbol{X}_{e}^{*} \right] \\ \boldsymbol{P}_e=\left[ \begin{matrix} & & & 1\\ & & 1& \\ & \cdots& & \\ 1& & & \\ \end{matrix} \right] \\ 因为range\left( \boldsymbol{X}_{ee} \right) =range\left( \boldsymbol{X}_e \right) =range\left( \boldsymbol{E}_L \right) Xee​=[Xe​,Pe​Xe∗​]Pe​= ​1​⋯​1​1​ ​因为range(Xee​)=range(Xe​)=range(EL​)

实验仿真

我们考虑这样的信号模型
$$\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:

其他问题

估计 A \boldsymbol{A} 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）：
X 是 A 的违逆，如果 X = V Ξ U H 其中酉矩阵 U , V 来自对 A 的 S V D 分解 A = U Σ V H 设 Σ = [ σ 1 ⋱ σ r O ] , 则 Ξ = [ 1 σ 1 ⋱ 1 σ r O ] ( A 可逆时， O 为 0 阶阵 ) X\text{是}A\text{的违逆，如果} \\ X=V\varXi U^H \\ \text{其中酉矩阵}U,V\text{来自对}A\text{的}SVD\text{分解}A=U\varSigma V^H \\ \text{设}\varSigma =\left[ \begin{matrix} \sigma _1& & & \\ & \ddots& & \\ & & \sigma _r& \\ & & & O\\ \end{matrix} \right] , \text{则}\varXi =\left[ \begin{matrix} \frac{1}{\sigma _1}& & & \\ & \ddots& & \\ & & \frac{1}{\sigma _r}& \\ & & & O\\ \end{matrix} \right] \left( A\text{可逆时，}O\text{为}0\text{阶阵} \right) X是A的违逆，如果X=VΞUH其中酉矩阵U,V来自对A的SVD分解A=UΣVH设Σ= ​σ1​​⋱​σr​​O​ ​,则Ξ= ​σ1​1​​⋱​σr​1​​O​ ​(A可逆时，O为0阶阵)
上述违逆内容来自于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逆。
设 r = r a n k ( A m × n ) 将 A 分解为 A m × n = P m × m [ I r O O O ] m × n Q n × n , 其中 P , Q 可逆 , I r 表示 r 阶单位阵 , A 可逆时 , O 为 0 阶阵 令 R m × r = P m × m [ I r O ] m × r , S = [ I r    O ] r × n Q n × n 则 X n × m = S H n × r ( R H r × m A m × n S H n × r ) − 1 r × r R H r × m 将 A 做 S V D ，有 A = U Σ V H Σ = [ σ 1 ⋱ σ r O ] , 则 Ξ = [ 1 σ 1 ⋱ 1 σ r O ] ( A 可逆时， O 为 0 阶阵 ) 现在证明一下 X n × m = S H n × r ( R H r × m A m × n S H n × r ) − 1 r × r R H r × m = V Ξ U H 这个等式成立 ( M o o r e − P e n r o s e 逆 = 违逆 ) p r o v e : 我们的目标是证明 X n × m = S H ( R H A S H ) − 1 R H = V Ξ U H 即可 左边 = S H ( R H U Σ V H S H ) − 1 R H = S H ( S H ) − 1 ( V H ) − 1 Ξ U − 1 ( R H ) − 1 R H = ( V H ) − 1 Ξ U − 1 = V Ξ U H ( 最后一步用到酉矩阵的性质 V H = V − 1 , U H = U H ) \text{设}r=rank\left( A_{m\times n} \right) \\ \text{将}A\text{分解为}A_{m\times n}=P_{m\times m}\left[ \begin{matrix} I_r& O\\ O& O\\ \end{matrix} \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\varSigma V^H \\ \varSigma =\left[ \begin{matrix} \sigma _1& & & \\ & \ddots& & \\ & & \sigma _r& \\ & & & O\\ \end{matrix} \right] , \text{则}\varXi =\left[ \begin{matrix} \frac{1}{\sigma _1}& & & \\ & \ddots& & \\ & & \frac{1}{\sigma _r}& \\ & & & O\\ \end{matrix} \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\varXi U^H\text{这个等式成立}\left( Moore-Penrose\text{逆}=\text{违逆} \right) \\ prove:\text{我们的目标是证明}X_{n\times m}=S^H\left( R^HAS^H \right) ^{-1}R^H=V\varXi U^H\text{即可} \\ \text{左边}=S^H\left( R^HU\varSigma V^HS^H \right) ^{-1}R^H=S^H\left( S^H \right) ^{-1}\left( V^H \right) ^{-1}\varXi U^{-1}\left( R^H \right) ^{-1}R^H=\left( V^H \right) ^{-1}\varXi U^{-1}=V\varXi U^H\left( \text{最后一步用到酉矩阵的性质}V^H=V^{-1},U^H=U^H \right) 设r=rank(Am×n​)将A分解为Am×n​=Pm×m​[Ir​O​OO​]m×n​Qn×n​,其中P,Q可逆,Ir​表示r阶单位阵,A可逆时,O为0阶阵令Rm×r​=Pm×m​[Ir​O​]m×r​,S=[Ir​O]r×n​Qn×n​则Xn×m​=SHn×r​(RHr×m​Am×n​SHn×r​)−1r×r​RHr×m​将A做SVD，有A=UΣVHΣ= ​σ1​​⋱​σr​​O​ ​,则Ξ= ​σ1​1​​⋱​σr​1​​O​ ​(A可逆时，O为0阶阵)现在证明一下Xn×m​=SHn×r​(RHr×m​Am×n​SHn×r​)−1r×r​RHr×m​=VΞUH这个等式成立(Moore−Penrose逆=违逆)prove:我们的目标是证明Xn×m​=SH(RHASH)−1RH=VΞUH即可左边=SH(RHUΣVHSH)−1RH=SH(SH)−1(VH)−1ΞU−1(RH)−1RH=(VH)−1ΞU−1=VΞUH(最后一步用到酉矩阵的性质VH=V−1,UH=UH)
违逆和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被估计出来，下面两个矩阵可以得到
Z L = [    z 1 0    z 2 0 ⋯ z I 0   z 1    z 2 ⋯ z I    ⋮    ⋮ ⋯    ⋮ z 1 L − 1 z 2 L − 1 ⋯ z I    L − 1 ] L × I Z R = [    z 1 0    z 1 ⋯ z 1 N − L   z 2 0    z 2 ⋯ z 2 N − L    ⋮    ⋮ ⋯    ⋮ z I 0 z I ⋯ z I    N − L ] I × N − L + 1 \boldsymbol{Z}_L=\left[ \begin{matrix} \,\,z_{1}^{0}& \,\,z_{2}^{0}& \cdots& z_{I}^{0}\\ \,z_1& \,\,z_2& \cdots& z_I\\ \,\,\vdots& \,\,\vdots& \cdots& \,\,\vdots\\ z_{1}^{L-1}& z_{2}^{L-1}& \cdots& z_{I\,\,}^{L-1}\\ \end{matrix} \right] _{L\times I} \\ \boldsymbol{Z}_R=\left[ \begin{matrix} \,\,z_{1}^{0}& \,\,z_1& \cdots& z_{1}^{N-L}\\ \,z_{2}^{0}& \,\,z_2& \cdots& z_{2}^{N-L}\\ \,\,\vdots& \,\,\vdots& \cdots& \,\,\vdots\\ z_{I}^{0}& z_I& \cdots& z_{I\,\,}^{N-L}\\ \end{matrix} \right] _{I\times N-L+1} ZL​= ​z10​z1​⋮z1L−1​​z20​z2​⋮z2L−1​​⋯⋯⋯⋯​zI0​zI​⋮zIL−1​​ ​L×I​ZR​= ​z10​z20​⋮zI0​​z1​z2​⋮zI​​⋯⋯⋯⋯​z1N−L​z2N−L​⋮zIN−L​​ ​I×N−L+1​

X e = E L A E R E L = [ Z L Z L Y d ⋯ Z L Y d K − 1 ] K L × I E R = [ Z R , Y d Z R , ⋯   , Y d M − K Z R ] I × ( N − L + 1 ) ( M − K + 1 ) \boldsymbol{X}_e=\boldsymbol{E}_L\boldsymbol{AE}_R \\ \boldsymbol{E}_L=\left[ \begin{array}{l} \boldsymbol{Z}_L\\ \boldsymbol{Z}_L\boldsymbol{Y}_d\\ \cdots\\ \boldsymbol{Z}_L\boldsymbol{Y}_{d}^{K-1}\\ \end{array} \right] _{KL\times I} \\ \boldsymbol{E}_R=\left[ \boldsymbol{Z}_R,\boldsymbol{Y}_d\boldsymbol{Z}_R,\cdots ,\boldsymbol{Y}_{d}^{M-K}\boldsymbol{Z}_R \right] _{I\times \left( N-L+1 \right) \left( M-K+1 \right)} Xe​=EL​AER​EL​= ​ZL​ZL​Yd​⋯ZL​YdK−1​​ ​KL×I​ER​=[ZR​,Yd​ZR​,⋯,YdM−K​ZR​]I×(N−L+1)(M−K+1)​
从${\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