Радиотехника и электроника, 2019, T. 64, № 2, стр. 127-135

ВВЕДЕНИЕ

Задача дифференцирования сигнала с ограниченным спектром при его корректной дискретизации имеет единственное точное решение (идеальный дифференциатор) в классе глобальных методов, которое может быть получено при дифференцировании интерполяционного ряда Котельникова (или, иным путем, на основе свойства дифференцирования для преобразования Фурье). Однако, при практической реализации такого численного дифференцирования возникают трудности, связанные с тем, что функция Котельникова и ее производные не ограничены во времени и вычисления требуют учета бесконечного количества отсчетов сигнала. Любые подходы, связанные с каким-либо ограничением количества отсчетов, являются поиском приближенных решений рассматриваемой задачи в классе уже локальных методов.

Наиболее распространенные традиционные подходы к задаче численного дифференцирования условно можно разделить на две группы. Первая предполагает предварительное решение задачи аппроксимации или интерполяции с последующим дифференцированием [1–6]. (Эта группа включает и конечные разности ввиду общности результатов.) Этот подход обобщен в данной работе на случай применения локальной интерполяции обобщенными многочленами, рассматриваются частные случаи, в рамках которых численное дифференцирование доводится до уровня практической реализации, показывается, что интерполяционный базис может быть выбран и неудачно в смысле обеспечения сходимости.

Ко второй группе можно отнести синтез оптимальных нерекурсивных широкополосных дифференциаторов [7, 8]. Оптимальность означает минимизацию ошибки аппроксимации частотной характеристики идеального дифференциатора на всем ее периоде. При этом само значение ошибки не нормируется, а приемлемые результаты достигаются при больших порядках дифференциаторов.

В данной работе рассмотрен численно-аналитический оптимальный синтез дифференциаторов небольшого порядка при заданных ограничениях на допустимые искажения, обеспечивающих максимальную граничную частоту. Такое численное дифференцирование требуется, например, при решении задач локальной сплайновой интерполяции сигналов [3, 5, 9, 10]. Полученные результаты апробированы автором при реализации программ и модулей, описанных в работе [11].

1. АНАЛИЗ СХОДИМОСТИ РЕГУЛЯРНОГО ЛОКАЛЬНОГО МЕТОДА ЧИСЛЕННОГО ДИФФЕРЕНЦИРОВАНИЯ СИГНАЛА

Будем полагать, что сетка дискретизации сигнала равномерна и бесконечна, узлы сетки ${{t}_{n}} = nT,$ $n = 0,$ $ \pm 1,\,\, \pm 2,...,$ где $T$ – период дискретизации сигнала. Под численным дифференцированием будем понимать задачу нахождения отсчетов k-й производной ${{s}^{{(k)}}}({{t}_{n}})$ в моменты дискретизации ${{t}_{n}}$ сигнала $s(t)$ со спектром, ограниченным частотой ${{\omega }_{m}},$ по его дискретным значениям $\left\{ {{{s}_{l}} = s({{t}_{l}})} \right\}_{{l = n - M}}^{{n + M}},$ M – порядок метода численного дифференцирования. Записывая общее выражение для локального метода численного дифференцирования, связывающее ${{s}^{{(k)}}}({{t}_{n}})$ и $\left\{ {{{s}_{l}}} \right\}_{{l = n - M}}^{{n + M}},$ сохраним обязательную линейность, нерекурсивность и регулярность, присущие точному (глобальному) решению задачи:

где $\left\{ {d_{m}^{{(k)}}} \right\}_{{m = - M}}^{M}$ – коэффициенты алгоритма численного дифференцирования сигнала, независимые от вида сигнала и момента дискретизации. Линейность здесь означает, что численной оценке производной линейной комбинации сигналов будет соответствовать такая же линейная комбинация из отдельных численных оценок производных каждого из сигналов. Нерекурсивность означает, что численная оценка отсчета производной в данный момент дискретизации не зависит от таковых оценок в другие моменты дискретизации. Регулярность означает, что в каждый момент дискретизации численное дифференцирование выполняется по одному и тому же правилу.

В ряде случаев коэффициенты алгоритма численного дифференцирования удобно находить для единичного периода дискретизации. Обозначив их при этом $\tilde {d}_{m}^{{(k)}},$ запишем

