Цифровая передискретизация сигналов на основе полиномиальной интерполяции. Фильтр Фарроу

Содержание

DSPL-2.0 — свободная библиотека алгоритмов цифровой обработки сигналов

Распространяется под лицензией LGPL v3

Страница библиотеки на GitHub

Обнаружили ошибку? Выделите ее мышью и нажмите ctrl+enter
Введение

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

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

Также бывает необходимо произвести преобразование частоты дискретизации в дробное число раз \frac{P}{Q}, где P и Q — целые числа.

В данном разделе мы рассмотрим один из эффективных методов решения указанных задач на основе полиномиальной Лагранжевой интерполяции [1, стр. 115]. Также будет рассмотрен фильтр Фарроу [2], реализующий метод полиномиальной Лагранжевой интерполяции для цифровой передискретизации сигналов [3, 4].

Постановка задачи

Пусть имеется входной сигнал s(n), n = 0,1,2 \ldots, отсчеты которого взяты с частотой дискретизации F_\textrm{s} Гц. Тогда отсчет с номером n соответствует моменту времени t_n = \frac{n}{F_\textrm{s}} с.

Необходимо произвести пересчет цифровых значений входного сигнала s(n) в цифровые значения сигнала y(k), k = 0,1,2 \ldots, взятые с частотой дискретизации F_\textrm{y} = \frac{P}{Q}  F_\textrm{s}, где P и Q целые числа[1]. Кроме того, первый отсчет y(0) должен быть задержан во времени относительно первого отсчета входного сигнала s(0) на величину \Delta t = \frac{x_0}{F_\textrm{s}} с, -1 \leq x_0 < 0.

Далее в этом разделе переменной n будут индексироваться отсчеты исходного сигнала s(n), взятые с частотой дискретизации F_\textrm{s}, а переменной k будут индексироваться отсчеты сигнала после передискретизации y(k).

Таким образом, k-ый отсчет y(k) соответствует моменту времени:

equation 1
(1)
Произведя нормировку к частоте дискретизации F_\textrm{s} = 1 Гц, можно записать:
equation 2
(2)
Временны́е соотношения моментов взятия отсчетов сигналов s(n) и y(k) при различных значениях P, Q и x_0, показаны на рисунке 1.

Временны́е соотношения моментов взятия
отсчетов сигналов  и , при различных значениях ,  и
Рисунок 1. Временны́е соотношения моментов взятия отсчетов сигналов s(k) и y(n), при различных значениях P, Q и x_0

На рисунке 1 показаны моменты дискретизации исходного сигнала s(n) (шкала n), а также моменты взятия отсчетов передискретизированного сигнала y(k) (шкала k).

Частный случай 1. P = Q = 1, x_0 = 0. В этом случае частота дискретизации F_\textrm{y} = F_\textrm{s}, и сигнал y(k) полностью повторяет s(n).

Частный случай 2. P = Q = 1, x_0 = 0.2. В этом случае частота дискретизации F_\textrm{y} = F_\textrm{s}, но сигнал y(k) имеет дробнную задержку относительно входного сигнала s(n). Заметим, что y(k) хранит задержанное на x_0 = 0.2 значение входного сигнала. Это означает, что шкала y(0) соответсвует моменту времени на x_0 = 0.2 отсчета раньше чем s(0). Это объясняет знак «минус» в выражениях (1) и (2).

Частный случай 3. P > 1,  Q = 1, x_0 = 0. В этом случае частота дискретизации сигнала y(k), F_\textrm{y} = P  F_\textrm{s}, в целое число раз P выше частоты дискретизации сигнала s(n). Таким образом, имеем цифровую интерполяцию сигнала s(n).

Частный случай 4. P > 1, Q > 1, P \neq Q, x_0 = 0. В этом случае имеем дробную передискретизацию сигнала s(n). Частота дискретизации сигнала y(k) равна F_\textrm{y} = \frac{P}{Q}  F_\textrm{s}.

Частный случай 5. P > 1, Q > 1, P \neq Q, x_0 > 0. В этом случае имеем дробную передискретизацию сигнала s(n) плюс добавление дробной задержки. Это наиболее общий из рассматриваемых частных случаев. Частота дискретизации сигнала y(k) равна F_\textrm{y} = \frac{P}{Q}  F_\textrm{s}.

Для решения задачи передискретизации необходимо по имеющемуся дискретному сигналу s(n), n = 0, 1,2, \ldots , произвести восстановление непрерывного сигнала s(t) и рассчитать его дискретные значения для моментов времени s(x_k), где x_k рассчитывается согласно (2).

Процесс передискретизации для частного случая 2 (компенсация дробной задержки) показан на рисунке 2.

Передискретизация сигнала
Рисунок 2. Передискретизация сигнала

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

Интерполяционный полином Лагранжа

