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

Содержание
Введение
В современных системах передачи информации необходимо обеспечить тактовую синхронизацию передатчика и приемника. При этом использование независимых задающих генераторов на передающей и принимающей стороне приводит к разнице частот дискретизации сигналов передатчика и приемника.
Помимо этого, часто встает задача компенсации, так называемой дробной задержки на величину меньше интервала дискретизации. Фильтр, способный компенсировать произвольную дробную задержку, можно использовать и для синхронизации задающих генераторов ввиду их медленной нестабильности. Для этого необходимо использовать дополнительный детектор временной ошибки для расчета дробной задержки сигналов.
Также бывает необходимо произвести преобразование частоты дискретизации в дробное число раз $\frac{P}{Q}$, где $P$ и $Q$ — целые числа.
В данном разделе мы рассмотрим один из эффективных методов решения указанных задач на основе полиномиальной Лагранжевой интерполяции [1]. Также будет рассмотрен фильтр Фарроу, реализующий метод полиномиальной Лагранжевой интерполяции для цифровой передискретизации сигналов [2-5].
Постановка задачи
Пусть имеется входной сигнал $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} \cdot F_\textrm{s}$, где $P$ и $Q$ целые числа. Кроме того, первый отсчет $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)$ соответствует моменту времени:
$$ t_k = k \cdot \frac{Q}{P \cdot F_\textrm{s}} - \frac{x_0}{F_\textrm{s}}. $$
(1)
Произведя нормировку к частоте дискретизации $F_\textrm{s} = 1$ Гц, можно записать:
$$ x_k = t_k \cdot F_\textrm{s} = k \cdot \frac{Q}{P} - x_0. $$
(2)
Временные соотношения моментов взятия отсчетов сигналов $s(n)$ и $y(k)$ при различных значениях $P$, $Q$ и $x_0$,показаны на рисунке 1.

