Python: 有效地计算 Cadzow 过滤器的非对角线元素的平均值
Python: Efficienctly calculate mean of off-diagonal elements for Cadzow filter
我目前在 Python 中实施了 Gadzow 过滤器。
放在一些上下文中。您从一个一维数组开始(让我们以 range(10) 为例)并从中构建一个类似 Hankel 的矩阵,例如
H= [[0, 1, 2, 3, 4, 5],
[1, 2, 3, 4, 5, 6],
[2, 3, 4, 5, 6, 7],
[3, 4, 5, 6, 7, 8],
[4, 5, 6, 7, 8, 9]])
然后你用这个矩阵做一些线性代数,这没问题。之后,最耗时的步骤是平均问题。
在新矩阵 B 中,您对结果矩阵的元素进行平均。在第一行中,您通过 H 中的准确度给出的路径对所有元素进行平均。所以类似于非对角线,但从右上角到左下角。在第二个切片中,您忽略第一行,依此类推。
矩阵 $H$ 在此分析步骤下是不变的,但例如矩阵
1 2 2 1
1 1 1 1
1 1 1 1
会变成
1 1.5 1.33 1
1 1 1 1
1 1 1 1
好的,我希望你能理解这个问题。我的(工作但效率低下)代码是
def av_diag(A,i,j):
dim = A.shape
# get the "borders" of A
lim = min((dim[0]-i,j+1))
# calculate the mean
return np.mean([A[i+it,j-it] for it in range(lim)])
def avHankel(A):
# get the mean for all elements by nested list comprehension
return np.array([[av_diag(A,i,j) for j in range(len(A[0]))] for i in range(len(A))])
我的数据需要一些时间,包含 2048 个数据点,生成 1024x1023 矩阵。
我很高兴有可能的技巧来加快速度。
谢谢
您可以将输入矩阵与 filter 矩阵进行卷积以加速您的代码。可以定义 filter 矩阵,以便在卷积的每个步骤中,它仅提取给定坐标处的反对角线。基本上,您的 filter 矩阵只是一个反恒等式矩阵。最后,由于卷积只会对反对角线的元素求和,您必须将输出除以正确的样本数以获得均值:
import numpy as np
from scipy.signal import fftconvolve
from time import time
def av_diag(A,i,j):
dim = A.shape
lim = min((dim[0]-i,j+1))
return np.mean([A[i+it,j-it] for it in range(lim)])
def avHankel(A):
return np.array([[av_diag(A,i,j) for j in range(len(A[0]))] for i in range(len(A))])
def fast_avHankel(A):
m, n = A.shape
filt = np.eye(m)[:,::-1]
Apad = np.pad(A, ((0, m-1), (m-1, 0)), mode = "constant", constant_values = 0)
Asum = fftconvolve(Apad, filt, mode = "valid")
Adiv = np.array([ [ min(m-i, j+1) for j in range(n) ] for i in range(m) ])
return Asum / Adiv
if __name__ == "__main__":
A = np.random.rand(500, 500)
starttime = time()
Hold = avHankel(A)
print(time() - starttime) # 10.6 seconds on a laptop
starttime = time()
Hnew = fast_avHankel(A)
print(time() - starttime) # 0.26 seconds on a laptop
我目前在 Python 中实施了 Gadzow 过滤器。
放在一些上下文中。您从一个一维数组开始(让我们以 range(10) 为例)并从中构建一个类似 Hankel 的矩阵,例如
H= [[0, 1, 2, 3, 4, 5],
[1, 2, 3, 4, 5, 6],
[2, 3, 4, 5, 6, 7],
[3, 4, 5, 6, 7, 8],
[4, 5, 6, 7, 8, 9]])
然后你用这个矩阵做一些线性代数,这没问题。之后,最耗时的步骤是平均问题。
在新矩阵 B 中,您对结果矩阵的元素进行平均。在第一行中,您通过 H 中的准确度给出的路径对所有元素进行平均。所以类似于非对角线,但从右上角到左下角。在第二个切片中,您忽略第一行,依此类推。
矩阵 $H$ 在此分析步骤下是不变的,但例如矩阵
1 2 2 1
1 1 1 1
1 1 1 1
会变成
1 1.5 1.33 1
1 1 1 1
1 1 1 1
好的,我希望你能理解这个问题。我的(工作但效率低下)代码是
def av_diag(A,i,j):
dim = A.shape
# get the "borders" of A
lim = min((dim[0]-i,j+1))
# calculate the mean
return np.mean([A[i+it,j-it] for it in range(lim)])
def avHankel(A):
# get the mean for all elements by nested list comprehension
return np.array([[av_diag(A,i,j) for j in range(len(A[0]))] for i in range(len(A))])
我的数据需要一些时间,包含 2048 个数据点,生成 1024x1023 矩阵。
我很高兴有可能的技巧来加快速度。
谢谢
您可以将输入矩阵与 filter 矩阵进行卷积以加速您的代码。可以定义 filter 矩阵,以便在卷积的每个步骤中,它仅提取给定坐标处的反对角线。基本上,您的 filter 矩阵只是一个反恒等式矩阵。最后,由于卷积只会对反对角线的元素求和,您必须将输出除以正确的样本数以获得均值:
import numpy as np
from scipy.signal import fftconvolve
from time import time
def av_diag(A,i,j):
dim = A.shape
lim = min((dim[0]-i,j+1))
return np.mean([A[i+it,j-it] for it in range(lim)])
def avHankel(A):
return np.array([[av_diag(A,i,j) for j in range(len(A[0]))] for i in range(len(A))])
def fast_avHankel(A):
m, n = A.shape
filt = np.eye(m)[:,::-1]
Apad = np.pad(A, ((0, m-1), (m-1, 0)), mode = "constant", constant_values = 0)
Asum = fftconvolve(Apad, filt, mode = "valid")
Adiv = np.array([ [ min(m-i, j+1) for j in range(n) ] for i in range(m) ])
return Asum / Adiv
if __name__ == "__main__":
A = np.random.rand(500, 500)
starttime = time()
Hold = avHankel(A)
print(time() - starttime) # 10.6 seconds on a laptop
starttime = time()
Hnew = fast_avHankel(A)
print(time() - starttime) # 0.26 seconds on a laptop