В курсе математического анализа доказывается, что через N точек проходит единственный полином степени N-1. Например, через две точки можно провести только одну прямую, через три точки можно провести только одну параболу. Единственный полином степени N-1, проходящий через N точек, называется интерполяционным полиномом Лагранжа:

equation 3
(3)
где a_n — коэффициенты полинома, которые рассчитываются на основе дискретных значений s(n), n = 0 \ldots N-1.

Для расчета коэффициентов a_n мы можем составить систему линейных уравнений (смотри рисунок 2):

equation 4
(4)

Данная система может быть записана в матричной форме как:

equation 5
(5)
где
equation 6
(6)

equation 7
(7)

Решение системы уравнений может быть получено в виде:

equation 8
(8)
где \mathbf{M^{-1}} — обратная матрица для матрицы \mathbf{M}.

Процедура вычисления обратной матрицы является одной из наиболее затратных, с точки зрения вычислительных ресурсов. В общем случае требуется количество операций O(N)\sim N^3, пропорциональное кубу размерности квадратной матрицы. Например, для расчета коэффициентов кубического полинома для (N=4) в общем случае требуется O(4)\sim 4^3 = 64 умножения, что трудно реализуемо на практике. В дальнейшем мы увидим, как можно сократить количество операций умножения до трех, причем две из трёх операций будут умножения на \frac{1}{2}, что легко реализуется в целочисленной арифметике путем битового сдвига на один разряд вправо.

Мы еще вернемся к данному вопросу ниже. Сейчас же мы рассмотрим эффективное представление полиномов по схеме Горнера (Horner structure).

Эффективное представление полиномов по схеме Горнера

Рассмотрим полином (3) более подробно:

equation 9
(9)
Вынесем t за скобки, получим:
equation 10
(10)
Снова вынесем t за скобки, получим:
equation 11
(11)
Таким образом, вынося t за скобки возможное число раз, получим множество вложенных скобок:
equation 12
(12)
Например, для N = 4 получим кубический полином:
equation 13
(13)
который может быть реализован в виде структуры, показанной на рисунке 3.

Расчет кубического полинома по схеме Горнера
Рисунок 3. Расчет кубического полинома по схеме Горнера

Данное представление называется схемой Горнера (Horner structure) и широко используется для эффективного расчета значений полиномов [5, стр. 126], [6, стр. 144]. Ее основное преимущество в том, что для расчета значения полинома степени N-1, при известных коэффициентах, требуется N-1 умножитель и сумматор.

Использование полиномиальной интерполяции для цифровой передискретизации сигналов

Построение полинома для интерполяции часто проводят при использовании ортогональных полиномов Лагранжа. Однако решение системы линейных уравнений (4) является более общим подходом, который может быть применен не только для Лагранжевой интерполяции, но и для других методов полиномиальной интерполяции (Эрмитова полиномиальная интерполяция, сплайн-интерполяция и других).

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

Наиболее часто используемым методом является кусочно-полиномиальная кубическая интерполяция ввиду компромисса между точностью и вычислительными затратами. Рассмотрим метод кусочно-полиномиальной кубической интерполяции более подробно.

Полином третьей степени может быть построен по четырем точкам согласно выражению (8).

Для любого 1 \leq x_k < N-2 можно получить два отсчета правее и левее x_k как это показано на рисунке 4a. Тогда значение y(k) может быть рассчитано на основе четырех входных отсчетов s(n), s(n-1), s(n-2), s(n-3).

Кусочно-кубическая интерполяция
Рисунок 4. Кусочно-кубическая интерполяция

Обратим внимание, что матрица \mathbf{M} зависит только от индексов входных отсчетов сигнала. Значит, чтобы не пересчитывать матрицу \mathbf{M} и обратную \mathbf{M^{-1}} для каждого нового k, мы можем зафиксировать индексы по оси абсцисс. Пусть x_k = n-1-\Delta_k (cм. рисунок 4а), тогда вместо x_k мы будем использовать значение  -1 \leq -\Delta_k < 0, как это показано на рисунке 4б. В этом случае мы можем один раз расчитать обратную матрицу \mathbf{M^{-1}} и использовать ее для любого x_k.

Система уравнений для расчета коэффициентов полинома при фиксированных индексах, согласно рисунку 4б, имеет вид:

equation 14
(14)
Тогда мы можем представить матрицу \mathbf{M} и получить обратную ей \mathbf{M^{-1}}:
equation 15
(15)

Мы не приводим процедуру расчета обратной матрицы. Вы можете самостоятельно убедиться в том, что \mathbf{M^{-1} M = I}, где \mathbf{I} — единичная матрица.

Коэффициенты кубического полинома могут быть получены как результат умножения обратной матрицы \mathbf{M^{-1}} на вектор отсчетов входного сигнала \mathbf{s}:

