Lossy cumsum в numy - PullRequest
       1

Lossy cumsum в numy

0 голосов
/ 11 июля 2019

У меня есть массив a длины N и мне нужно выполнить следующую операцию:

b[n] = sum_{i=0}^{n} a[i]p^{n-i}

С p в [0..1].Это уравнение представляет собой сумму с потерями, где первые индексы в сумме взвешиваются с большей потерей (p ^ {ni}), чем последние.Последний индекс (i = n) всегда взвешивается на 1. Если p = 1, то операция является простой суммой.

b = np.cumsum(a)

Если, если p! = 1, я могу реализовать эту операцию вНеэффективный для процессора способ:

b = np.empty(np.shape(a))
# I'm using the (-1,-1,-1) idiom for reversed ranges
p_vec = np.power(p, np.arange(N-1, 0-1, -1))
# p_vec[0] = p^{N-1}, p_vec[-1] = 1

for n in range(N):
   b[n] = np.sum(a[:n+1]*p_vec[-(n+1):])

Или неэффективным с точки зрения памяти, но векторизованным способом (ИМО тоже неэффективен с точки зрения использования процессора, поскольку тратится много времени):

a_idx = np.reshape(np.arange(N+1), (1, N+1)) - np.reshape(np.arange(N-1, 0-1, -1), (N, 1))
a_idx = np.maximum(0, a_idx)

# For N=4, a_idx looks like this:
# [[0, 0, 0, 0, 1],
#  [0, 0, 0, 1, 2],
#  [0, 0, 1, 2, 3],
#  [0, 1, 2, 3, 4]]

a_ext = np.concatenate(([0], a,), axis=0) # len(a_ext) = N + 1
p_vec = np.power(p, np.arange(N, 0-1, -1)) # len(p_vec) = N + 1

b = np.dot(a_ext[a_idx], p_vec)

Есть лилучший способ достичь этого «с потерями» cumum?

1 Ответ

1 голос
/ 11 июля 2019

То, что вы хотите, это БИХ-фильтр, вы можете использовать scipy.signal.lfilter(), вот код:

Ваш код:

import numpy as np
N = 10
p = 0.8
np.random.seed(0)
x = np.random.randn(N)
y = np.empty_like(x)
p_vec = np.power(p, np.arange(N-1, 0-1, -1))
for n in range(N):
    y[n] = np.sum(x[:n+1]*p_vec[-(n+1):])
y    

выход:

array([1.76405235, 1.81139909, 2.42785725, 4.183179  , 5.21410119,
       3.19400307, 3.50529088, 2.65287549, 2.01908154, 2.02586374])

При использовании lfilter():

from scipy import signal
y = signal.lfilter([1], [1, -p], x)
print(y)

вывод:

array([1.76405235, 1.81139909, 2.42785725, 4.183179  , 5.21410119,
       3.19400307, 3.50529088, 2.65287549, 2.01908154, 2.02586374])
Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...