Рисунок 1. Временные соотношения моментов взятия отсчетов сигналов $s(k)$ и $y(n)$,
при различных значениях $P$, $Q$ и $x_0$
Черным показаны моменты дискретизации исходного сигнала $s(n)$, а красным — моменты взятия отсчетов передискретизированного сигнала $y(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)$.
Частный случай 3. $P > 1$, $ Q = 1$, $x_0 = 0$. В этом случае частота дискретизации сигнала $y(k)$, $F_\textrm{y} = P \cdot F_\textrm{s}$, в целое число раз $P$ выше частоты дискретизации сигнала $s(n)$. Таким образом,б имеем цифровую интерполяцию сигнала $s(n)$.
Рисунок 2. Передискретизация сигнала
Частный случай 4. $P > 1$, $Q > 1$, $P \neq Q,$ $x_0 = 0$. В этом случае имеем дробную передискретизацию сигнала $s(n)$. Частота дискретизации сигнала $y(k)$ равна $F_\textrm{y} = \frac{P}{Q} \cdot F_\textrm{s}$.
Частный случай 5. $P > 1$, $Q > 1$, $P \neq Q,$ $x_0 > 0$. В этом случае имеем дробную передискретизацию сигнала $s(n)$ плюс добавление дробной задержки. Это наиболее общий из рассматриваемых частных случаев. Частота дискретизации сигнала $y(k)$ равна $F_\textrm{y} = \frac{P}{Q} \cdot F_\textrm{s}$.
Для решения задачи передискретизации необходимо по имеющемуся дискретному сигналу $s(n)$, $n = 0, 1,2, \ldots $, произвести восстановление непрерывного сигнала $s(t)$ и рассчитать его дискретные значения для моментов времени $s(x_k)$, где $x_k$ рассчитывается согласно (2).
Процесс передискретизации для частного случая 2 (компенсация дробной задержки) показан на рисунке 2.
Поскольку мы оперируем только с индексами отсчетов сигнала, то и значения $x_k$ также задаются не в абсолютных значениях времени, а нормированными к частоте дискретизации.
Интерполяционный полином Лагранжа
В курсе математического анализа доказывается, что через $N$ точек проходит единственный полином степени $N-1$. Например, через две точки можно провести только одну прямую, через три точки можно провести только одну параболу. Единственный полином степени $N-1$, проходящий через $N$ точек, называется интерполяционным полиномом Лагранжа:
$$ s(t) = \sum_{n = 0}^{N-1}a_n \cdot t^n, $$
(3)
где $a_n$ — коэффициенты полинома, которые рассчитываются на основе дискретных значений $s(n)$, $n = 0 \ldots N-1$.
Для расчета коэффициентов $a_n$ мы можем составить систему линейных уравнений (смотри рисунок 2):
\begin{equation*} \begin{cases} a_0 + a_1 \cdot 0 + a_2 \cdot 0^2 + \ldots + a_{N-1} \cdot 0^{N-1} = s(0), \\ a_0 + a_1 \cdot 1 + a_2 \cdot 1^2 + \ldots + a_{N-1} \cdot 1^{N-1} = s(1), \\ a_0 + a_1 \cdot 2 + a_2 \cdot 2^2 + \ldots + a_{N-1} \cdot 2^{N-1} = s(2), \\ \vdots \\ a_0 + a_1 \cdot (N-1) + a_2 \cdot (N-1)^2 + \ldots + a_{N-1} \cdot (N-1)^{N-1} = s(N-1). \end{cases} \end{equation*}
(4)
Данная система может быть записана в матричной форме как:
\begin{equation*} \mathbf{M \cdot a = s}, \end{equation*}
(5)
где
\begin{equation*} \begin{array} & \mathbf{M} = \left[ \begin{array}{ccccc} 1 & 0 & 0^2 & \ldots & 0^{N-1}\\ 1 & 1 & 1^2 & \ldots & 1^{N-1}\\ 1 & 2 & 2^2 & \ldots & 2^{N-1}\\ \vdots& \vdots & \vdots &\ddots & \vdots\\ 1 & N-1 & (N-1)^2 &\ldots & (N-1)^{N-1} \end{array} \right], & \mathbf{a} = \left[ \begin{array}{ccccc} a_0 \\ a_1 \\ a_2 \\ \vdots \\ a_{N-1} \end{array} \right], & \mathbf{s} = \left[ \begin{array}{ccccc} s(0) \\ s(1) \\ s(2) \\ \vdots \\ s(N-1) \end{array} \right]. \end{array} \end{equation*}
(6)
Решение системы уравнений может быть получено в виде:
\begin{equation*} \mathbf{a = M^{-1} \cdot s}, \end{equation*}
(7)
где $\mathbf{M^{-1}}$ — обратная матрица для матрицы $\mathbf{M}$.
Процедура вычисления обратной матрицы является одной из наиболее затратных, с точки зрения вычислительных ресурсов. В общем случае требуется количество операций $O(N)\sim N^3$, пропорциональное кубу размерности квадратной матрицы. Например, для расчета коэффициентов кубического полинома для $(N=4)$ в общем случае требуется $O(4)\sim 4^3 = 64$ умножения, что трудно реализуемо на практике. В дальнейшем мы увидим, как можно сократить количество операций умножения до трех, причем две из них будут умножения на $\frac{1}{2}$, что легко реализуется в целочисленной арифметике путем битового сдвига на один разряд вправо.
Мы еще вернемся к данному вопросу ниже. Сейчас же мы рассмотрим эффективное представление полиномов по схеме Горнера (Horner structure).
Эффективное представление полиномов по схеме Горнера
Рассмотрим полином (3) более подробно:
\begin{equation*} s(t) = a_0 + a_1 \cdot t + a_2 \cdot t^2 + \ldots + a_{N-1} \cdot t^{N-1}. \end{equation*}
(8)
Вынесем $t$ за скобки, получим:
\begin{equation*} s(t) = a_0 + t \cdot \left(a_1 + a_2 \cdot t + \ldots + a_{N-1} \cdot t^{N-2} \right). \end{equation*}
(9)
Снова вынесем $t$ за скобки, получим:
\begin{equation*} s(t) = a_0 + t \cdot \left(a_1 + t \cdot \left(a_2 + \ldots + a_{N-1} \cdot t^{N-3} \right)\right). \end{equation*}
(10)
Таким образом, вынося $t$ за скобки возможное число раз, получим множество вложенных скобок:
$$ s(t) = a_0 + t \cdot \left( a_1 + t \cdot \left( a_2 + \ldots + t \cdot \left(a_{N-2}+ a_{N-1} \cdot t \right) \right)\right) = $$ $$ = t \cdot \left( t \cdot \left( \ldots t \cdot \left( t \cdot a_{N-1} +a_{N-2} \right) \ldots + a_2 \right) + a_1 \right) + a_0. $$
(11)
Например, для $N = 4$ получим кубический полином:
$$ s(t)= t \cdot \left( t \cdot \left( t \cdot a_3 + a_2 \right) + a_1 \right) + a_0, $$
(12)
который может быть реализован в виде структуры, показанной на рисунке 3.
Рисунок 3. Расчет кубического полинома по схеме Горнера.
Данное представление называется схемой Горнера (Horner structure) и широко используется для эффективного расчета значений полиномов [6]. Ее основное преимущество в том, что для расчета значения полинома степени $N-1$, при известных коэффициентах, требуется $N-1$ умножитель и сумматор.
Использование полиномиальной интерполяции для цифровой передискретизации сигналов
Построение полинома для интерполяции часто проводят при использовании ортогональных полиномов Лагранжа. Однако решение системы линейных уравнений (4) является более общим подходом, который может быть применен не только для Лагранжевой интерполяции, но и для других методов полиномиальной интерполяции (Эрмитова полиномиальная интерполяция, сплайн-интерполяция и других).
Мы уже говорили о том, что операция вычисления обратной матрицы является вычислительно сложной задачей. Кроме того, входной поток поток данных может быть очень длинным, и для обработки мы можем применить кусочную полиномиальную интерполяцию. При этом встает вопрос выбора порядка полинома (количество отсчетов входного сигнала для интерполяции).
Наиболее часто используемым методом является кусочно-полиномиальная кубическая интерполяция ввиду компромисса между точностью и вычислительными затратами. Рассмотрим метод кусочно-полиномиальной кубической интерполяции более подробно.
Полином третьей степени может быть построен по четырем точкам согласно выражению (7). Для любого $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м. рисунок 4a), тогда вместо $x_k$ мы будем использовать значение $ 0 \leq \Delta_k < 1,$ как это показано на рисунке 4b. В этом случае мы можем один раз расчитать обратную матрицу $\mathbf{M^{-1}}$ и использовать ее для любого $x_k$.
Система уравнений для расчета коэффициентов полинома при фиксированных индексах, согласно рисунку 4b, имеет вид:
\begin{equation*} \begin{cases} a_0 + a_1 \cdot (-2) + a_2 \cdot (-2)^2 + a_3 \cdot (-2)^3 = s(n-3),\\ a_0 + a_1 \cdot (-1) + a_2 \cdot (-1)^2 + a_3 \cdot (-1)^3 = s(n-2),\\ a_0 + a_1 \cdot 0 + a_2 \cdot 0^2 + a_3 \cdot 0^3 = s(n-1),\\ a_0 + a_1 \cdot 1 + a_2 \cdot 1^2 + a_3 \cdot 1^3 = s(n).\\ \end{cases} \end{equation*}
(13)
Тогда мы можем представить матрицу $\mathbf{M}$ и получить обратную ей $\mathbf{M^{-1}}$:
\begin{equation*} \begin{array} & \mathbf{M} = \left[ \begin{array}{cccc} 1 & -2 & 4 & -8\\ 1 & -1 & 1 & -1\\ 1 & 0 & 0 & 0\\ 1 & 1 & 1 & 1\\ \end{array} \right], & \mathbf{M^{-1}} = \left[ \begin{array}{cccc} 0 & 0 & 1 & 0\\ \frac{1}{6} & -1 & \frac{1}{2} & \frac{1}{3}\\ 0 & \frac{1}{2} & -1 & \frac{1}{2}\\ -\frac{1}{6} & \frac{1}{2} & -\frac{1}{2}& \frac{1}{6}\\ \end{array} \right]. \end{array} \end{equation*}
(14)
Мы не приводим процедуру расчета обратной матрицы. Вы можете самостоятельно убедиться в том, что $\mathbf{M \cdot M^{-1} = I}$, где $\mathbf{I} - $ единичная матрица.
Коэффициенты полинома, согласно (7), могут быть получены как результат умножения обратной матрицы $\mathbf{M^{-1}}$ на вектор отсчетов входного сигнала $\mathbf{s}$:
\begin{equation*} \begin{array} & \mathbf{a} = \mathbf{M^{-1}} \cdot \mathbf{s} = \left[ \begin{array}{cccc} 0 & 0 & 1 & 0\\ \frac{1}{6} & -1 & \frac{1}{2} & \frac{1}{3}\\ 0 & \frac{1}{2} & -1 & \frac{1}{2}\\ -\frac{1}{6} & \frac{1}{2} & -\frac{1}{2}& \frac{1}{6}\\ \end{array} \right] \cdot \left[ \begin{array}{c} s(n-3) \\ s(n-2) \\ s(n-1) \\ s(n)\\ \end{array} \right]. \end{array} \end{equation*}
(15)
Произведем матричное умножение (15) и получим выражения для коэффициентов полинома:
\begin{equation*} \begin{cases} a_0 = s(n-1),\\ a_1 = \frac{1}{6} \cdot s(n-3) - s(n-2) + \frac{1}{2} \cdot s(n-1) + \frac{1}{3} \cdot s(n),\\ a_2 = \frac{1}{2} \cdot s(n-2) - s(n-1) + \frac{1}{2} \cdot s(n),\\ a_3= -\frac{1}{6} \cdot s(n-3) + \frac{1}{2} \cdot s(n-2) -\frac{1}{2} \cdot s(n-1) +\frac{1}{6} \cdot s(n).\\ \end{cases} \end{equation*}
(16)
Из выражения (16) можно заметить, что коэффициенты полинома зависят от задержанных значений входного сигнала $s(n)$. Мы можем трактовать формулы (16) для расчета коэффициентов как разностные уравнения КИХ-фильтра. Тогда расчет каждого коэффициента будет выполнять свой КИХ-фильтр, и структурная схема фильтра передискретизации может быть представлена, как это показано на рисунке 5.
Каждый КИХ-фильтр рассчитывает один коэффициент полинома, после чего они подаются на схему Горнера для расчета полинома при текущем значении дробной задержки $\Delta_k$. Структура, приведенная на рисунке 5, носит название фильтра Фарроу (Farrow filter) в честь ученого, впервые применившего подобный фильтр для изменения задержки сигнала [2-5].
Рисунок 5. Структура фильтра Фарроу для цифровой передискретизации сигналов на основе
полиномиальной Лагранжевой интерполяции
Важно отметить, что выход фильтра Фарроу $y(k)$ отстает от входного сигнала $s(n)$ на один отсчет в единицах отсчетов входного сигнала (при частоте дискретизации $F_{\textrm{s}}$).
Использование фиксированных значений оси абсцисс позволяет избавиться от необходимости обращения матрицы в реальном времени. Однако фильтр Фарроу также требует значительного количества умножителей на $\frac{1}{6}$, $\frac{1}{3}$ и $\frac{1}{2}$. При этом, в случае реализации в целочисленной арифметике, умножения на $\frac{1}{2}$ можно считать тривиальными, потому что они реализуются как сдвиг на один разряд вправо.
В следующем параграфе мы модифицируем фильтр Фарроу для уменьшения количества умножителей при расчете коэффициентов полинома.
Минимизация умножителей при расчете интерполяционного полинома Лагранжа
При реализации фильтров передискретизации сигналов мы должны стремиться минимизировать количество умножителей при расчете коэффициентов полинома.
Посмотрим внимательно на выражение (16). Коэффицент $a_3$ можно представить в виде:
\begin{equation*} a_3= \frac{1}{6} \cdot \left( s(n) - s(n-3)\right) + \frac{1}{2} \cdot \left( s(n-2) - s(n-1)\right). \end{equation*}
(17)
В выражении (16) также можно заметить, что
\begin{equation*} a_1 + a_3= -\frac{1}{2} \cdot s(n-2) + \frac{1}{2} \cdot s(n), \end{equation*}
(18)
откуда можно рассчитать коэффициент $a_1$ как:
\begin{equation*} a_1 = \frac{1}{2} \cdot \left( s(n) - s(n-2) \right)-a_3. \end{equation*}
(19)
Для расчета коэффициента $a_2$ можно использовать следующее следствие из (16):
\begin{equation*} a_1 + a_2 + a_3 = s(n) - s(n-1), \end{equation*}
(20)
тогда:
\begin{equation*} a_2 = s(n) - s(n-1) - a_1 - a_3. \end{equation*}
(21)
Таким образом, мы можем переписать (16) в оптимизированном виде всего с тремя умножителями (причем два из них равны $\frac{1}{2}$):
\begin{equation*} \begin{cases} a_0 = s(n-1),\\ a_3= \frac{1}{6} \cdot \left( s(n) - s(n-3)\right) + \frac{1}{2} \cdot \left( s(n-2) - s(n-1)\right), \\ a_1 = \frac{1}{2} \cdot \left( s(n) - s(n-2) \right)-a_3,\\ a_2 = s(n) - s(n-1) - a_1 - a_3.\\ \end{cases} \end{equation*}
(22)
Структурная схема оптимизированного фильтра цифровой передискретизации сигналов, соответствующая (22), приведена на рисунке 6.
Рисунок 6. Структурная схема оптимизированного фильтра цифровой передискретизации
Таким образом, мы получили алгоритм, который рассчитывает коэффициенты интерполяционного полинома Лагранжа с использованием всего одного умножителя на $\frac{1}{6}$ и двух тривиальных умножителей на $\frac{1}{2}$. В следующих разделах мы рассмотрим примеры использования фильтра цифровой передискретизации на основе полиномиальной Лагранжевой интерполяции.
Выводы
В данном разделе мы рассмотрели структурные схемы фильтров цифровой передискретизации на основе полиномиальной Лагранжевой интерполяции.
Была рассмотрена схема Горнера для расчета значений полинома. Количество умножителей и сумматоров схемы Горнера равно порядку полинома.
Получена структура фильтра Фарроу, и произведена его оптимизация. В результате коэффициенты интерполяционного полинома Лагранжа могут быть рассчитаны при использовании всего трех умножителей. Требуется одно умножение на $\frac{1}{6}$ и два тривиальных умножения на $\frac{1}{2}$.
В следующем разделе мы рассмотрим характеристики фильтра цифровой передискретизации и примеры его использования для компенсации дробной задержки, интерполяции сигналов и дробного изменения частоты дискретизации сигнала.

Ваши комментарии, вопросы и пожелания вы можете оставить на форуме или в гостевой книге.
Список литературы
[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 Floyd M. Interpolation in Digital Modems-Part I: Fundamentals: IEEE Transactions on Communications, Vol. 41, No. 3, March 1993, P. 501-507.

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

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

[6] Макконелл Дж. Основы современных алгоритмов. Москва, Техносфера, 2004.


Oбнаружили ошибку в тексте? Выделите ее мышкой и нажмите
Описание (необязательно)
Закрыть Х