Алгоритм (1) описывает нерекурсивный цифровой фильтр с нарушенной причинностью с входным сигналом $\left\{ {{{s}_{n}}} \right\}_{{n = - \infty }}^{{ + \infty }},$ выходным сигналом $\left\{ {\psi _{n}^{{(k)}}} \right\}_{{n = - \infty }}^{{ + \infty }}$ и коэффициентами $\left\{ {d_{m}^{{(k)}}} \right\}_{{m = - M}}^{M}.$ При физической реализации алгоритма в реальном времени нарушение причинности может быть устранено введением соответствующего запаздывания. В рамках данной работы это не дает никаких преимуществ, поэтому запаздывание не вводится, что, однако, не запрещает рассматривать частотную характеристику $H_{{\text{д }}}^{{{\text{(}}k{\text{)}}}}(\omega )$ для алгоритма (1) как отношение преобразования Фурье дискретного сигнала $\left\{ {\psi _{n}^{{(k)}}} \right\}_{{n = - \infty }}^{{ + \infty }}$ к преобразованию Фурье дискретного сигнала $\left\{ {{{s}_{n}}} \right\}_{{n = - \infty }}^{{ + \infty }}\,:$

В дальнейшем, по аналогии с идеальной, частотная характеристика дифференциатора при нечетном порядке производной будет полагаться чисто мнимой с нечетно-симметричной мнимой частью, а при четном порядке производной – действительной четно-симметричной, а также $H_{{\text{д }}}^{{(k)}}(0) = 0,$ при этом и сами коэффициенты дифференциатора (2) удовлетворяют условиям симметрии:

Под сходимостью метода численного дифференцирования будем понимать его свойство осуществлять точное дифференцирование сигнала при неограниченном уменьшении периода дискретизации:

В рассматриваемом классе сигналов с учетом периодичности спектральных функций дискретных сигналов и частотной характеристики, а также свойства дифференцирования для преобразования Фурье, записанное условие сходимости можно переформулировать следующим образом:

Разложив экспоненту в (3) в ряд Маклорена, запишем:

Члены ряда с номерами $n > k$ никак не повлияют на значение предела, а для выполнения (5) требуется, чтобы

Полученные условия сходимости (6) дают систему из $k + 1$ линейных алгебраических уравнений c $2M + 1$ неизвестными. В частном случае $k = 2M$ эта система может быть единственным образом разрешена. Значения коэффициентов алгоритма численного дифференцирования, полученные в результате решения системы, при $M = 1,\,\,k = 2;$ $M = 2,\,\,k = 4;$ $M = 3,\,\,k = 6$ приведены в соответствующих строках табл. 1.

В другом частном случае, считая частотную характеристику дифференциатора чисто мнимой с нечетно-симметричной мнимой частью, положим $\tilde {d}_{0}^{{(k)}} = 0$ и количество неизвестных в (6) сократится до $2M,$ а система будет однозначно разрешаемой при $k = 2M - 1.$ Значения коэффициентов дифференциатора, полученные при решении системы в этом случае, при $M = 1,\,\,k = 1;$ $M = 2,\,\,k = 3;$ $M = 3,\,\,k = 5$ приведены в соответствующих строках табл. 1.

2. АНАЛИЗ ТОЧНОСТИ ПРИ ЧИСЛЕННОМ ДИФФЕРЕНЦИРОВАНИИ СИГНАЛА

В виду линейности численного дифференцирования можно отдельно рассматривать преобразование каждой гармонической составляющей сигнала, задаваемое (1).

При воздействии дискретного гармонического сигнала ${{s}_{n}} = U\cos (\omega {{t}_{n}} + \varphi )$ на идеальный дифференциатор откликом является гармонический сигнал той же частоты, что и воздействующий, с амплитудой пропорциональной частоте воздействия $U{{\left| {j\omega } \right|}^{k}}.$ В случае воздействия гармонического сигнала на дифференциатор с частотной характеристикой, отличной от идеальной, амплитуда отклика $U\left| {H_{{\text{д }}}^{{(k)}}(\omega )} \right|.$ В обоих случаях начальные фазы результирующих гармонических составляющих либо одинаковы, либо отличаются на π. Поэтому формально можно считать, что нежелательная разность между откликом дифференциатора и откликом идеального дифференциатора является гармоническим сигналом с амплитудой $U\left| {{{{(j\omega )}}^{k}} - H_{{\text{д }}}^{{(k)}}(\omega )} \right|.$

В качестве характеристики точности (показателя качества) дифференциатора будем рассматривать коэффициент искажений, определяемый как отношение действующего значения нежелательной разности между откликами идеального и неидеального дифференциаторов к действующему значению гармоники, получаемой при точном дифференцировании:

