python实现 线性卷积用Toeplitz 矩阵运算

python实现 线性卷积用Toeplitz 矩阵运算

前言

在看论文的时候，发现Toeplitz 矩阵和线性卷积有关系，于是翻了程佩青老师的数字信号处理课本，发现是有讲过这点的。

Toeplitz 矩阵：从左上到右下的斜对角线都相同，如下

引子

以这道题为例子编程

python代码

一维情况

import numpy as np
from scipy.linalg import toeplitz

a=[3,7,5,-1,2]
b=[4,-1,2,3]
c_len=len(a)+len(b)-1

#把向量b变成Toeplitz 矩阵
col=[b[0] if i==0 else 0 for i in range(len(a))]
row=b+[0]*(c_len-len(b))
t = toeplitz(col,row)

c=np.dot(np.array(a),t)
print(c)

二维情况

import numpy as np

def convolution_to_matmul_2d(a,b,mode='same'):
    ''' 
    a: 1*n 
    b: m*k
    Toeplitz is (n+m-1)*m ,which is transforms from a.
    so a convolved with b becomes Toeplitz multiplied by b
    
    比如在地震反演中，反射系数与子波相卷积，子波是一维的雷克子波a，反射系数是二维的矩阵b
    exp:
    a = np.array([4,-1, 2, 3])
    b = np.array([[3,7,5,-1,2],
                [3,7,5,-1,2],
                [3,7,5,-1,2]]).T
    toeplitz_matrix.T: array([[ 4., -1.,  2.,  3.,  0.,  0.,  0.,  0.]
,
                              [ 0.,  4., -1.,  2.,  3.,  0.,  0.,  0.],
                              [ 0.,  0.,  4., -1.,  2.,  3.,  0.,  0.],
                              [ 0.,  0.,  0.,  4., -1.,  2.,  3.,  0.],
                              [ 0.,  0.,  0.,  0.,  4., -1.,  2.,  3.]])
    convolution_result: array([[12., 12., 12.],
                               [25., 25., 25.],
                               [19., 19., 19.],
                               [14., 14., 14.],
                               [40., 40., 40.],
                               [11., 11., 11.],
                               [ 1.,  1.,  1.],
                               [ 6.,  6.,  6.]])
    '''
    n=a.shape[0]
    m=b.shape[0]
    k=b.shape[1]

    ic(n,m)

    # 创建 Toeplitz 矩阵
    toeplitz_matrix=np.zeros((n+m-1,m))
    for i in range(len(b)):
        toeplitz_matrix[i:i+len(a), i] = a

    convolution_result=np.matmul(toeplitz_matrix,b)
    if mode =="same":
        convolution_result=convolution_result[n//2:n//2+m,:]
        # convolution_result=convolution_result[n//2:n//2+m,:]
    
    ic(convolution_result.shape)

    return convolution_result

record=convolution_to_matmul_2d(ricker_wavelet,reflection)
plt.imshow(record,cmap=plt.cm.seismic)
plt.title("record")

原理解释

书上P18把这个线性卷积如何变成矩阵推导的很详细，就不赘述了。最后可见H矩阵是依次循环右移一位后下移一行，使得形成各对角线相同的矩阵，即Toeplitz矩阵。

依次循环右移一位后下移一行，使得形成各对角线相同的矩阵，即Toeplitz矩阵。