equation 16
(16)
Произведем матричное умножение (16) и получим выражения для коэффициентов полинома:
equation 17
(17)
Из выражения (17) можно заметить, что коэффициенты полинома зависят от задержанных значений входного сигнала s(n). Мы можем трактовать формулы (17) для расчета коэффициентов как разностные уравнения КИХ-фильтра. Тогда расчет каждого коэффициента будет выполнять свой КИХ-фильтр, и структурная схема фильтра передискретизации может быть представлена, как это показано на рисунке 5.

Каждый КИХ-фильтр рассчитывает один коэффициент полинома, после чего они подаются на схему Горнера для расчета полинома при текущем значении дробной задержки \Delta_k. Структура, приведенная на рисунке 5, носит название фильтра Фарроу (Farrow filter) [2].

Структура фильтра Фарроу для цифровой передискретизации  сигналов на основе 
полиномиальной интерполяции Лагранжа
Рисунок 5. Структура фильтра Фарроу для цифровой передискретизации сигналов на основе полиномиальной интерполяции Лагранжа

Важно отметить, что выход фильтра Фарроу y(k) отстает от входного сигнала s(n) на один отсчет в единицах отсчетов входного сигнала (при частоте дискретизации F_{\textrm{s}}).

Использование фиксированных значений оси абсцисс позволяет избавиться от необходимости обращения матрицы в реальном времени. Однако фильтр Фарроу также требует значительного количества умножителей на \frac{1}{6}, \frac{1}{3} и \frac{1}{2}. При этом, в случае реализации в целочисленной арифметике, умножения на \frac{1}{2} можно считать тривиальными, потому что они реализуются как сдвиг на один разряд вправо.

В следующем параграфе мы модифицируем фильтр Фарроу для уменьшения количества умножителей при расчете коэффициентов полинома.

Минимизация умножителей при расчете интерполяционного полинома Лагранжа

При реализации фильтров передискретизации сигналов мы должны стремиться минимизировать количество умножителей при расчете коэффициентов полинома.

Преобразуем выражение (17). Коэффицент a_3 можно представить в виде:

equation 18
(18)
В выражении (17) также можно заметить, что
equation 19
(19)
откуда можно рассчитать коэффициент a_1 как:
equation 20
(20)
Для расчета коэффициента a_2 можно использовать следующее следствие из (17):
equation 21
(21)
тогда:
equation 22
(22)

Таким образом, мы можем переписать (17) в оптимизированном виде всего с тремя умножителями (причем два из них равны \frac{1}{2}):

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

Структурная схема оптимизированного фильтра цифровой передискретизации
Рисунок 6. Структурная схема оптимизированного фильтра цифровой передискретизации

Таким образом, мы получили алгоритм, который рассчитывает коэффициенты интерполяционного полинома Лагранжа с использованием всего одного умножителя на \frac{1}{6} и двух тривиальных умножителей на \frac{1}{2}.

Выводы

В данном разделе мы рассмотрели структурные схемы фильтров цифровой передискретизации на основе полиномиальной Лагранжевой интерполяции.

Была рассмотрена схема Горнера для расчета значений полинома. Количество умножителей и сумматоров схемы Горнера равно порядку полинома.

Получена структура фильтра Фарроу, и произведена его оптимизация. В результате коэффициенты интерполяционного полинома Лагранжа могут быть рассчитаны при использовании всего трех умножителей. Требуется одно умножение на \frac{1}{6} и два тривиальных умножения на \frac{1}{2}.

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

Примечания

[1] Все приведенные выражения также справедливы для произвольных отличных от нуля вещественных значений P и Q, в том числе иррациональных.

Список литературы
[1] Каханер Д., Моулер К., Нэш С. Численные методы и математическое обеспечение. Москва, Мир, 1998.

[2] Farrow C.W. A Continuously Variable Digital Delay Element. Circuits and Systems, IEEE International Symposium. 1988, p. 2641–2645. vol. 3

[3] Gardner F. Interpolation in Digital Modems-Part I: Fundamentals: IEEE Transactions on Communications, Vol. 41, No. 3, March 1993, P. 501-507.

[4] Erup L., Gardner F. Interpolation in Digital Modems-Part II: Implementation and Performance: IEEE Transactions on Communications, Vol. 41, No. 6, June 1993, p.998-1008.

[5] Макконелл Дж. Основы современных алгоритмов. Москва, Техносфера, 2004, 368с. ISBN 5-94836-005-9.

[6] Блейхут Р. Быстрые алгоритмы цифровой обработки сигналов. Москва, Мир, 1989, 448 c.

[7] Franck A. Efficient Algorithms for Arbitrary Sample Rate Conversion with Application to Wave Field Synthesis. PhD thesis. Universitätsverlag Ilmenau, Ilmenau, 2012.

Последнее изменение страницы: 14.11.2020 (11:21:52)
Страница создана Latex to HTML translator ver. 5.20.11.14