Как найти период, доминирующий в серии timestamp'ов (с шумом)?

Question

[Сергей Соколов]

Набор данных – timestamp'ы неких событий. Некоторые из них приблизительно ложатся на периодическую сетку, а некоторые просто шум.

Как выделить примерно-периодические серии событий и отделить их от шума?

Например, несмотря на небольшие вариации, здесь период 5:
3, 8, 13.1, 17.9, 23.1, 28, 33

Добавим немного шума, который надо игнорировать:
3, 4, 8, 13.1, 17.9, 19, 23.1, 28, 33

Линейная регрессия помогла бы искать линейную ф-ю y = a + bx, но я хочу найти период. Как это сделать?

Upd. для Matlab нашёлся релевантный пример в описании ф-ии fft() – пример с зашумленным сигналом (Noisy Signal). Пытаюсь его применить к своему случаю. В матлабе полный нуб, поэтому со скрипом.

Answer 1

[Андрей]

Перевести в частотную область (преобразование Фурье) - доминирующий период даст Вам большое значение на соответствующей частоте.

Comments

[Сергей Соколов]

как регулировать чувствительность к отклонениям от точной сетки? В данных выше допускать отклонения в 0.1, а в 0.2 уже считать шумом, например.

[Сергей Соколов]

Преобразование Фурье от первого примера, без шума, даёт такой набор.Как из него получить "5" ?

[[126.1,0],[-17.522252093395643,36.328408396621214],[-17.590096886790246,13.842806382860214],[-17.437651019814087,4.267429532510812],[-17.43765101981412,-4.267429532510807],[-17.590096886790082,-13.842806382860239],[-17.522252093395473,-36.328408396621256]]

[Сергей]

Сергей Соколов: Вам на самом деле надо преобразовывать не последовательность.
А скажем делать преобразование фурье для единичных всплесков, которые происходят в моменты времени согласно вашей последовательности.
Тогда вы сможете построить АЧХ и увидеть главный период последовательности.

т.е. вы должны определить некую функцию f(t) которая принимает 1 в моменты времени обозначенные в вашей последовательности.

[Сергей Соколов]

Сергей: аа, т.е. надо добить нулями с шагом 0.1 !

[Сергей]

можно и так

[Сергей Соколов]

Сергей: похоже, FFT "не работает" на единичных пульсах.. Вот мои попытки. Может, надо их как-то блюрить до состояния, похожего на единичные синусоидальные wavelet'ы?

[Сергей]

Сергей Соколов: Можно и блюрить, просто гармоник тогда будет меньше.
Ваша гармоника самая первая в спектре в любом случае

Answer 2

[xmoonlight]

Задача решается достаточно просто:

SOURCE: Array ( [0] => 13.1 [1] => 3 [2] => 4 [3] => 8 [4] => 17.9 [5] => 19 [6] => 23.1 [7] => 28 [8] => 33 )
------------
SORTED: Array ( [0] => 3 [1] => 4 [2] => 8 [3] => 13.1 [4] => 17.9 [5] => 19 [6] => 23.1 [7] => 28 [8] => 33 )

DISTANCE: Array ( [0] => 1 [1] => 4 [2] => 5.1 [3] => 4.8 [4] => 1.1 [5] => 4.1 [6] => 4.9 [7] => 5 )
AVERAGE DISTANCE:3.75

CLUSTERS:

Array
(
    [-] => Array
        (
            [data] => Array
                (
                    [0] => 1
                    [1] => 1.1
                )

            [max] => 1.1
            [min] => 1
        )

    [+] => Array
        (
            [data] => Array
                (
                    [0] => 4
                    [1] => 5.1
                    [2] => 4.8
                    [3] => 4.1
                    [4] => 4.9
                    [5] => 5
                )

            [max] => 5.1
            [min] => 4
        )

)

AVERAGE PERIOD POSITIVE [+]: 4.65
AVERAGE PERIOD NEGATIVE [-]: 1.05
RANGE [+]: [4.55; 5.1]

3: 0, 1, 5, 10.1, 14.9, 16, 20.1, 25, 30, => 1
4: 1, 0, 4, 9.1, 13.9, 15, 19.1, 24, 29, => 0
8: 5, 4, 0, 5.1, 9.9, 11, 15.1, 20, 25, => 2
13.1: 10.1, 9.1, 5.1, 0, 4.8, 5.9, 10, 14.9, 19.9, => 2
17.9: 14.9, 13.9, 9.9, 4.8, 0, 1.1, 5.2, 10.1, 15.1, => 1
19: 16, 15, 11, 5.9, 1.1, 0, 4.1, 9, 14, => 0
23.1: 20.1, 19.1, 15.1, 10, 5.2, 4.1, 0, 4.9, 9.9, => 1
28: 25, 24, 20, 14.9, 10.1, 9, 4.9, 0, 5, => 2
33: 30, 29, 25, 19.9, 15.1, 14, 9.9, 5, 0, => 1
------------
CHAINS:

Array
(
    [0] => 1
    [1] => 0
    [2] => 2
    [3] => 2
    [4] => 1
    [5] => 0
    [6] => 1
    [7] => 2
    [8] => 1
)

------------
RESULT:

Array
(
    [0] => 3
    [1] => 8
    [2] => 13.1
    [3] => 17.9
    [4] => 23.1
    [5] => 28
    [6] => 33
)

DISTANCE RESULT: Array ( [0] => 5 [1] => 5.1 [2] => 4.8 [3] => 5.2 [4] => 4.9 [5] => 5 )

PERIOD (average): 5

Comments

[Сергей Соколов]

рассматривая только интервалы между соседями, а не через одного (и более), упускаем верные гипотезы: интервал 3–8 напр. Разве нет?

[Сергей Соколов]

xmoonlight: интуиция подсказывает, что может существовать подход, похожий на линейную регрессию, только вместо линейной функции угадывать период (и фазу). Находить ошибку – отклонение реальных данных от гипотезы – и стараться её минимизировать. Для пары точки, напр. 3 и 8, offset=3, period=5 будет оптимальной. Для них же offset=1 и period=2.5 тоже верная гипотеза, но более высокие гармоники должны иметь, видимо, вес меньший, чем низкие.

[xmoonlight]

Сергей Соколов: вот, я сделал то, что написал выше (в комменте): 5-ку - алгоритм получил автоматом! Работает!

[Сергей Соколов]

xmoonlight: но как? Мой хрустальный шар сломался и по дамрам php ход решения ге восстатавливантся полностью.

[Сергей Соколов]

Опишите, пожалуйста, алгоритм.