Боротьба з млинами
У цій статті ліричний герой кидає виклик оптимальної реалізації класичного поліноміального інтерполятора Лагранжа (Фарроу), в процесі битви випадково відкриває і доводить тривіальне нікому не потрібне математгічне заклинання, за допомогою якого намагається потіснити супротивника, але за результатами всіх раундів бою рішення.
— Де ви бачите велетнів? — спитав Санчо Панса. — Та ось вони, з величезними руками,— відповів пан. — У деяких із них довжина рук сягає майже двох миль. — Помилуйте, сеньйоре, — заперечив Санчо, — те, що там видніється, зовсім не велетні, а вітряки; те, що ви приймаєте за їхні руки, - це крила: вони кружляють від вітру і приводять у рух млинові жорна. — Зараз видно недосвідченого шукача пригод, — зауважив Дон Кіхот, — то велетні. І якщо ти боїшся, то від'їжджай убік і помолися, а я тим часом вступлю з ними в жорстокий і нерівний бій.
Одного разу мій друг попросив мене промоделювати йому простий і маловитратний варіант ресемплінгу (зміни частоти дискретизації) цифрового сигналу. Особливих вимог до лінійності АЧХ і ФЧХ не було, а бажання заощаджувати час і пам'ять було, тому частина методів типу поліфазних інтерполюючих банків FIR-фільтрів відразу була виключена з розгляду, а основну увагу приділено методам локальної поліноміальної інтерполяції. Поліноми швидко і легко рахувати, для зберігання їх коефіцієнтів потрібно мало пам'яті, і за допомогою одного набору коефіцієнтів можна вважати будь-які значення всередині вихідного інтервалу дискретизації.
Знайомимося із противником — реалізація Фарроу
Відомо, що через будь-які N точок (звідси і надалі для визначеності мибудемо розглядати рівномірну сітку - з однаковим кроком, хоча для багатьох описуваних технік і методів це не є обов'язковою умовою), так ось, через будь-які N точок з різними координатами абсцис можна провести єдиний поліном N-1 порядку, що однозначно визначається набором N коефіцієнтів (почитав , що LaTex на хабрі не безпосередньо, а через картинки, які важкі і недовговічні, тому обходжуся тільки верхніми-нижніми індексами в html-тегах):
Для зручності розрахунків нормуємо нашу рівномірну сітку аргументу до цілих чисел, на прикладі чотирьох точок виберемо зручний інтервал аргументу (Фарроу вибрав від [-2, 1] з центральним інтервалом від [-1, 0]) і запишемо систему лінійних рівнянь з умови проходження полінома через наші 4 точки:
Вирішувати можна як завгодно, хоча Гауссом матрицю звертати, але Фарроу оптимізував кількість операцій для розрахунку коефіцієнтів цього полінома:
Разом операцій:
- 1 множення на константу (1/6)
- 2 поділу на 2 (вважається швидко і в цілісній арифметиці та у форматі з плаваючою точкою)
- 8 додавань/віднімань
Перший претендент - сплайн Катмулла-Рома
Можна помітити, що як шуканий поліном, так і всі його похідні є лінійними комбінаціями невідомихкоефіцієнтів полінома Тому ми можемо задавати умови рівності похідних будь-якого порядку у вузлах сітки заданим значенням, і одержувати лінійну систему рівнянь коефіцієнти полінома. Локальний кубічний сплайн Ерміта будується за чотирма умовами: значення полінома та його перша похідна в краях інтервалу. А отримати наближені значення цих похідних у потрібних точках ми можемо, продиференціювавши інтерполюючий поліном Лагранжа, який проходить через них. Наприклад, нормуємо аргументи чотирьох точок до інтервалу [-1, 2], візьмемо три точки з аргументами, проведемо через них параболу, і розрахуємо її коефіцієнти. А оскільки нам потрібно лише значення першої похідної цієї параболи в центральній точці 0, то достатньо розрахувати її коефіцієнт при першому ступені аргументу, він і буде використовуватися як похідна для побудови сплайна. Похідна у правій межі інтервалу розраховується аналогічно. В результаті нехитрого розрахунку отримуємо таку систему рівнянь для коефіцієнтів сплайну:
рішення якої можна отримати у наступному вигляді:
Можемо помітити, що при послідовній (потоковій) обробці вхідних даних похідна лівої межі інтервалу та сама, що й похідна правого кордону попереднього інтервалу, і отже на кожному інтервалі нам необхідно розраховувати тільки одну похідну — у правій межі інтервалу, а значення в лівій межі інтервалу із запам'ятованої при попередньому розрахунку у правій. За допомогою такої економії на сірниках отримуємо таку кількість операцій для розрахунку коефіцієнтів полінома:
- 2 множення/поділу на 2
- 7 додавань/віднімань
Другий претендент - сплайн Ерміта по 6 точках
Якщо ми будемооцінювати похідні у краях центрального інтервалу за трьома, а, по п'яти точках, їх розрахунок трохи ускладниться, та інші формули не поменяются:
З урахуванням попереднього зауваження про необхідність розрахунку лише однієї похідної на кожному інтервалі кількість операцій:
- 1 множення на константу (2/3)
- 2 множення/поділу на ступінь 2
- 9 додавань/віднімань
Порівняння варіантів
За кількістю операцій сплайн Катмулла-Рома суттєво виграє у Фарроу, Ерміт по 6 точках трохи програє (на одне додавання). Також варто зазначити, що обидва сплайни мають перший порядок гладкості (безперервні нульову та першу похідні) внаслідок своєї побудови, на відміну від Фарроу, який не забезпечує безперервність похідної. Але залишається питання точності апроксимації функцій цими сплайнами. В якості тестової функції був обраний синус, проведена його апроксимація наведеними вище варіантами сплайнів для різного кроку сітки, який нормується у відносну величину кількості точок на період. Вимірювалася величина максимального відхилення від апроксимованої функції. Результати представлені на графіку, при виборі логарифмічного масштабу по обох осях залежності виходять майже лінійні, що сходяться в одну точку при зменшенні частоти дискретизації та наближенні її до частоти Найквіста (2 точки на період тестового синуса, що апроксимується).
На графіку видно, що сплайн Катмулла-Рома поступається Фарроу за точністю, а останній поступається Ерміту по 6 точках, при практично однаковій кількості необхідних розрахунку операцій. Здобути чисту перемогу не вдалося, але в принципі результат їм непоганий.
Заклинання
Ми можемо розвинути спільну ідею сплайну Ерміта на похідні будь-якихпорядків. Наприклад, побудуємо сплайн за такими умовами: його значення в краях інтервалу збігаються зі значеннями відліків у цих точках, ідругіпохідні в цих точках збігаються з їх оцінками, зробленими по параболах через трійки точок - тобто, все як у Катмулла-Рома тільки замість перших братимемо другі похідні — парабола, як поліном другого порядку, цілком дозволить нам їх розрахувати. І тут ми з подивом виявимо, що побудований таким чином сплайн повністю збігається з нашим шановним Фарроу, який є насправді поліном Лагранжа. Бажаючі можуть перевірити це, я навіть напишу формулу оцінки другою похідною параболою по 3 точках (хоча і ця тривіальна вправа):
Вдалий збіг? Випадково, бо мало точок і малий порядок полінома? Чисельно перевіривши цей факт на різній кількості точок та порядку поліномів, я переконався в тому, що заклинання працює:
Нехай ми маємо інтерполяційний поліном, побудований за деякими значеннями на парній кількості точок рівномірної сітки
Тоді всі парні похідні цього полінома в точці x0 збігаються з похідними тих самих порядків полінома інтерполяційного, побудованого по точках
а в точці x1 - полінома, побудованого за точками
Доведення. Розглянемо функцію на рівномірній сітці з непарної кількості точок - зробимо зміщення аргументу, що наводить аргумент центрального вузла в 0
Існує її єдиний інтерполяційний поліном P0(x). Додамо до цього набору точок ще одну xk+1 (для доказу не важливо, додаємо ми її праворуч, ліворуч або взагалі в середину і не до вузлів вихідної сітки). Інтерполяційний поліном на розширеному наборі крапок за умови рівномірної сітки можна подати у формі Ньютона
де A – якась константа.Дійсно, від цієї адитивної добавки потрібно нульові значення в точках вихідного полінома, а константа A визначається з умови рівності P1(xk+1) = yk+1 – необхідне значення полінома у доданій точці. Через рівномірність сітки -x-i = xi, значить
Очевидно, що у правій частині стоїть сума P0(x) та полінома з ненульовими коефіцієнтами тільки при непарних ступенях аргументу. Отже, коефіцієнти при парних ступенях у поліномів P1(x) і P0(x) збігаються, проте парні похідні будь-якого полінома в нулі визначаються лише відповідним його коефіцієнтом при парному ступені аргументу. Таким чином, всі парні похідні P1(x) та P0(x) збігаються. Наведені вище міркування справедливі для будь-якого розташування доданої точки xk+1, зокрема при додаванні її праворуч або ліворуч по вихідній сітці. Таким чином, похідні парних порядків у нулі у вихідного полінома і полінома розширеного набору точок збігаються. Ч.т.д.
Останній раунд - Лагранж по Фарроу з Лагранжем через похідні.
Оскільки інтерполяційний поліном мінімального порядку є єдиним, то ми можемо з урахуванням щойно доведеного заклинання будувати його через умови на значення парних похідних у краях центрального інтервалу — наприклад, поліном 5-го порядку через 6 точок за 6 умовами: значення відліків і других та четвертих похідних , розрахованих за кожними п'ятьма точками. При цьому користуватися тією ж оптимізацією, що вище описана, — розрахунок похідних тільки в правій межі кожного інтервалу. У початковій постановці проходження полінома через усі точки ці інваріанти були очевидні. Застосуємо це до нашого вихідного Лагранжу по 4 точках та оптимізуємо кількість операцій. Наведу одразу код Matlab свого варіанту, та варіанта учасника з ніком _Anatoliy з форумуелектронікса, який незалежно та самостійно перемагав Фарроу ще до моїх експериментів:
Як видно, обидва коди містять меншу кількість операцій для розрахунку коефіцієнтів полінома, тож у плані спорту я вважаю, що ми обидва здобули перемоги. І хоча висловлюються думки, що це економія на сірниках, що в наш час, коли космічні кораблі борознять скрізь є математичні співпроцесори з інтегрованими операціями з числами з плаваючою точкою, швидкість операцій можна порівняти зі швидкістю переміщення вмісту регістрів і тим більше читання/запису пам'яті для збереження /відновлення розрахованих значень - подібні оптимізації нічого не дають, у мене все одно є якесь почуття задоволеності результатами, незважаючи на їх слабку практичну значимість.
А у нас тут можна отримати грант на тестовий період Яндекс.Хмари. Варто лише у полі «секретний пароль» запровадити «Хабр»