Введенный показатель качества (7) может быть интерпретирован и как относительная ошибка аппроксимации частотной характеристики идеального дифференциатора.

Коэффициент искажений удобно рассматривать как функцию нормированной частоты $\omega T.$ В большинстве случаев для методов дифференцирования, удовлетворяющих условиям сходимости (6), можно рассматривать граничную частоту ${{(\omega T)}_{{{\text{г р }}}}},$ характеризующую интервал $[0,{{(\omega T)}_{{{\text{г р }}}}}],$ в пределах которого коэффициент искажений не превышает некоторого порогового значения:

Анализ коэффициента искажений (8) может быть положен в основу выбора метода численного дифференцирования и/или запаса по частоте дискретизации: следует обеспечить, чтобы значение коэффициента искажений не превышало некоторого допустимого уровня в пределах полосы частот, где локализован спектр сигнала ${{\omega }_{m}}T < {{(\omega T)}_{{{\text{г р }}}}}.$

3. ЧИСЛЕННОЕ ДИФФЕРЕНЦИРОВАНИЕ НА ОСНОВЕ ЛОКАЛЬНОЙ РЕГУЛЯРНОЙ ИНТЕРПОЛЯЦИИ ОБОБЩЕННЫМИ МНОГОЧЛЕНАМИ

Интерполяцию будем рассматривать как преобразование $\psi (t) = {{I}_{n}}\left\{ {s(nT)} \right\},$ которое отображает множество дискретных сигналов вида $\left\{ {{{s}_{n}} = s(nT)} \right\}_{{n = - \infty }}^{{ + \infty }},$ полученных в результате дискретизации непрерывных сигналов $s(t)$, на множество непрерывных сигналов (интерполирующих функций) вида $\psi (t)$ и обладает свойствами согласованности и линейности, то есть:

где ${{k}_{1}},{{k}_{2}}$ – произвольные, в общем случае комплексные, числа.

С учетом (10) для интерполирующей функции можно получить [10]:

где ${{\varphi }_{n}}(t) = {{\varphi }_{0}}(t - nT) = {{I}_{k}}\left\{ {{{\delta }_{{k - n}}}} \right\}$ представляет собой результат интерполяции единичного отсчета, соответствующего моменту дискретизации ${{t}_{n}} = nT$ и называется базисной функцией, ${{\varphi }_{0}}(t)$ – порождающая функция, которую можно получить, построив непосредственно интерполирующую функцию для единичного отсчета, соответствующего нулевому моменту дискретизации:

В виду условия (9) интерполяционный базис $\left\{ {{{\varphi }_{n}}(t)} \right\}_{{n = - \infty }}^{{ + \infty }}$ фундаментален [10].

Рассмотрим случай, когда при локальной интерполяции на каждом частном интервале $t \in [{{t}_{n}} - 0.5T,{{t}_{n}} + 0.5T]$ фрагмент ${{\psi }_{n}}(t)$ интерполирующей функции совпадает с многочленом, график которого проходит через M предшествующие и M последующие относительно момента дискретизации ${{t}_{n}}$ отсчеты ${{s}_{{n - M}}},...,{{s}_{{n - 1}}},{\text{ }}$ ${{s}_{n}},{{s}_{{n + 1}}},{{s}_{{n + 2}}},...,{{s}_{{n + M}}},$ где $M \geqslant 0$ – порядок интерполяции. Интерполирующий фрагмент описывается выражением

где $\left\{ {{{w}_{k}}(t)} \right\}_{{k = 0}}^{{N - 1}}$ – линейно-независимая система функций (первичный базис), ${{C}_{{n,k}}}$ – коэффициенты, которые ввиду (10) линейно связанны с отсчетами сигнала: ${{a}_{{m,k}}}$ – постоянные, определяемые методом интерполяции и одинаковые для каждого частного интервала дискретизации.

На рис. 1 показан пример формирования интерполяционного фрагмента при $M = 2$ (штриховая кривая соответствует вспомогательным построениям, сплошная кривая – полученный результат).

Рассматривая восстановление единичного отсчета в соответствии с (12) и подставляя при этом (14) в (13), для порождающей функции получим:

а коэффициенты для (14) и (15) являются элементами матрицы вида

определяемой как:

где E – единичная матрица. Указанные коэффициенты существуют и единственны, когда $N = 2M + 1,$ а матрица W имеет отличный от нуля определитель.

