Matrix Pencil Method

将一个复时间序列写为有限个指数求和的形式，指数中的参数$s_i$为复数：
$$E(t)=\sum _i{A}_i{e}^{s_{i}t}$$

而且，各个组分的贡献是按从大到小依次排布的。

如果按照Z变换（$\simeq$离散拉普拉斯变换）的视角，Z[E(t)]=E(s), 则上式中的$s_i$其实就对应于E(s)的一阶极点，而$A_i$与留数相关。

Matrix Pencil Method原则上是一种PCA的方法, 针对离散的数据$E_1,E_2,...$。具有最小二乘意义下的最优。

我们参考了这一篇文章的方法:
Tapan K. Sakar and Odilon Pereira, Using the Matrix Pencil Method to Estimate the Parameters of Sum of Complex Exponetials, IEEE Antennas and Propagation Magazine, Vol. 37, No. 1, February 1995.

例子

拟合三个e指数组成的带噪声数据(后面代码的默认例子)：

t在0-1区间，信号的$A_i$和$s_i$分别是：
$A_0=1,A_1=6+3i,A_2=0.2i,$
$s_0=-5.2-12i,s_1=-3+23J,s_2=-3+90i$,
数据伴随5%左右的误差：$E(t)={\sum }_i^{}{A}_i{e}^{s_{i}t}+noise$
拟合效果是：

拟合一个带噪声的指数衰减振荡曲线：

Fit damping curves of an oscillator (raw data with ~5% amplitude noise and 10degree phase uncertainty, raw-blue, fit-red):

程序

python3 程序

'''
Matrix Pencil Method

@author: XXX
@date: 20190811
@updates:
@ref: Tapan K. Sakar and Odilon Pereira, Using the Matrix Pencil Method to Estimate the Parameters of Sum of Complex Exponetials, 
IEEE Antennas and Propagation Magazine, Vol. 37, No. 1, February 1995.
'''

import numpy as np 
import matplotlib.pyplot as plt 

def constructY(y, N, L):
    '''
    y: complex signal sequence.
    N: len(y)
    L: L<N, pencil parameter, N/3 < L < N/2 recommended. 
    return: constructed Y matrix, 
    [
        [y[0], y[1], ..., y[L-1]],
        [y[1], y[2], ..., y[L]], 
        ...
        [y[N-L-1], y[N-L], ..., y[N-1]]
    ]
    (N-L)*(L+1) matrix. 
    '''
    Y = np.zeros((N-L, L+1), dtype=np.complex_)
    for k in range(N-L):
        Y[k, :] = y[k:(k+L+1)]
    return Y 

def constructZM(Z, N):
    '''
    Z: 1-D complex array.
    return N*M complex matrix (M=len(Z)):
    [
        [1,  1,  1, ..., 1 ],
        [z[0], z[1], .., z[M-1]],
        ...
        [z[0]**(N-1), z[1]**(N-1), ..., z[M-1]**(N-1)]
    ]
    '''
    M = len(Z)
    ZM = np.zeros( (N, M), dtype=np.complex_) 
    for k in range(N):
        ZM[k, :] = Z**k 
    return ZM 

def SVDFilter(Sp, p=3.0):
    '''
    Sp: 1-D normed eigenvalues of Y after SVD, 1-st the biggest
    p: precise ditigits, default 3.0. 
    return: M, M is the first integer that S[M]/S_max <= 10**(-p)
    '''
    Sm = np.max(Sp) 
    pp = 10.0**(-p)
    for m in range(len(Sp)):
        if Sp[m]/Sm <= pp:
            return m+1 
    return m+1 


def pencil(N, y, dt, M=None, p=8.0, Lfactor=0.40):
    '''
    Purpose:
      Complex exponential fit of a sampled complex waveform by Matrix Pencil Method.
    Authors: 
      XXX
    Arguments:
      N    - number of data samples. ==len(y)       [INPUT]
      y    - 1-D complex array of sampled signal.   [INPUT]
      dt   - sample interval.                       [INPUT]
      M    - pencil parameter. 
             if None: use p to determin M.
             if given in range(0, Lfractor*N), then use it
             if given out of range, then use p to determin M.
      p    - precise digits of signal, default 8.0 or 10**(-8.0).
    Returns: (Z, R, M, (residuals, rank, s))
      Z    - 1-D Z array. 
      R    - 1-D R array.
      M    - M in use. 
      (residuals, rank, s)   - np.linalg.lstsq further results. 
    Method:
      y[k] = y(k*dt) = sum{i=0--> M} R[i]*( Z[i]**k ) 
      Z[i] = exp(si*dt)
    '''
    # better between N/3~N/2, pencil parameter:
    L = int(N*Lfactor)  
    # construct Y matrix (Hankel data matrix) from signal y[i], shape=(N-L, L+1):
    Y = constructY(y, N, L)
    # SVD of Y: 
    U, S, V = np.linalg.svd(Y, full_matrices=True)
    #results: U.shape=(N-L, N-L), S.shape=(L+1, ), V.shape=(L+1, L+1)

    # find M: 
    if M is None:
        M = SVDFilter(np.abs(S), p=p)
    elif M not in range(0, L+1):
        M = SVDFilter(np.abs(S), p=p) 
    else: 
        pass
    # matrix primes based on M:
    #Vprime = V[0:M, :] # remove [M:] data set. only 0, 1, 2, ..., M-1 remains
    #Sprime = S[0:M]
    V1prime = V[0:M, 0:-1] # remove last column
    V2prime = V[0:M, 1:] # remove first column
    #smat = np.zeros((U.shape[0], M), dtype=np.complex_)
    #smat[:M, :M] = np.diag(Sprime)
    #Y1 = np.dot(U[:-1, :], np.dot(smat, V1prime))

    V1prime_H_MPinv = np.linalg.pinv(V1prime.T) # find V1'^+ , Moore-Penrose pseudoinverse 
    V1V2 = np.dot(V1prime_H_MPinv, V2prime.T) # construct V1V2 = np.dot(V1'^+, V2') 
    Z = np.linalg.eigvals(V1V2) # find eigenvalues of V1V2. Zs.shape=(M,)
    #print(V1V2.shape, Z)

    # find R by solving least-square problem: Y = np.dot(ZM, R)
    ZM = np.row_stack([Z**k for k in range(N)]) # N*M 
    R, residuals, rank, s = np.linalg.lstsq(ZM, y)
    return (Z, R, M, (residuals, rank, s))


