Choliy.NumericaMethods.2011Dec01

10. Часовi ряди

У цьому роздiлi розглядаються деякi окремi питання аналiзу часових рядiв. Часовий ряд це сукупнiсть вимiрiв деякої змiнної h(t) зi зростанням часу. Такi ряди звичайнi для фiзичного експерименту, як правило, скiнченого в часi. Однак, спостережну величину, що є функцiєю часу, неможливо спостерiгати неперервно. Тому її дискретизують, тобто знiмають (фiксують, записують) у моменти t1, t2 . . . , якi надалi вважаються рiвновiддаленими. Якщо iнтервал мiж вiдлiками позначити через , то отриманi в результатi дискретизацiї данi це ряд:

hn = h(nΔ), n = . . . , −2, −1, 0, 1, 2, . . . ,

або, якщо так зручнiше, n = 0, 1, . . . . Величину називають iнтервалом вибiрки. Символ i у цьому роздiлi завжди позначає комплексну одиницю, а f - частоту.

Перетворення Фур’є. Фiзичний процес можна описувати не тiльки у часовiй областi, але й у частотнiй. Теоретичну базу для цього дає перетворення Фур’є:

Загальноприйнято Фур’є перетворення функцiї позначати такою ж буквою, але у верхньому регiстрi.

Маючи iнтервал вибiрки можна визначити частоту Найквiста:

1 fN = 2Δ,

це найбiльша частота, яку можна розрiзнити в даних при такому iнтервалi вибiрки. Дуже зручно вважати, що = 1, тодi частота Найквiста = 1/2. Якщо функцiя h(t) не мiстить частот, вищих, нiж fN , то при заданому жодна частота втрачена не буде, а h(t) буде повнiстю визначатися своїми hn. Однак, якщо fN недостатньо висока, то вищi частоти f > fN не вдасться

91

зберегти, навiть гiрше, вони появляться в областi [0, fN ] як наведенi частоти, яких насправдi у оригiнальному рядi немає. Ця проблема називається накладанням частот (aliasing) i за смислом означає маскування високих частот низькими. На рис. 10.1, тоненькою лiнiєю показано графiк функцiї, яка насправдi послужила для утворення часового ряду, що вибирався у моменти, якi вiдповiдають квадратним мiткам. Очевидно, що цi квадратики лягають на зовсiм iнше коливання, нiж те, що показане тоненькими лiнiями. Таким чином висока частота замаскувалася пiд нижчу частоту.

Рис. 10.1: Накладання частот.

Надалi вважаємо, що ми маємо N експериментальних точок hk = h(tk), k = 0, . . . , N − 1, tk = k , N - парне число. Розглянемо перетворення Фур’є на дискретному наборi частот

92

ми поклали = 1, тому в кiнцевiй формулi його немає. Отже, дискретнi версiї перетворень Фур’є:

У англомовнiй лiтературi часто використовується скорочення DFT = Discrete Fourier Transform, дискретне перетворення Фур’є.

Швидке перетворення Фур’є. Це не нове математичне перетворення, це лише алгоритм, який суттєво пришвидшує обчислення. В англомовнiй лiтературi вiн вiдомий як FFT (Fast Fourier Transform). Вiзьмемо вираз для прямого перетворення i суму, що стоїть в ньому розiб’ємо на суму парних та непарних доданкiв:

Перший доданок є N/2 точковим перетворенням Фур’є лише парних точок ряду, другий - теж N/2 точковим перетворенням, але лише непарних точок. Множник, що виявився за сумою, визначається частотою, для якої знаходиться повне перетворення, та початковою кiлькiстю точок. Важливим у цьому алгоритмi є той факт, що кiлькiсть точок парна. Ще краще, якщо кiлькiсть точок вибрати рiвною степенi двiйки N = 2p, тодi вказану формулу можна рекурсивно використовувати ”вглиб“ аж p − 1 разiв. На найглибшому рiвнi залишається лише двi точки, а перетворення стає тривiальним, оскiльки при N = 2, k може приймати лише два значення 0 та 1:

93

Алгоритм FFT для N = 2p потребує лише N log2 N операцiй, що робить його одним з найшвидших комп’ютерних алгоритмiв. Його використання надзвичайно широке. Алгоритм легко паралелiзується, що дає додатковi переваги. Вiдомi також алгоритми FFT для довжини послiдовностi, що виражається добутком степенiв iнших простих чисел, але вони знайшли застосування лише у спецiальних задачах. Такi алгоритми називаються алгоритмами Винограда.

