CC BY
98
Труды Карельского научного центра РАН № 4. 2014. С. 106-115
УДК 519.652
ПРИМЕНЕНИЕ КВАЗИЭРМИТОВЫХ КУБИЧЕСКИХ СПЛАЙНОВ ДЛЯ ПЕРЕДИСКРЕТИЗАЦИИ ЗВУКОВЫХ ФАЙЛОВ
С. Е. Михеев, П. Д. Морозов
Санкт-Петербургский государственный университет
Проведены исследования передискретизации с увеличением частоты сэмплирования цифровых звуковых файлов в К раз. Для сохранения частотных характеристик звука переход на большую частоту воспроизведения сопровождается вставкой между каждыми двумя исходными соседними сэмплами К — 1 добавочного сэмпла. Последние получались интерполяцией сплайнами двух разных типов — линейными и квазиэрмитовыми кубическими (КЭК-сплайнами). На концах звена КЭК-сплайна взамен обычных, применяемых для построения эрмитовых сплайнов, вычислялись центральные разностные производные по сэмплам исходного файла.Если под уровнем гладкости понимать максимум величины, обратной модулю второй производной в точках ее непрерывности, то доказано, что передискретизация с увеличением частоты дискретизации в К раз с помощью линейной интерполяции снижает в К раз в окрестностях исходных узлов уровень гладкости выходного сигнала, но при передискретизации с помощью КЭК-сплайнов уровень гладкости восстанавливаемого сигнала снизится не более чем в 3 раза для любых К.
Ключевые слова: интерполяция, сплайн, сэмпл, разностная производная, выпуклость, гладкость.
S. Е. Mikheev, P. D. Morozov. APPLICATION OF QUASI-HERMITIAN CUBIC SPLINES FOR OVERSAMPLING OF AUDIO FILES
Oversampling with a АГ-fold increase in the rate of sampling digital audio files was studied. To save the sound frequency characteristics, the transition to a higher frequency should be accompanied by the insertion of К — 1 additional samples between each two adjacent original samples. The additional samples were generated by interpolation with two different splines — quasihermitian cubic (QHC-spline) and linear ones. At the end of the QHC-spline link central difference derivatives were calculated from the source file samples instead of the genuine ones used to construct hermitian splines. If the smoothness level is defined as the maximum of inverted second derivative modulo in the points of its continuity, it is proved that the oversampling with the sampling frequency increasing up by К times and with linear interpolation reduces the output signal smoothness level АГ-fold in the neighborhood of the original node, but the oversampling by QHC-spline would drop the smoothness level of the recovered signal no more than 3 times for all K.
Key words: interpolation, spline, sample, difference derivative, convexity, smoothness.
Введение
Неизбежные потери при передаче и хранении аналогового звука и, что, наверное, самое существенное, при копировании аудиозаписей открыли путь для все большего распространения цифрового звука, лишенного этих недостатков. Однако в цифровом звуке есть свои недостатки, не характерные для аналогового, а именно утеря информации при дискретизации и появление паразитных шумов при восстановлении аналогового сигнала по дискретным значениям. Переводом аналогового звука в цифровой, или оцифровкой, занимаются аналого-цифровые преобразователи (АЦП). Какова бы ни была логика работы ныне существующих АЦП, каждое число Sj в выходном потоке {sj}o, именуемое сэмплом, интерпретируется как некое значение аналогового сигнала S, относимое К интервалу времени далее сэмплъный интервал. Таков, например, формат WAV. (Сигнал s — функция времени t.)
Длина сэмпльных интервалов времени назначается одинаковой в течение одного сеанса звукозаписи. Число сэмпльных интервалов в одной секунде называется частотой дискретизации. Это важнейшая характеристика качества работы АЦП. Чем больше частота — тем больше адекватность цифрового сигнала аналоговому. Однако увеличение частоты приводит к повышению: интенсивности потока цифровых данных при ограниченной пропускной возможности аппаратуры, вычислительной нагрузки на цифровые процессоры, объема необходимой для хранения цифрового сигнала памяти. В зависимости от конкретной ситуации выбирается лучшее соотношение между уровнем адекватности и объемом хранения. Для этого в АЦП аппаратно реализуется несколько стандартных частот дискретизации: 11025, 22050, 32000, 44100, 48000, 88200, 96000 Гц.
Восстановлением аналогового звука из цифрового занимаются цифро-аналоговые преобразователи (ЦАП), работающие на тех же стандартных частотах дискретизации. Работу ЦАП можно описать примерно следующим образом: поток сэмплов {sj}o преобразуется в аналоговый ступенчатый сигнал, который подается на аналоговый фильтр, выполняющий аппаратное сглаживание. Что влияет на уровень паразитных шумов? ЦАПом с идеальной выходной характеристикой считается такой, который выдает ступенчатый сигнал-напряжение прямоугольной формы, величина
которого строго равна сэмплу, а длительность
— сэмпльному интервалу дискретизации, при этом величина напряжения строго пропорциональна сэмплу. В действительности сигнал имеет не совсем прямоугольную форму, связанную с тем, что существует некоторое время установления, которое «размазывает» строго вертикальное начало ступеньки и добавляет еще некоторые затухающие колебания. При этом, к сожалению, данные отклонения от строгой ступенчатости нелинейно связаны с перепадом напряжения от ступеньки к ступеньке. Кроме того, величине этих перепадов соответствуют амплитуды гармоник, в разложении ступенчатого сигнала, с частотами выше частот дискретизации. Конечно, при частоте дискретизации более 20 кГц эти гармоники лежат вне диапазона восприятия человеческого уха, и на данную разновидность шумов можно не обращать внимание при конструировании аналогового фильтра. Но в низкочастотном диапазоне — диапазоне слышимости человеческого уха — от них могут возникнуть обертоны, которые, частично проникая сквозь аналоговый фильтр, отрицательно сказываются на качестве воспроизведения выходного сигнала.
Простейшим способом уменьшить время установления и облегчить работу аналогового фильтра является передискретизация с кратным увеличением частоты с помощью интерполяции. Она позволяет уменьшить перепады напряжения между ступеньками напряжения на выходе из ЦАПа. А уменьшение таких амплитуд облегчает дальнейшую работу включенного на выходе из ЦАПа аналогового фильтра, сглаживающего ступенчатый сигнал. Благоприятные последствия: снижение стоимости аналогового фильтра и повышение качества окончательного выходного сигнала. Кратная передискретизация в К раз заключается в том, что в поток сэмплов цифрового звука с частотой дискретизации Р между каждой парой соседних исходных сэмплов вставляется К — 1 новый сэмпл, после чего пополненный поток сэмплов воспроизводится как цифровой звук с частотой дискретизации КР. Вставку можно делать различными способами. Если, например, использовать линейную интерполяцию, то перепады напряжения между ступеньками в ступенчатом сигнале уменьшаются в К раз, с соответствующим снижением паразитных шумов на выходе из ЦАПа и аналогового фильтра.
Интерполирование разными типами сплайнов для передискретизации (а значит, и для
повышения качества воспроизведения) и составило исследовательский интерес. В частности, сплайнами вида (порядка 3, дефекта 2), сглаженными в соответствии с алгоритмом в [3] до сплайнов вида й'зд (дефекта 1), а также самыми простыми линейными сплайнами 5^1 (порядка 1, дефекта 1).
Описание задачи
Оцифрованный сигнал представим в виде последовательности СЭМПЛОВ — чисел {«г}о, привязанных к сэмпльным интервалам между равномерно расположенными с шагом дискретизации к моментами времени {4г}о+1 и соответствующих усредненным значениям ИСХОДНОГО аналогового сигнала (звукового, электромагнитного) на этих интервалах либо на более узких подынтервалах:
[и, и + е)с [и, и+1], г = О, ЛГ, £<к.
В том же случае, когда е -С /г, говорят, что од? - - - > являются отсечками входного сигнала 3 В моменты соответственно Тре-
буется найти интерполяционный алгоритм, который приближал бы значения входного сигнала в промежуточных моментах времени (дополнительных узлах), в совокупности с исходными узлами образующих тоже равномерную сетку. Требования к алгоритму:
1. В исходных узлах 4о, Ь,... возвращает значения входного сигнала.
2. Для получения дополнительных сэмплов Зц на дополнительных узлах из сэмпльно-го интервала [и,и+1] используются только значения входного сигнала в узлах, непосредственно примыкающих к Ь\ и ^+1. Количество таких узлов желательно минимизировать.
3. Аналоговый сигнал после прохождения выходного цифрового сигнала через ЦАП должен иметь наименьшее количество паразитных шумов.
Изменение частот дискретизации
Рассмотрим учетверение частоты. При сохранении частотных характеристик входного сигнала на выходе с каждым входным сэмплом сопоставляются 4 сэмпла выходного сигнала, то есть требуется назначить дополнительно 3 сэмпла. В связи с чем и производилась интерполяция на каждом сэмпльном интервале. Были рассмотрены линейная интерполяция и кубическая квазиэрмитовая. Оба варианта порождают локальные сплайны, то есть не требующие одновременной обработки всего множества узлов.
Линейная интерполяция позволяла обрабатывать неограниченную последовательность входных сэмплов. Результирующая сту-
пенчатая функция S(t) := Sij, t G [U,ti+1), при подаче в ЦАП достаточно хорошо сглаживалась аппаратно на дополнительных узлах и несколько хуже на исходных, то есть на каждом четвертом. Иными словами, если построить линейный сплайн по новой совокупности сэмплов, то он имеет разрывы производной на старых узлах, и им соответствует снижение качества аппаратного сглаживания.
В связи с этим возникает предположение, что если, сравнительно с линейным сплайном, пожертвовать немного гладкостью на дополнительных узлах, но уменьшить разрывы производной на старых узлах, то можно ожидать общего улучшения качества воспроизведения. Такое перераспределение скачков производной с исходных узлов на дополнительные возможно при получении дополнительных сэмплов не линейной интерполяции, а такой, которая обеспечивает непрерывность производной сплайна, получаемого «склейкой» звеньев — интерполяционных многочленов между исходными узлами. Среди различных вариантов таких сплайнов выберем наиболее простой.
Квазиэрмитовы кубические сплайны
Классический эрмитов сплайн состоит из звеньев, являющихся интерполяционными полиномами Эрмита, строящимися по значениям интерполируемой функции и ее производной на двух узлах, ограничивающих участок интерполяции данного звена [1]. «Склейка» таких звеньев дает эрмитов сплайн. Его тип — S3,2, то есть он гладкий.
В нашем случае входной сигнал представляет собой отсечки значений звукового давления в моменты {U}o+1, но информация о значениях производных давления в какие-либо моменты отсутствует. И имеющиеся теория и алгоритмы позволяют строить кубический сплайн (типа б'зд), то есть еще глаже эрмито-вого. Однако такой сплайн уже не относится к категории локальных, и в вычислении четырех параметров каждого звена используется информация по всем N + 2 узлам. При большом N (сотни тысяч) задача становится практически неразрешимой. Путь для выхода из коллизии между вычислительной трудоемкостью и гладкостью сплайна вполне естествен: получить приближенные значения производных в узлах интерполяции по некоторому ограниченному шаблону, а потом построить полином Эрмита. И уже в 1880-м году Т. В. Sprague [8] фактически предложил формулы для звеньев сплайна 3, которые строились как полиномы Эрмита по значениям на концах сэмпльного интервала звена 1-й и 2-й про-
0
изводных интерполируемой функции, вычисляемых приближенно с помощью полинома Лагранжа по шеститочечному шаблону, и значениям интерполируемой функции. В 1899 году Каруп (J. Karup) [5] предложил формулы для звеньев сплайна S^2 с вычислением производных на концах сэмпльного интервала звена по четырехточечному шаблону. С введением понятия сплайна в 1946 году [7] и соответствующей терминологии сплайны, конструируемые согласно формулам Sprague-Karup, естественно называть квазиэрмитовыми кубическими сплайнами (КЭК-сплайнами). Однако и до конца двадцатого века в некоторых публикациях можно обнаружить название «касательная интерполяция» (osculatory interpolation) для метода Карупа. В частности, Акима в [4, см. (10)] под этим названием привел небольшую модификацию вычислительных формул Карупа, в конечном итоге приводящую к тому же сплайну, что и метод Карупа. Мы здесь построим КЭК-сплайны, исходя из других начальных соображений, и получим расчетные формулы, несколько отличные от формул Карупа и Акимы. И хотя итоговый сплайн оказывается таким же, предлагаемая методика его построения позволила провести анализ уровня гладкости (см. опр. 1) выходного сигнала, который появляется после сглаживания в ЦАПе.
Существующая методика передискретизации довольно-таки далека от идеи интерполяции. Так, согласно [6, гл. 10.2], все дополнительные сэмплы Sij на сэмпльном интервале назначаются равными 0 с дальнейшим сглаживанием цифровым фильтром нижних частот, то есть даже без линейной интерполяции, с применением «разрывных функций». Поскольку предлагаемая методика передискретизации нацелена на улучшение существующей методики, здесь будут исследованы только линейные сплайны Бц, обеспечивающие «непрерывность», и квазиэрмитовы кубические S3 2, обеспечивающие «гладкость».
Построение КЭК-сплайна.
С каждым сэмплом, начиная с третьего и до предпоследнего, сопоставим центральную разностную производную:
si+l — Sj_i .
ггц =-----—-----, г = 3,4,...
Со вторым сэмплом сопоставим левую разностную производную:
Начальное звено сплайна получим линейной интерполяцией по первому и второму сэмплам. Впрочем, можно было бы интерполировать между первым узлом и вторым как угодно — последствия этого произвола, как будет видно в дальнейшем, закончатся к третьему узлу. Остальные звенья от этого не зависят.
Прочие звенья между моментами и и ^+1, г ^ 1, получаются в результате решения задачи кратного интерполирования с помощью полинома третьей степени Нз(Ь):
— sji — rrij,
(1)
где ^ = і, і + 1. Поскольку т-) не являются значениями производной интерполируемой функции в, такой полином назовем квазиэр-митовым кубическим, а КЭК-сплайны и есть сплайны, составленные из таких полиномов.
Свойства КЭК-сплайнов
Решение задачи (1) с вычислительной точки зрения удобнее описать не через глобальную переменную і, а через локальную г следующим образом:
г =
t — entier(i//i)/i h
(2)
для всех Ь е [епЫег(Ь/к)Ь, епЫег(1/к)Ь + к\. Здесь еп^ег — функция, выдающая наибольшее целое, не превосходящее аргумента. Из (2) видно, что г пробегает значения от 0 до 1, когда Ь пробегает значения между двумя соседними узлами. Формула (2), очевидно, не позволяет определить 4 по г, однако при построении интерполяционного сплайна этого не потребуется.
Решение задачи (1) в локальных переменных доставляется известной формулой:
С[г,1+1](т) = (l — т)2(l + 2т)si + т2(3 — 2т)si+l +
+ т( 1 — т)2кгПг — т2( 1 — т)кгПг+(3)
Не умаляя общности, можно положить единицу времени равной шагу дискретизации (Ь = 1) и Ьг = г.
Соберем коэффициенты при одинаковых степенях гв (3):
С[г, 1+1] (г) = а*Г3 + ЬгТ2 +СгТ +
где
т2 =
S2 - Si
Щ — 2(^+2 + 3sj Sj_i),
h = ^(“si+2 + 4si+i - 5Si + 2sj_i),
s%-\-1 s%—i ,
5 di —
Иллюстрация поведения г-го звена сплайна тесно связана с параллелограммом звена, который образуется касательной к звену в узле г и параллельной ему прямой, содержащей отрезок совместно с касательной к
звену в узле г + 1 и параллельной ему прямой, содержащей отрезок Х^Х^, где Если упомянутые отрезки не параллельны, то параллелограмм звена, очевидно, существует. Поскольку у параллелограмма звена нет вертикальных сторон, он всегда будет иметь две боковые вершины (левую и правую). Оставшиеся две вершины можно именовать как верхнюю и нижнюю.
Рис. 1. Боковая вершина (правая) параллелограмма звена образована касательными
На рисунке 1 черными жирными линиями обозначены содержащие отрезки Х^\Х^\ и ХгХг+2 прямые, тонкими серыми — касательные к сплайну в узлах г и г + 1, которые параллельны соответствующим отрезкам.
Сравнение разных способов получения дополнительных сэмплов
Выясним, какие изменения можно ожидать в уровне паразитных шумов при разных типах интерполяции.
В стандартных способах кратной передискретизации обычно применяется цифровой фильтр нижних частот (ФНЧ), предшествующий ЦАПу. Так, в работе [6] при повышении частоты дискретизации в К раз предлагается вставка к := К — 1 дополнительных нулевых сэмплов между двумя соседними сэмплами входного сигнала, после чего, до поступления в ЦАП, расширенный поток сэмплов пропускается через ФНЧ.
Можно пойти еще дальше. Если уж есть возможность производить предварительную логическую обработку цифрового сигнала до ЦАПа, то следует изменить входной поток сэмплов так, чтобы аналоговый сигнал на выходе из ЦАПа оказался наиболее удобным для фильтрации нижних частот выходным аналоговым фильтром.
Рассмотрим разные типы интерполяции с точки зрения порождаемых ими шумов. Что происходит на стыках звеньев интерполяции?
Пусть {sj}o — сэмплы входного сигнала, sn,..., Site — дополнительные сэмплы в моменты (в дополнительных узлах), соответственно,
tij=i + Sj, j = (4)
где S — шаг новой дискретизации, S =
h/K !„_!= 1/К.
Помимо левых разностных производных полезны будут правые разностные произ-
водные
sf(h) :=
®г+1
h ’
< N
и центральные разностные производные вто-
+-,и^ 4(Н)~4-ЛН)
рого порядка [п) := ----------—-------, где
к — шаг исходной дискретизации.
Конечную разность первого порядка определим классическим образом: Д(г) = Яг+1 — (то есть правая конечная разность), а конечную разность второго порядка будем использовать левую:
А2(г) = Д(г) - Д(г - 1) = в*+1 - 2si + Яг-ь (5)
Линейная интерполяция
При кратном увеличении частоты дискретизации с помощью линейной интерполяции разностные производные остаются неизменными. Действительно, в^к/К) =
К--------, а числитель последней дроби в К
п
раз меньше величины — Si, г < N.
Но центральные разностные производные второго порядка увеличиваются в К раз. Для г ^ 1 было и стало соответственно
8+-(h) :=
зГ(6) :=
4(V -
= К
28
4(h) - 4- i(h) 2 h '
Числители двух последних дробей при линейной интерполяции равны. Следовательно, 4~(Ь/К) = Кз^~(к). Отметим, что в новой дискретизации разностные производные второго порядка в дополнительных узлах обнуляются: = 0, ] = ъ =
1,..., N — 1.Такая ситуация со вторыми разностными производными соответствует появлению шумов в сигнале с периодами, близкими к 46 [2], то есть четырехкратно пре-
восходящими периоды шумов от нового шага
дискретизации. Отсеивание низкочастотных шумов предъявляет дополнительные требования к ФНЧ.
Все конечные разности первого порядка на одном звене уменьшаются в К раз относительно исходной конечной разности, соответствующей этому звену, поэтому новые конечные разности второго порядка в исходном узле
д;ю = л2»/*=дл2«<« (6)
тоже уменьшаются в К раз, но в дополнительных узлах они равны нулю.
Не вдаваясь в детали работы ФНЧ по входному ступенчатому сигналу, можно считать, что о ее высоком качестве будет свидетельствовать близость отфильтрованного сигнала к некоторому аналоговому сигналу, который соответствует гладкой кривой, с кусочнонепрерывной второй производной и проходящей через все начала ступенек входного сигнала. Передискретизация может значительно понизить уровень гладкости этого ассоциированного аналогового сигнала.
Определение 1. Под уровнем гладкости будем понимать максимум величины, обратной модулю второй производной в точках ее непрерывности.
Ступенчатому сигналу соответствует бесконечное множество ассоциированных аналоговых сигналов различного уровня гладкости. На нем можно выделить подмножество с наибольшим уровнем гладкости. Множество ассоциированных сигналов передискретизирован-ного ступенчатого сигнала, разумеется, является подмножеством ассоциированных сигналов входного ступенчатого сигнала. Поэтому его уровень гладкости может только понизиться. Приведенные выше выкладки устанавливают предел такого снижения.
Лемма 1. Пусть непрерывная вещественная функция д имеет кусочно-непрерывную производную на сегменте \а, /3\, тогда найдется точка 7 на этом сегменте, такая, что
д(/3) - д(а)
/3 — а
=: ръ
Доказательство. Пусть точки разрыва {ск®}! производной функции д удовлетворяют неравенствам ос < си\ < а.2 • • • < щ < /3. Введем обозначения ао := а и :=
/3. Применяя последовательно на каждом интервале непрерывности производной формулу Ньютона-Лейбница и затем теорему о среднем, имеем:
\д{Р)-д{а)\ = \ [ д'{г])Аг]
9 9
ГОЦ+1 ГОЦ+1
^2 ЇШг, <Е/ №(л)\&Л =
__п 04 • п 04
г=о *
і=0 *
2^(оі+і - ач) \д'ы\ .
і=О
Если утверждение леммы было бы ложно, то есть \д'(Ь)\<р\ для всех Ь из (а, (3)\{оц}{, то из последней цепи неравенств следовало бы абсолютно ложное утверждение
I
\д(/3) - д{а) \ < ^ (04+1 - ац) /ц = (/3 - а) р1.
і=о
□
Для кусочно-гладкого сигнала я в терминах конечных разностей, согласно лемме 1, справедливо утверждение
|Д(г)
(37* Є (и,и + к)) |в'(7»)| ^
к
і = 0,1,...
Теорема 1. Пусть вещественная гладкая функция в имеет кусочно-непрерывную вторую производную на сегменте [^-1,^+1], тогда существует т* € (и, и + К) такое, что
|Д2(»)|
к2
=: р2, г = 0,1,...
Доказательство. Доопределяя производную в' в точках и произвольными конечными значениями, распишем конечную разность второго порядка в момент
|Д2(г)| = |я*+1 - 2Si + в»_1| =
= — / я/(4)си+ / я/(4)(И
I «'йг — 1 ''ti
/ [/(і +/і) —/(і)] (ІІ ** Іі — і
'І/і — 1 гі
[ І/(і +/і) —/(і)! (І і
**Іі — 1
— 1
где 7(і) есть некоторый момент из интервала (і,і + К) (см. лемму 1).
Пусть теорема 1 неверна, то есть ^(і)! < Р2 ДЛЯ всех І ИЗ (іг_і,іг+і) За ИСКЛЮЧЄНИЄМ точек разрыва второй производной. Тогда из последней цепи неравенств следует абсолютно ложное утверждение
IД2(г)| < / Ьр2АЬ = Н2Р2 = |д2(г)| .
**Іі — 1
□
0
Теорема 2. Передискретизация с увеличением частоты дискретизации в К раз с помощью линейной интерполяции снижает в К раз в окрестностях исходных узлов уровень гладкости высшего по уровню гладкости ассоциированного сигнала.
Доказательство. Действительно, в силу (6) конечная разность второго порядка после дискретизации А2(г) меньше в К раз, чем А2(г). Но, согласно теореме 1, восстановленный сигнал будет иметь в некоторой точке интервала [и~г,*-, вторую производную, по модулю не меньшую, чем
\А1(г)\/62 = \А1(г)\К2 = \А2(г)\К.
□
Оптимизация дополнительных сэмплов
От упомянутых низкочастотных шумов, появляющихся при линейной интерполяции, можно попытаться избавиться более равномерным распределением разностных производных второго порядка по дополнительным узлам. Для этих целей здесь будут рассмотрены КЭК-сплайны дефекта 1.
Поскольку они имеют непрерывную производную в исходном узле, при получении конечных разностей второго порядка в окрестности исходного узла линейные части прилегающих звеньев взаимно уничтожаются, и достаточно найти отклонения этих двух звеньев от касательной в ближайшем справа дополнительном узле. Отклонение прилегающего правого звена Ог (6) в узле Ьц =1 + 5 можно получить решением задачи Коши следующего вида:
- С[М+1]~б + С[М+1](°)-9Г “
Щ' = С'" = ЗЯг+г - ЭЯг-Ц + Эя* - Зя*_1, £>"(+0) = С"(+0) =
— — йг+2 + + 2$г_1,
Щ(0) = Щ0) = 0.
Оно элементарно:
А(*)= [6С[^+1](0«И =
и 0
=а
=[ГЛ[с'^+с,^о))а™=
= ГГ (с[м+чг + с[м+ч(°)) <^ =
= 1о (С^+1] ~2 + С?.<+1](°)*) =
— ($г+2 — Зв*-(-1 + ЗЯг — в*—1) — +
62
+ (—вг+2 + + 2вг-1) — .
Последнее выражение будет удобно представить с помощью конечных разностей:
= (Д2(г + 1) - Д2(>)) |+
+ (-д2(; + 1) + гд2(!))
Для получения аналогичным образом отклонения влево -ОД—6) вычислим значение второй производной левого прилегающего звена в г-м узле:
Д?(-0) = Д£_1(+0)+ [1о"1<и =
= —23г+1 + 5Яг — 4$г_1 + 2 + (ЗЗг+1 — 93г +
+ 93г_1 - Зв^2)Ь
= ^г+1 — + 53г_1 — вг-2-
Отсюда, решая аналогичную задачу Коши с другими производными второго и третьего порядков на обратном времени, получаем
53
^(—6) = (—+ Звг — ЗЯг-! + +
82
+ {2вг+1 - 5Si + - 3*-2)у •
Следовательно, конечная разность второго порядка в г-м узле такова:
А2(г)=А(5) + А(-5) =
г3
= (я*+2 - + бвг - Авг-! + вг-2) — +
82
+ (-3*+2 + 6йг+1 - 10в* + бз*-! - вг-2)у =
= (Д2(г + 1) - 2Д2(г) + Д2(» - 1)) | +
+ (-Д2(; + 1)+4Д2й-Д2(г-1))|. (7)
Допустим теперь, что пополненному потоку сэмплов соответствует имеющая непрерывную вторую производную функция /, тогда, как известно, на интервале [г — 1 + >с8, Ьц] = [г —
6, г + 5] должен найтись момент г, такой, что
/Г(т) =
А2(г)
52
= (Д2(г+1)-2Д2(г) + Д2(г-1))| + + (-Д2(г + 1) + 4Д2(г) - Д2(г - 1)) ^ = Д3(г - 1)5 - Д2(г + 1) + 4Д2(г) - Д2(г - 1)
“ 2 ' Аналогичным образом поступаем в дополнительных узлах. Так же, как и на стыке звеньев, получим уклонения от касательной за время 6 при прямом и обратном направлении. Поскольку теперь имеется только один полином третьей степени, то вклады в уклонения его члена третьей степени взаимно уничтожаются, а уклонения от члена второй степени равны друг другу и в сумме равны
Ан(Ъ') = 2 [ [ (6ад'(5 + 2Ъг) с1т (14 =
7о Jo = (6а^6+2Ь^62 = 3(Зг+2-ЗЗг+1+ЗЗг-Зг_1)^'53 +
+ (—йг+2 + — 5Яг + 2з1—{)52 =
= 3(Д2(г + 1)-Д2(г))^3 +
+ (—Д2(г + 1) + 2Д2(г)) 62. (8)
Оценки (7) и (8) помогают установить следующий результат:
Теорема 3. При передискретизации с увеличением частоты дискретизации в К раз с помощью интерполяции КЭК-сплайнами уровень гладкости восстанавливаемого сигнала снизится всюду не более чем в 3 раза независимо от К.
Доказательство. Конечные разности второго порядка (5) в исходных узлах с точностью до величин второго порядка малости1 относительно нового шага дискретизации 6 согласно (7) имеют вид
Д2(*)«(-Д2(* +1) + 4Д2(>) - д2(; -1))
отсюда
|Д»«|
<
<
(|Д2(г + 1)|+4|Д2(г)| + |Д2(г-1)|)|-«
^3 тах |Д2(г+7')| 82 = *€{-1,0,1}' 1
= -^2 таХ |Д2(^+^’)|-к *€{-1,0,1}
Согласно теореме 1, это означает, что уровень гладкости восстанавливаемого сигнала снизится не более чем в 3 раза на сегменте [и — 5,и + 5\. Проделаем аналогичную оценку уровня гладкости в окрестностях дополнительных узлов (4). Согласно (8),
|Д2(Ь)| и | —Д2(г + 1) + 2Д2(г)| 52 ^
^3 тах IД2 (г +/с) I <52
/г€{0,1} П
3 = 1,^,
то есть в окрестности дополнительных узлов конечная разность второго порядка увеличивается не более чем в 3 раза, что влечет снижение уровня гладкости там тоже не более чем в 3 раза. □
Теоремы 2 и 3 наглядно показывают, что интерполяция КЭК-сплайнами существенно выгоднее линейной интерполяции, поскольку предъявляет меньшие требования к ФНЧ.
Пример 1. Построим два последовательных звена КЭК-сплайнов. Пусть /1=1, = 1,8;
1,6, 0,3, 5^-(-2 О, 1 И ^^-(-3
— 1,5. Тогда на участке [г, г + 1] разностные производные равны пц = ггц+\ = —0,75, а сплайн получится С'[^+1](т) = 1,1т3—1,65т2—
0,75т + 1,6. На участке [г + 1, г + 2] имеем ГПг+2 = —0,9 И С[;+м+2](т) = -1,25т3 + 1,8т2-0,75т+ 0,3.
С помощью этих двух сплайнов производится передискретизация с повышением частоты в К = 4 раза. На рисунке 2 представлены два последовательных звена сплайнов, полученных линейной и квазиэрмитовой кубической интерполяцией на отрезке [г, г+2].
------- линейная интерполяция
кубичвская интбрполяция
], 75 г + 1 г -I- 1, 25 г + 1, 5 з Ч- 1, 75
Рис. 2. Линейная и квазиэрмитова кубическая интерполяции.
1Конкретно с погрешностью третьего порядка малости: 3 (Д2(г + 1) — Д2(г)) jS3
@
Вертикальные пунктирные линии соответствуют отсечкам сигнала после передискретизации.
На рис. 2 хорошо видно, что доставляемые КЭК-сплайном б'з^ значения выходного сигнала, в сравнении с линейным сплайном Бц, в исходных узлах дадут меньшие конечные разности второго порядка. В силу теоремы 1 это означает более высокий уровень гладкости при использовании 5^2 сравнительно с который благоприятнее для ФНЧ.
Исследование выпуклости и вогнутости КЭК-сплайнов
Возможно ли изменение выпуклости на вогнутость и наоборот в г-м звене? Ясно, что такому изменению соответствует точка перегиба
^г 1 &г+2 “Ь 2йг_1
За% 3 5^_|_2 Зз^_|_ х “Ь Зв^
Отсюда получаем необходимые условия смены выпуклости и вогнутости в г-м звене:
О < Яг+2 - 43г+1 + бя* - 2з*_1 <
<С 3(з^-(-2 35^_|_х “Ь Зв^ й^-х)
либо
О < -Яг+2 + 43г+1 - 5з* + 2з*_1 <
<С 3(3^_|_2 35^-(-1 “Ь Зв^
Что эквивалентно
/0 <С &г+2 А^г+1 “Ь 2Яг-1 /д\
10 < 23г+2 - 53г+1 + 4вг - ' '
либо
ГО > Бг+2 - 43г+1 + 5з* - 2Бг-\ ^
10 > 2з*+2 - бЯг-ц + ^ ' '
Объединяя системы (9) и (10), получаем (Зг+2 - 4з*+1 + - 2в*_!)•
•(23г+2 - 53г+1 + 4з* - 3*_1) > 0. (11)
1. Вырожденный случай ггц = т*+1 ■<=>■
'Ь,г+2 г г 1 г—1 ■ Тогда
1 &г+2 “Ь 2йг_1
^ ®г+2 Зй^-!-! “Ь Зв^
_ 1 $1+2 4в*_|_1 ~Ь 2з<1—1 __
^ ®г+2 35^_|_х "Ь
_ 1 -35г+1 + 6Si - 35г_1 _ 1
3 2й^-(-1 ~Ь 42—1 2
2. Невырожденный случай и выполняется условие (11).
Это означает, что тр Е (0,1) \ {1/2} и параллелограмм звена существует. Звено сплайна одной частью проходит внутри него, другой
— снаружи.
3. Условие (11) не выполняется.
Это влечет невырожденность. Следовательно, параллелограмм звена существует. Поскольку также верно тре(0,1), то звено сплайна либо строго выпукло, либо строго вогнуто для г € [0,1]. Что возможно только когда сплайн полностью находится внутри параллелограмма звена.
Заключение
Теорема 2 показывает, что получение дополнительных сэмплов в дополнительных узлах посредством линейной интерполяции приводит к увеличению вторых производных в окрестностях исходных узлов, в то время как интерполяция с помощью КЭК-сплайнов не порождает таких возмущений. Программная реализация интерполяционных алгоритмов на основе линейных и КЭК-сплайнов показала, что отсутствие таких возмущений благоприятно для качественного воспроизведения звука.
«Родные» частоты АЦП-ЦАП ZET 230 давали следующие возможности умножения частот: двух- и четырехкратное от 22050 Гц к 44100 Гц и к 88200 Гц, и трехкратное с 32000 Гц к 96000 Гц. Во всех трех случаях дополнительные сэмплы (1, 3, 2 соответственно) вставлялись в выходной файл с помощью интерполяции двух видов: линейной и квазиэр-митовой кубической. По сравнению с линейной интерполяцией при интерполяции КЭК-сплайнами получался более чистый звук, то есть в нем было снижено дребезжание и прочие шумы, что позволяло распознавать на слух выходной сигнал с субъективно меньшим напряжением и концентрацией внимания. Что и ожидалось, поскольку КЭК-сплайны в сравнении с линейными подготавливают выходной сигнал, лучше подходящий для дальнейшего аппаратного сглаживания в ЦАП.
Вполне возможно, что идея передискретизации с применением КЭК-сплайнов будет продуктивна и для повышения качества воспроизведения видеофайлов. Следует отметить, однако, что трудности ее реализации на порядки выше, чем в случае аудиофайла, поскольку интерполяция тогда должна быть векторной величины двумерного аргумента.
0
Литература
1. Завьялов Ю. С., Квасов Б. И., Мирошниченко В. Л. Методы сплайн-функций. М.: Наука, 1980. 352 с.
2. Камачкин А. М., Михеев С. Е., Евстафьева В. В. Модели колебаний в нелинейных системах. СПб.: Изд-во СПбГУ, 2004. 194 с.
3. Михеев С. Е. Об одном парадоксе в теоремах о методе Ньютона // Вестник СПбГУ. Сер. 10. 2013. № 1. С. 22-36.
4. Akima Н. A New Method of Interpolation and Smooth Curve Fitting Based on Local Procedures // Journal of the Association for Computing Machinery. Vol. 17, 1970. P. 589-602.
СВЕДЕНИЯ ОБ АВТОРАХ:
Михеев Сергей Евгеньевич
профессор, д. ф.-м. н.
Санкт-Петербургский государственный университет факультет прикладной математики - процессов управления,
Университетский пр., 35, Старый Петергоф, Санкт-Петербург, Россия, 198504 эл. почта: him2@mail.ru тел.: 89602405627
Морозов Петр Дмитриевич аспирант каф. информационных систем Санкт-Петербургский государственный университет факультет прикладной математики - процессов управления,
Университетский пр., 35, Старый Петергоф, Санкт-Петербург, Россия, 198504 эл. почта: pm-morozovpd@yandex.ru тел.: 89046328358
5. Karup J. On a new mechanical method of graduation. Transactions of the Second International Actuarial Congress, London, 1899.
6. Lyons R. Understanding Digital Signal Processing, 2nd Edition. Prentice Hall PTR Upper Saddle River, NJ, USA, 2004. 688 p.
7. Schoenberg I. J. Contributions to the Problem of Approximation of Equidistant Data by Analytic Functions // Quarterly of Applied Mathematics. 1946. Vol. 4, nos. 1, 2. P. 45-99, 112-141.
8. Sprague T. B. Explanation of a new formula for interpolation // Journal of the Institute of Actuaries. 1880. Vol. 22. P. 270-285.
Mikheev, Sergey
Saint-Petersburg State University
Faculty of Applied Mathematics & Control Theory
35 Universitetskiy St., 198504 Saint-Petersburg, Russia
e-mail: him2@mail.ru
tel.: 89602405627
Morozov, Petr
Saint-Petersburg State University
Faculty of Applied Mathematics & Control Theory
35 Universitetskiy St., 198504 Saint-Petersburg, Russia
e-mail: pm-morozovpd@yandex.ru
tel.: 89046328358
