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矩阵。