Две несложные головоломки с временными рядами.

tengiz · Post by **tengiz** » 17 Jan 2002 07:18

Имеем вещественный временной ряд: X(n) для n = 0..N-1, 0 для n < 0.

Задача 1. простая:

нужно быстро посчитать взвешенное среднее вещественного временного ряда X(n) для случая, когда окно W(k) = w для –K+1 < k < K-1, 0 для всех других случаев. N >> K.

Задача 2. немного сложнее:

нужно быстро посчитать взвешенное среднее вещественного временного ряда X(n), для случая, когда окно W(k) = w * (1 + cos (pi * k / K)) / 2 для –K+1 < k < K-1, 0 для всех других случаев.

N >> K. Что значит быстро: асимптотическая сложность должна быть заведомо лучше, чем O (N*K).

Vovka · Post by **Vovka** » 17 Jan 2002 13:13

tengiz,
можно поинтересоваться, что такое окно и как считается взвешенное среднее? И правильно ли я понял, что "временной ряд" это просто чистовая последовательность?

Или смысл этой задачи именно в знании статистики?

Melkor · Post by **Melkor** » 17 Jan 2002 14:08

<blockquote>quote:<hr>Originally posted by tengiz:
Имеем вещественный временной ряд: X(n) для n = 0..N-1, 0 для n < 0.

Задача 1. простая:

нужно быстро посчитать взвешенное среднее вещественного временного ряда X(n) для случая, когда окно W(k) = w для –K+1 < k < K-1, 0 для всех других случаев. N >> K.

Задача 2. немного сложнее:

нужно быстро посчитать взвешенное среднее вещественного временного ряда X(n), для случая, когда окно W(k) = w * (1 + cos (pi * k / K)) / 2 для –K+1 < k < K-1, 0 для всех других случаев.

N >> K. Что значит быстро: асимптотическая сложность должна быть заведомо лучше, чем O (N*K).<hr></blockquote>

Как я понял, речь идет о конволюции (иначе все нижеизложенное не в тему). Задачу 1 я знал, поэтому приводить ответ не буду. Задача 2, в принципе, решается по аналогии, так что особой чести ее решить, зная решение первой, нет, но идея такова: для окна храним суммы Nk, Nk*cos(k*pi/K), Nk*sin(k*pi/K), при переносе окна на следующий элемент вычитаем хвост добавляем голову (как в задаче 1) для всех 3х сумм, и пересчитываем суммы косинусов и синусов где-то так:
- new sum(cos) = sum(cos)*cos(pi/K) + sum(sin)*sin(pi/K)
- new sum(sin) = sum(sin)*cos(pi/K) - sum(cos)*sin(pi/K).
При этом надо достаточно аккуратно вести себя с головой и хвостом окна.

tengiz · Post by **tengiz** » 17 Jan 2002 19:21

Melkor, в Вашем решении есть ошибка.

tengiz · Post by **tengiz** » 17 Jan 2002 19:25

Vovka, взвешенное среднее - это новый временной ряд Y(n), каждый элемент которого равен сумме значений X(n+k) * W (k). В нашем случае суммирование делается по k (-K+1..K-1).

Melkor · Post by **Melkor** » 17 Jan 2002 20:35

<blockquote>quote:<hr>Originally posted by tengiz:
Melkor, в Вашем решении есть ошибка.<hr></blockquote>

Принципиальная, или типа того, что сначала надо вычесть хвост, потом обновить суммы, а потом добавить к ним голову, ну или знак где не тот?

tengiz · Post by **tengiz** » 17 Jan 2002 22:29

<blockquote>quote:<hr>Originally posted by Melkor:
...или типа того, что сначала надо вычесть хвост, потом обновить суммы, а потом добавить к ним голову, ну или знак где не тот?<hr></blockquote>Нет, дело не в этом. С этим вообще можно просто обойтйсь и сказать, что на вход алгоритма сразу поступает разность X(n+K-1) - X(n-K+1). Hint - откуда у Вас вообще появилось вращение (или умножение на комплексную экспоненту)? Какие суммы и как нужно обновлять - должно ли быть всё настолько симметрично?

Melkor · Post by **Melkor** » 17 Jan 2002 23:24

<blockquote>quote:<hr>Originally posted by tengiz:
Нет, дело не в этом. С этим вообще можно просто обойтйсь и сказать, что на вход алгоритма сразу поступает разность X(n+K-1) - X(n-K+1). Hint - откуда у Вас вообще появилось вращение (или умножение на комплексную экспоненту)? Какие суммы и как нужно обновлять - должно ли быть всё настолько симметрично?<hr></blockquote>