Определив матрицу (16) и порождающую функцию (15), интерполирующий фрагмент (13) сможем представлять в виде:

Порождающая функция и интерполирующая кривая в рассматриваемом случае могут иметь разрывы, абсциссы которых соответствуют серединам интервалов дискретизации, однако в моменты дискретизации являются непрерывными и дифференцируемыми, что может быть использовано для оценки производных сигнала. Производная сигнала k-го порядка в момент дискретизации ${{t}_{n}}$ полагается совпадающей с производной интерполирующей функции:

где $d_{m}^{{(k)}} = \varphi _{0}^{{(k)}}( - mT)$ – коэффициенты алгоритма численного дифференцирования.

Коэффициенты алгоритма оценки производной, соответствующие ${{w}_{k}}(t) = {{t}^{k}},k = 0,...,N - 1,$ приведены в табл. 1 для $M = 1,2,3.$ Можно убедиться, что приведенные коэффициенты удовлетворяют условиям сходимости (6).

Коэффициенты, соответствующие первичному тригонометрическому базису

приведены в табл. 2 для $M = 1,2,3.$ Условия сходимости для них не выполняются, в некоторых случаях равенства в (6) имеют место с точностью не хуже 10^–2.

4. ОПТИМАЛЬНЫЕ ЛОКАЛЬНЫЕ РЕГУЛЯРНЫЕ МЕТОДЫ ЧИСЛЕННОГО ДИФФЕРЕНЦИРОВАНИЯ СИГНАЛОВ

В задаче оптимизации численного дифференцирования будем предполагать поиск коэффициентов дифференциатора $\left\{ {d_{m}^{{(k)}}} \right\}_{{m = - M}}^{M},$ удовлетворяющих условиям симметрии (4), сходимости (6) и обеспечивающих наибольшую граничную частоту ${{\left( {\omega T} \right)}_{{{\text{г р }}}}}$ при заданном допустимом уровне искажений ${{K}_{{{\text{и }}\,\,{\text{д о п }}}}},$ порядке метода численного дифференцирования M и порядке производной k.

Поставленная задача является многопараметрической однокритериальной задачей оптимизации. Ее решение потребовало разработки программного обеспечение большого объема, осуществлялось методами численной оптимизации, с использованием алгоритмов случайного поиска, в том числе генетического алгоритма. Поскольку решение получено численно, то в дальнейшем речь идет о квазиоптимальных дифференциаторах. Полученные результаты представлены в табл. 3, 4. В табл. 5 приведены граничные частоты квазиоптимальных дифференциаторов и дифференциаторов, полученных на основе локальной полиномиальной интерполяции.

Зависимости коэффициента искажений от нормированной частоты показаны на рис. 2 для первой, второй и третьей производной. Штриховые кривые соответствуют алгоритмам на основе локальной полиномиальной интерполяции (табл. 1), сплошные – квазиоптимальным алгоритмам (табл. 3, 4). Цифры возле кривых обозначают порядок алгоритма численного дифференцирования M.

ЗАКЛЮЧЕНИЕ

Сформулированы условия сходимости алгоритма численного дифференцирования: сходимость алгоритма численного дифференцирования, понимаемая как его свойство осуществлять точное дифференцирование сигнала при неограниченном уменьшении периода дискретизации, имеет место при выполнении условий (6).

Качество численного дифференцирования сигнала может быть охарактеризовано коэффициентом искажений (7), определяемым как отношение действующего значения нежелательной разности между откликами идеального и неидеального дифференциаторов к действующему значению гармоники, получаемой при точном дифференцировании.

Установлено, что свойства дифференциаторов сигналов, получаемых на основе методов локальной интерполяции зависят от удачного выбора базиса. Так, например, при выборе локального кусочно-полиномиального интерполяционного базиса условия сходимости обеспечиваются, а при выборе локального тригонометрического базиса – нет.

Показано, что оптимальный синтез дифференциаторов сигналов может осуществляться напрямую, путем численного решения многопараметрической оптимальной задачи, инкапсулирующей и требования сходимости.

Задача численного дифференцирования сигнала с ограниченным спектром может быть решена локально с наперед заданной точностью в некотором ограниченном диапазоне частот. Получены коэффициенты квазиоптимальных дифференциаторов для оценки 1–4-й производных сигнала при порядке алгоритмов 1–4. Найденные дифференциаторы обладают максимальной граничной частотой коэффициента искажений при заданном допустимом уровне искажений. Численное дифференцирование сигнала, в указанных случаях, доведено до уровня практической реализации.