if __name__ == '__main__':
    # testing
    dt = 0.002
    x = np.arange(0.0, 1.0, dt,  dtype=np.float_) 
    N = len(x) # sample number 

    import random 
    # sampled function
    y = (np.exp(-x*(5.2+12.0j)) + \
        (6.0+3J)*np.exp(-x*(3. - 23.0J)) + \
        (0.20J)*np.exp(x*(3. - 90.0J)) )*(
            1.0 + 0.05*(np.random.randn(N,)+1J*np.random.randn(N,))
            )

    # test here: 
    # p could be determined approximately by SNR level. 
    Z, R, M, _ =  pencil(N, y, dt, M=None, p=-np.log10(0.05), Lfactor=0.33)
    print(np.log(Z)/dt)
    print(R)
    print('M=', M) 
    
    svals = np.log(Z)/dt 
    avals = R 
    ypred = x*0. + 0J 
    
    for si, ai in zip(svals, avals):
        ypred += ai*np.exp(si*x) 
    
    fig = plt.figure() 
    ax = fig.add_subplot(111) 
    ax.set_xlabel("t") 
    ax.set_ylabel("y(t)") 
    ax.plot(x, np.real(ypred), "r-", label="fit", lw=1.5)  
    ax.scatter(x, np.real(y), 6, label="raw", color='b')  
    ax.legend()
    plt.show()
    '''
    L = int(N/3) # better between N/3~N/2, pencil parameter. 
    Y = constructY(y, N, L)
    U, S, V = np.linalg.svd(Y, full_matrices=True)
    print(Y.shape, '=', U.shape, S.shape, V.shape)
    print(S)
    M = SVDFilter(np.abs(S), p=10.0) 
    print(M)

    #smat = np.zeros((U.shape[0], V.shape[0]), dtype=np.complex_)
    #smat[:V.shape[0], :V.shape[0]] = np.diag(S)
    #Yp = np.dot(U, np.dot(smat, V))
    #print(np.allclose(Y, Yp) )
    
    # matrix primes based on M:
    Vprime = V[0:(M), :] 
    Sprime = S[0:(M)]
    V1prime = Vprime[:, 0:-1]
    V2prime = Vprime[:, 1:]
    
    #smat = np.zeros((U.shape[0], M), dtype=np.complex_)
    #smat[:M, :M] = np.diag(Sprime)
    #Y1 = np.dot(U[:-1, :], np.dot(smat, V1prime))

    V1prime_H_MPinv = np.linalg.pinv(V1prime.T) 
    V1V2 = np.dot(V1prime_H_MPinv, V2prime.T ) 
    Zs = np.linalg.eigvals(V1V2)
    print(V1V2.shape, Zs)
    '''

程序可能的输出是：

[-11.16113418-766.10559942j  -3.01301208 +23.02930587j
  -5.26331942 -12.41123864j   2.98811039 -89.99467627j]
[-3.26109356e-02-0.15529529j  6.03994228e+00+2.92902757j
  9.86343936e-01+0.06872874j  1.96607132e-03+0.20251491j]
M= 4

其他应用

一方面，可以由系统响应数据找出系统的极点。
另一方面，可以作傍轴光束线性传播行为的分析(模式分解为一系列不同大小的高斯模)：
傍轴柱对称的光束，可以利用此方法，把z=0处的横场模式写成几个高斯型的求和的形式，从而得到任意z处的横场分布。 E(x), 将X=x^2视为新的x, E(X)可以写为有限e指数和的形式，
而每一个成分$A_i\exp (s_{i}X)\equiv A_i\exp (s_i{x}^2)$ 就是高斯光束的场分布。

相关话题

stackexchange:1428566/fit-sum-of-exponentials
mexpfit: Multi-exponential Fitting
C++ modified fast ESPRIT algorithm
一些文章：

• Fast ESPRIT algorithms based on partial singular value decompositions
• Coding Prony’s method in MATLAB and applying it to biomedical signal filtering
• Fast fitting of multi-exponential decay curves
• Modern Characterization of Electromagnetic Systems and Its Associated Metrology - 2021 - Sarkar - Matrix Pencil Method MPM
• mxpfit- A library for finding optimal multi-exponential approximations
• On approximation of functions by exponential sums
• Matrix Pencil Method for Estimating Parameters of Exponentially Damped/Undamped Sinusoids in Noise