[img:7361d5059e]images/smiles/icon_sad.gif[/img:7361d5059e] Что-то до сих пор не пойму, где ошибка. Вращение - оттого, что функция под окном вращается (только края обрезаются). Вот мой ход решения более подробно:

Пусть для шага i имеются sum_cos(i) = sum(k=-K+1..K-1)(cos(k*pi/K)*N(i+k)), sum_sin(i) = sum(k=-K+1..K-1)(sin(k*pi/K)*N(i+k) и sum(i) = sum(k=-K+1..K-1)*N(i+k). Нам нужны значения sum_cos(i+1), sum_sin(i+1) и sum(i+1).

sum_cos(i+1) = sum(k=-K+1..K-1)(cos(k*pi/K)*N(i+1+k) =
sum(k=-K+2..K)(cos((k-1)*pi/K)*N(i+k) =
sum(k=-K+1..K-1)(cos((k-1)*pi/K)*N(i+k) - cos((-K+1)*pi/K)*N(i-K+2) + cos((K-1)*pi/K)*N(i+K)

sum(k=-K+1..K-1)(cos((k-1)*pi/K)*N(i+k) = sum(...)((cos(k*pi/K)*cos(pi/K) + sin(k*pi/K)*sin(pi/K))*N(i+k)) = sum_cos(i)*cos(pi/K) + sum_sin(i)*sin(pi/K).

Аналогично считаем sum_sin(i+1). Y(i) = w/2*(sum(i) + sum_cos(i))

Drom · Post by **Drom** » 17 Jan 2002 23:46

tengiz,
1) можно ссылку на определение терминов и какой-нубудь introduction для чайников типа меня? (что ето еще за конволюция, блин [img:1082ae4dfc]images/smiles/icon_sad.gif[/img:1082ae4dfc] )
2) сложность должна получиться О(N+K) (то есть =О(N) при N>>K) ?

tengiz · Post by **tengiz** » 18 Jan 2002 03:22

Drom:

1) Convolution - это свёртка. См., например, Discrete Convolution.
2) Да.

tengiz · Post by **tengiz** » 18 Jan 2002 03:33

Melkor, если хотите, могу Вам в привате написать, хотя Вам, на мой взгляд, всего лишь нужно внимательнее посмотреть на Ваше решение - и всё.

Hint #2: У вас есть вектор (x, 0) - как посчитать его X координату после поворота на f радиан?

Hint #3: Прогоните через Ваш алгоритм последовательность X(n), где X = 1 для n = K-1, и 0 во всех остальных случаях.

Melkor · Post by **Melkor** » 18 Jan 2002 13:52

<blockquote>quote:<hr>Originally posted by tengiz:
Melkor, если хотите, могу Вам в привате написать, хотя Вам, на мой взгляд, всего лишь нужно внимательнее посмотреть на Ваше решение - и всё.

Hint #2: У вас есть вектор (x, 0) - как посчитать его X координату после поворота на f радиан?

Hint #3: Прогоните через Ваш алгоритм последовательность X(n), где X = 1 для n = K-1, и 0 во всех остальных случаях.<hr></blockquote>

tengiz, Вы надо мной издеваетесь [img:f723263c20]images/smiles/icon_smile.gif[/img:f723263c20] Я написал даже программу, которая считает конволюцию грубой силой и итеративно сдвигая окно - разница в 16 знаке... (правда, я проверял для W(k) = cos(pi*k/K), а не W(k) = w * (1 + cos (pi * k / K)) / 2). Вот код:

<blockquote>code:<hr><pre>
const int N=30,K=5;
double X[N],Y[N];
int i,k;
srand(0);
for(i=0;i<N;i++)
X[i]=rand()*(1.0/RAND_MAX), Y[i]=0;

for(i=K-1;i<=N-K;i++)
for(k=-K+1;k<=K-1;k++) Y[i]+=cos(pi/K*k)*X[i+k];

double sumsin=0,sumcos=0,sumcosnew;
for(k=-K+1;k<=K-1;k++) {
sumcos += cos(pi/K*k)*X[K-1+k];
sumsin += sin(pi/K*k)*X[K-1+k];
}

for(i=K-1;i<=N-K;i++) {
Y[i] -= sumcos;
sumcosnew = sumcos*cos(pi/K) + sumsin*sin(pi/K) + X[i-K+1] + cos((K-1)*pi/K)*X[i+K];
sumsin = sumsin*cos(pi/K) - sumcos*sin(pi/K) + sin((K-1)*pi/K)*X[i+K];
sumcos = sumcosnew;
}
</pre><hr></blockquote>

tengiz · Post by **tengiz** » 18 Jan 2002 23:38

Melkor, да ну бросьте - это Вы надо мной издеваетесь [img:c215524330]images/smiles/icon_smile.gif[/img:c215524330] .

Вы же алгоритм пишете, а не формулу – поэтому пара лишних умножений и сложение do make difference (не говоря уже о вычислении sin и cos, хотя это уже на самом деле непринципиально). Принципиально здесь другое - Вы на каждом шаге итерации доворачиваете Ваш вектор (X(i), 0) на угол pi/K, хотя этого на самом деле не нужно делать. Способ избежать лишних вычислений - сдвинуть всё на один шаг. Т.е. Ваш последний цикл лучше смотрится следующим образом (конечно, граничные условия цикла при этом тоже нужно сдвинуть – я позволю себе не ломать голову над правильными сдвигами индекса в X[i]):<blockquote>code:<hr><pre>sumcosnew = sumcos * cos(pi/K) + sumsin * sin(pi/K) + X[i];
sumsin = sumsin * cos(pi/K) - sumcos * sin(pi/K);
sumcos = sumcosnew;</pre><hr></blockquote>При этом, как Вы сами видите, после поворота нужно обновлять только одну сумму, а не все. И вместо 6 умножений и 5 сложений, мы имеем 4 умножения и 3 сложения на шаг. Ну и на десерт – комплексное умножение можно ещё несколько оптимизировать и вместо 4 умножений и двух сложений сделать 3 умножения и 4 сложения, что имеет смысл на аппаратных платформах, где умножение заметно дороже сложения.

Melkor · Post by **Melkor** » 19 Jan 2002 00:24

<blockquote>quote:<hr>Originally posted by tengiz:

<blockquote>code:<hr><pre>sumcosnew = sumcos * cos(pi/K) + sumsin * sin(pi/K) + X[i];
sumsin = sumsin * cos(pi/K) - sumcos * sin(pi/K);
sumcos = sumcosnew;</pre><hr></blockquote><hr></blockquote>

Так лучше, конечно. Только я ж не алгоритм приводил, а так, общую концепцию... А за задачу спасибо - может, пригодится как-нибудь.

tengiz · Post by **tengiz** » 19 Jan 2002 03:52

<blockquote>quote:<hr>Originally posted by Melkor:
...может, пригодится как-нибудь.<hr></blockquote>Я это использовал при реализации FIR фильтра реального времени - можно каскадировать (или распараллеливать) такие алгоритмы для приближения конечного окна любой формы, просто разложив функцию окна в ряд Фурье. Так что пригодиться такая фигня действительно может.

Valeus · Post by **Valeus** » 12 Feb 2002 19:02

<blockquote>quote:<hr>Originally posted by tengiz:
Я это использовал при реализации FIR фильтра реального времени - можно каскадировать (или распараллеливать) такие алгоритмы для приближения конечного окна любой формы, просто разложив функцию окна в ряд Фурье. Так что пригодиться такая фигня действительно может.<hr></blockquote>

Все это очень хорошо, только на длинных рядах не будет ли потеря точности помехой (FP)?

tengiz · Post by **tengiz** » 12 Feb 2002 19:44

<blockquote>quote:<hr>Originally posted by Valeus:
Все это очень хорошо, только на длинных рядах не будет ли потеря точности помехой (FP)?<hr></blockquote>Да, конечно будет. Накапливание ошибки будет особенно заметно на DSP процессорах с 32 битной или более короткой плавающей точкой. Но на DSP процессоре прямое вычисление и так оказывается довольно быстрым. Фокус в том, что вся эта муть была использована на Intel при точности отсчётов временного ряда 16 бит, поэтому разрядность в 64 бита (много фильтров) или 80 бит (один фильтр - аккомуляторы фильра и все коэффициетны в регистрах FPU) вполне достаточна даже для длинных рядов.

tengiz · Post by **tengiz** » 12 Feb 2002 19:53

Valeus, спасибо за Ваш вопрос - похоже, есть простой способ избавиться даже от таких ошибок даже на бесконечных рядах - это требует несложной коррекции алкоритма. Откуда ещё одна несложная задача - как можно вообще избавиться от накапливания ошибки, при условии, что разрядность входного и выходного рядов заметно меньше, чем разрядность вычислителя.