Перiодограма. Якщо представити собi окремi вiдлiки часового ряду як вiдхилення вiд умовного нуля, то сума квадратiв усiх вiдлiкiв пропорцiйна енергiї системи, що згенерувала ряд. Перетворення Фур’є переводить енергiю у частотну область i розподiляє її по частотах. Важливе значення має теорема Парсеваля, згiдно з якою

що є ще одним записом для закону збереження енергiї. Амплiтуда кожного з доданкiв Hk показує як багато енергiї зосередилось на частотi kN , а

залежнiсть Hk2 вiд k показує розподiл енергiї по частотах. Таку залежнiсть називають перiодограмою.

Щоб отримати перiодограму, знаходимо Фур’є перетворення ряду, вiдлiки перiодограми будуть такими:

P (0) = H02/N2,

P (fk) = (|Hk|2 + |HN−k|2)/N2, k = 1, . . . , N/2 − 1

P (fN ) = |HN/2|2/N2.

Бiльш глибокi питання, як-от статистичнi властивостi перiодограми, оцiнка правомiрностi та якiсть отриманих оцiнок виходять за межi нашого курсу.

Кореляцiйнi функцiї. Коварiацiйною функцiєю двох рядiв, g(tk) = gk та h(tk) = hk з однаковим iнтервалом вибiрки називається:

ˆ

де g,ˆ h - вибiрковi середнi значення процесiв. Коварiацiйна функцiя ряду з собою називається автоковарiацiйною функцiєю:

94

Як видно з означення, коварiацiя залежить вiд j, що по смислу є зсувом одного ряду вiдносно другого. Автоковарiацiя при j = 0 дає варiацiю ряду. Очевидно, що j [−N + 1, N − 1]

Аналогiчно вводяться кореляцiйнi функцiї, якi є коварiацiями, нормованими на варiацiю: автокореляцiйна функцiя:

c(j) = rr(0)(j),

чи взаємна кореляцiйна функцiя рядiв gk, hk:

c(j) = r(j) . σg · σh

Кореляцiйнi функцiї використовують тодi, коли потрiбно вияснити наскiльки ряди ”схожi“ мiж собою, чи може, зсунута версiя ряду ”схожа“ на сам ряд. Значення кореляцiйної функцiї за модулем не перевищують одиницi, одиниця вказує на повну тотожнiсть рядiв, мiнус одиниця - на їх цiлковиту несхожiсть. Ряд, отриманий дискретизацiєю синуса буде схожий на ряд, отриманий з косинуса, зсунутий на чверть перiода вперед i зовсiм не схожий на нього ж, зсунутого на чверть перiода назад. Кореляцiйна функцiя повнiстю випадкового ряду рiвна нулю всюди, окрiм нульового значення. Такий ряд не корелює нi з яким iншим рядом. Перетворення Фур’є вiд автокореляцiйної функцiї називають спектром.

Оскiльки фiзичнi ряди рiдко коли є цiлком випадковими, то кожна точка ряду hk якось пов’язана (корелює) з попереднiми, наприклад, з p попереднiми, що можна описати таким лiнiйним спiввiдношнням:

hk = a1hk−1 + a2hk−2 + a3hk−3 + · · · + aphk−p + ek,

де ai - коефiцiєнти залежностi (регресiї), ek - додаткова випадкова похибка в момент tk. Домноживши рiвняння на hk−1 i взявши суму по всiх точках, отримуємо:

95

Проробивши цю процедуру p разiв, кожен раз з iншим множником hk−2, hk−3, . . . , складаємо систему p рiвнянь:

з якої знаходимо ai, i = 1, . . . , p. При вiдомих значеннях регресiйних коефiцiєнтiв ai авторегресiйний спектр послiдовностi hk для будь-якої частоти f [0, fN ] визначається виразом:

Особливiстю цiєї оцiнки спектра є те, що в областi максимумiв спектр зростає експоненцiйно, що робить його пiки гострими, на вiдмiну вiд перетворення Фур’є, пiки якого квадратичнi по частотi i виглядають значно пласкiшими. Роздiльна здатнiсть по частотi у такого спектру висока, однак, амплiтуди залежать вiд енергiї нелiнiйно. Зазвичай, спектри використовують для знаходження частот, чи гармонiк, що присутнi в рядi. Амплiтуду та фазу гармонiк визначають методом найменших квадратiв.

96

Задачi. Всi задачi цього роздiлу виконуються на штучних часових рядах, якi потрiбно попередньо самостiйно згенерувати. Зазвичай, це сума одногодвох синусiв та випадкової розподiленої нормально компоненти, з заданою дисперсiєю:

X

hj = Ak sin (ωktj + ϕk) + εj.

k

Результат виконання задачi - таблиця для всiх допустимих значеннях аргумента.

1.Згенеруйте 1024 точки часового ряду: 5 sin (πt/18). Знайдiть його DFT.

2.Знайти перетворення Фур’є ряду з задачi 1, скориставшись FFT.

3.Побудувати перiодограму ряду з задачi 1.

4.Знайти автокореляцiйну функцiю ряду з задачi 1, при j = 0, . . . , 30.

5.Побудувати спектр ряду з задачi 1.

6.Побудувати авторегресiйний спектр ряду з задачi 1.

7.Згенеруйте 1024 точки часового ряду: 5 sin (πt/18) + 1.5ε, де ε - випадкове рiвномiрно розподiлене на [0, 1] число. Знайдiть його DFT.

8.Знайти перетворення Фур’є ряду з задачi 7, скориставшись FFT.

9.Побудувати перiодограму ряду з задачi 7.

10.Знайти автокореляцiйну функцiю ряду з задачi 7, при j = 0, . . . , 30.

11.Побудувати спектр ряду з задачi 7.

12.Побудувати авторегресiйний спектр ряду з задачi 7.

13.Згенеруйте 1024 точки часового ряду: 10 sin (πt/18 + π) + 5 sin (πt/80). Знайдiть його DFT.

14.Знайти перетворення Фур’є ряду з задачi 13, скориставшись FFT.

15.Побудувати перiодограму ряду з задачi 13.

16.Знайти автокореляцiйну функцiю ряду з задачi 13, при j = 0, . . . , 30.

17.Побудувати спектр ряду з задачi 13.

18.Побудувати авторегресiйний спектр ряду з задачi 13.

19.Згенеруйте 1024 точки часового ряду: 5 sin (πt/18 + π) + 10 sin (πt/80) + 0.1t+ε, де ε - випадкове рiвномiрно розподiлене на [0, 1] число. Знайдiть його DFT.

20.Знайти перетворення Фур’є ряду з задачi 19, скориставшись FFT.

21.Побудувати перiодограму ряду з задачi 19.

22.Знайти автокореляцiйну функцiю ряду з задачi 19, при j = 0, . . . , 30.

23.Побудувати спектр ряду з задачi 19.

24.Побудувати авторегресiйний спектр ряду з задачi 19.

97

25.Згенеруйте 1024 точки часового ряду: 5 sin (πt/18 + π) + sin (πt/150) + ε, де ε - випадкове рiвномiрно розподiлене на [0, 1] число. Знайдiть його DFT.

26.Знайти перетворення Фур’є ряду з задачi 19, скориставшись FFT.

27.Побудувати перiодограму ряду з задачi 19.

28.Знайти автокореляцiйну функцiю ряду з задачi 19, при j = 0, . . . , 30.

29.Побудувати спектр ряду з задачi 19.

30.Побудувати авторегресiйний спектр ряду з задачi 19.

31.Знайти FFT ряду, що складається з 1024 випадкових, рiвномiрно розподiлених на [0, 1] точок.

32.Знайти автокореляцiйну функцiю ряду, що складається з 1024 випадкових, рiвномiрно розподiлених на [0, 1] точок.

98

11. Диференцiювання

99

Якщо потрiбно знайти похiднi в точках таблицi, для яких q є цiлим числом,

Диференцiювання за Лагранжем. Будемо дiяти аналогiчно, але скористаємося тепер полiномом Лагранжа:

де Π(x) = (x − x0)(x − x1) . . . (x − xn).

Трохи перепишемо цей вираз. Оскiльки вузли рiвновiддаленi, hq = x − x0, hi = xi − x0, поперереднi спiввiдношення приймають виду:

Π(x) = hn+1q(q − 1)(q − 2) . . . (q − n),

Π0(xi) = (xi−x0)(xi−x1) . . . (xi−xi−1)(xi−xi+1) . . . (xi−xn) = hni!(n−i)!(−1)n−i.

Вираз для полiному Лагранжа, записаний з використанням q:

Якщо задатися якимось значенням n i взяти похiднi пiд знаком суми, отримаємо формули, придатнi для обчислень.

100