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

Содержание
Введение
В предыдущем разделе мы рассмотрели фильтр Фарроу для цифровой передискретизации сигналов. Получена структурная схема эффективного фильтра цифровой передискретизации (ресэмплинга сигналов) на основе кусочно-кубической полиномиальной интерполяции (рисунок 1).

Рисунок 1. Структурная схема оптимизированного фильтра Фарроу
на основе кусочно-кубической полиномиальной интерполяции.
Расчет коэффициентов интерполяционного полинома требует одно умножение на $\frac{1}{6}$ и два тривиальных умножения на $\frac{1}{2}$.
В данном разделе мы рассмотрим примеры использования фильтра передискретизации для компенсации дробной задержки, цифровой интерполяции и дробного изменения частоты дискретизации в $\frac{P}{Q}$ раз.
Пересчет индексов отсчетов входного сигнала при кусочно-кубической интерполяции
Высокая эффективность цифрового передискретизатора достигается за счет использования кусочно-кубической интерполяции и пересчета индексов $n$ отсчетов входного сигнала $s(n)$, $n =0,1,2, \ldots$ в интервал $[-2, \ldots 1]$, как это показано на рисунке 2.
Рисунок 2. Пересчет индексов отсчетов входного сигнала в интервал $[-2, \ldots 1].$
Пусть частота дискретизации входного сигнала $s(n)$, $n =0,1,2, \ldots$ равна $F_\textrm{s}$. Тогда сигнал $y(k)$, $k =0,1,2, \ldots$ на выходе фильтра передискретизатора будет иметь частоту дискретизации $F_\textrm{y} = \frac{P}{Q} \cdot F_\textrm{s}$. Кроме того первый отсчет выходного сигнала $y(k)$ задержан относительно входного $s(n)$ на величину $0 \leq x_0 <1$.
Везде далее в данном разделе, так же как и в предыдущем разделе, переменная $n$ индексирует отсчеты входного сигнала $s(n)$, а переменная $k$ индексирует отсчеты выходного сигнала $y(k)$.
Поскольку частоты дискретизации входного и выходного сигналов различны, то одни и те же индексы входного и выходного отсчетов соответствуют различным моментам времени. Например, $5-$ый отсчет выходного сигнала $y(5)$ после передискретизации не совпадает с $5-$м отсчетом входного сигнала $s(5)$. Поэтому мы должны указывать, по какой шкале мы индексируем отсчет. Мы будем использовать три шкалы: шкала $n$ для моментов времени относительно отсчетов входного сигнала $s(n)$, шкала $k$ для моментов времени относительно выходного сигнала и шкала $t$ в интервале $[-2 \ldots 1]$ для расчета коэффициентов кубического полинома.
Для $k-$го отсчета выходного сигнала $y(k)$ по шкале $k$ необходимо рассчитать соответствующие индексы $n$ входного сигнала $s(n)$ по шкале $n$ (рисунок 2a), а затем перейти к шкале $t$ в интервале $[-2 \ldots 1]$ для расчета коэффициентов кубического полинома (рисунок 2b). При этом момент взятия $k-$го выходного отсчета по шкале $n$ должен быть пересчитан в значение $0 \leq \Delta_k < 1$ по шкале $t$.
Момент $x_k$ взятия $k-$го выходного отсчета после передискретизации по шкале $n$ входного сигнала равен:
$$ x_k = k \cdot \frac{Q}{P} - x_0. $$
(1)
Для расчета кубического полинома нам требуется четыре отсчета входного сигнала $s(n)$, $s(n-1)$,$s(n-2)$ и $s(n-3)$, как это показано на рисунке 2a. Причем два отсчета $s(n)$ и $s(n-1)$ должны располагаться правее $x_k$, а два других отсчета $s(n-2)$ и $s(n-3)$ - левее.
Тогда индекс $n$ отсчета $s(n)$, соответствующего моменту времени $t = 1$ по шкале $t$ от -2 до 1, равен:
$$ n = \lfloor x_k \rfloor + 2, $$
(2)
где оператор $\lfloor x_k \rfloor$ означает округление к ближайшему меньшему (floor).
Значение $\Delta_k$ по шкале $t$ равно:
$$ \Delta_k = \lfloor x_k \rfloor + 1 - x_k. $$
(3)
На рисунке 3 показан пример пересчета индексов входного и выходного сигналов для $5-$го выходного отсчета $y(5)$ при $P = 4$, $Q = 3$, $x_0 = 0.2$.
Рисунок 3. Пример пересчета индексов при $P = 4$, $Q = 3$, $x_0 = 0.2$
Момент $x_5$ по шкале $n$, согласно (1), равен:
$$ x_5 = 5 \cdot \frac{3}{4} - 0.2 = 3.55. $$
(4)
Тогда индекс $n$, соответственно (2), равен:
$$ n = \lfloor 3.55 \rfloor + 2 = 5, $$
(5)
и параметр $\Delta_5$ по шкале $t$ равен:
$$ \Delta_5 = \lfloor 3.55 \rfloor + 1 - 3.55 = 0.45. $$
(6)
Таким образом, необходимо рассчитать коэффициенты кубического полинома $P(t) = a_0 + a_1 \cdot t +a_2 \cdot t^2 + a_3 \cdot t^3 $ согласно формулам:
\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*}
(7)
после чего значение $y(k)$ можно получить как значение интерполяционного полинома $P(t)$ для $t=-\Delta_k$, т.е. $y(k) = P \left(-\Delta_k \right)$.
Использование фильтра Фарроу на основе полиномиальной интерполяции для компенсации дробной задержки сигнала
Пусть входной сигнал $s(n)$ содержит $N=8$ отсчетов, как это показано на рисунке 4.

Рисунок 4. Входной сигнал фильтра Фарроу.
Произведем компенсацию дробной задержки данного сигнала на величину $x_0 = 0.25$ периода дискретизации. При изменении дробной задержки частота дискретизации остается неизменной, параметры $P$ и $Q$ равны $P = Q = 1$, и количество отсчетов на выходе фильтра Фарроу равно количеству отсчетов на входе $N=8$.
В таблице 1 приведен процесс компенсации дробной задержки сигнала $s(n),$ $n = 0 \ldots 7,$ показанного на рисунке 4.
Таблица 1. Компенсация дробной задержки сигнала $s(n)$.
Значения $x_k,$ $n$ и $\Delta_k$ рассчитываются согласно (1), (2) и (3) соотвественно. Далее приведены значения $s(n-3) \ldots s(n)$ для каждого значения $k$. Красным отмечены нулевые значения $s(n-3) \ldots s(n)$ для $n<3$ и $n>N-1$, выходящие за пределы индексации сигнала $s(n),$ $n = 0 \ldots N-1$. В следующих столбцах таблицы 1 приведены значения коэффициентов кубического полинома $a_0 \ldots a_3$, рассчитанные согласно (7) для текущего момента $k$. Наконец, в последнем столбце значение сигнала $y(k) = P(-\Delta_k)$ на выходе фильтра Фарроу.
На рисунке 5 показан исходный сигнал $s(n),$ $n = 0 \ldots 7 $ (черный) и сигнал $y(k),$ $n = 0 \ldots 7 $ после компенсации дробной задержки $x_0$ (красный).

Рисунок 5. Входной сигнал $s(n)$ фильтра Фарроу и сигнал после компенсации дробной задержки $y(k)$.
Из таблицы 1 следует, что при компенсации дробной задержки параметры $\Delta_k = x_0$ для всех $k = 0 \ldots N-1$. В этом случае фильтр Фарроу можно трактовать как КИХ-фильтр 3-го порядка с изменяющейся групповой задержкой $\tau(\omega)$ путем задания требуемого значения $x_0$.
На рисунке 6 показано семейство импульсных характеристик фильтра Фарроу для различного значения дробной задержки $x_0$.

Рисунок 6. Семейство импульсных характеристик фильтра Фарроу для различного значения дробной задержки $x_0$
Амплитудно-частотная характеристика (АЧХ) $\left| H \left(\textrm{e}^{j\cdot \omega} \right) \right|^2$ и групповая задержка $\tau(\omega)$ полученных фильтров показаны на рисунке 7.

Рисунок 7. АЧХ и групповая задержка фильтров Фарроу для различного значения дробной задержки $x_0$
Из графиков АЧХ и групповой задержки можно сделать вывод, что коррекция задержки возможна только в полосе до $0.4\pi$ рад/с. Для более широкой полосы коррекции фильтр Фарроу на основе полиномиальной кусочно-кубической интерполяции не подходит ввиду недопустимо высокого искажения характеристик фильтра. Если требуется использовать компенсацию дробной задержки в более широкой полосе, то необходимо применить более сложные фильтры, использующие полиномиальную интерполяцию более высокого порядка, или альтернативные КИХ-фильтры, или БИХ-фильтры коррекции групповой задержки сигналов.
Программа, реализующая расчет таблицы 1 и данных для построения рисунка 5, написана с использованием библиотеки DSPL. Исходный код программы можно получить в разделе примеров библиотеки DSPL.
Пользователи Matlab и GNU Octave могут произвести расчет таблицы 1 выполнив скрипт resample_lagrange_fd.m.
						
clear all;
close all;
clc;

p = 1;
q = 1;
x0 = 0.25;



s0 = [1 2 2 1 -0.5 -1 -2 -0.5];

tin = 0:length(s0)-1;

y = zeros(1, floor(length(s0)*q/p));

s = [0, 0, s0, 0, 0];


fprintf('---------------------------------------------------------');
fprintf('-------------------------------------------------------\n');
fprintf('k     x_k    n        d    s(n-3)    s(n-2)   s(n-1)    s(n)     ');
fprintf('a_0       a_1       a_2       a_3       y(k)\n');
fprintf('---------------------------------------------------------');
fprintf('-------------------------------------------------------\n');
for k = 0:length(y)-1
  x =  k*q/p - x0       ;
  n =  floor(x) + 2 + 3 ;
  d =  floor(x) + 1 - x ;

  tout(k+1) = x;

  a0 = s(n-1);
  a3 = 1/6 * (s(n) - s(n-3)) + 0.5*(s(n-2) - s(n-1));
  a1 = 0.5 * (s(n) - s(n-2)) - a3;
  a2 = s(n) - s(n-1) - a3 - a1;
    
  y(k+1) = a0 - a1 * d + a2*d^2 - a3*d^3;  
    
  fprintf('%d  %7.2f   %d   %7.2f', k, x, n, d);
  fprintf('%8.1f %8.1f %8.1f %8.1f ', s(n-3), s(n-2), s(n-1), s(n));
  fprintf('%9.3f %9.3f %9.3f %9.3f ', a0, a1,a2,a3);
  fprintf('%9.3f\n', y(k+1));

end

fprintf('---------------------------------------------------------');
fprintf('-------------------------------------------------------\n');


figure; stem(tin, s0); 
hold on;  
stem(tout, y, 'r'); 
axis([-0.5, 7.5, -2.5, 2.5]);
grid on;


Расчет семейства импульсных характеристик, АЧХ и групповой задержки фильтров Фарроу для различного значения дробной задержки $x_0$, показанная на рисунках 6 и 7, также производилась при использовании библиотеки DSPL. Исходный код программы расчета АЧХ и групповой задержки фильтров Фарроу также можно получить в соответствующем разделе примеров библиотеки DSPL.
Пользователи Matlab и GNU Octave могут произвести расчет АЧХ и групповой задержки фильтров Фарроу для различного значения дробной задержки $x_0$ выполнив скрипт resample_lagrange_fd_filter.m.
						
clear all;
close all;
clc;

addpath('functions');
 
% Check Matlab or GNU Octave
isOctave = exist('OCTAVE_VERSION', 'builtin');
if(isOctave)
  pkg load signal;
end

x0 = 0:0.1:0.9;

s = [0, 1, 0, 0];

figure(1); hold on;
figure(2); hold on;
figure(3); hold on;

w  = 0:0.001:pi;

for k = 1:length(x0)
  h = resample_lagrange(s, 1, 1, x0(k));
  figure(1); subplot(2,5,k); stem(0:3, h);
  axis([0, 3, -0.2, 1.2]); grid on;
  
  H = freqz(h, 1, w);
  tau = grpdelay(h,1,w);
    
  figure(2); plot(w/pi, 20*log10(abs(H))); 
  axis([0, 1, -20, 2]); grid on;
  xlabel('w/pi');
  ylabel('|H(e^{jw})|^2, dB');
  
  figure(3); plot(w/pi, tau); 
  axis([0, 1, 0, 3]); grid on;
  xlabel('w/pi');
  ylabel('\tau(w)');
end
Для выполнения скрипта resample_lagrange_fd_filter.m и других скриптов, приведенных ниже, требуется функция resample_lagrange.m.
						
function y = resample_lagrange(s, p, q, x0)
% y = resample_lagrange(s, p, q, x0)
% Digital resampling by polynomial Lagrange interpolation.
% Function changes input signal s samplerate to p/q times and adds fractional
% delay.
%
% Input parameters
%  s   - input signal vector [N x 1];
%  p   - p paramter of samplarate conversion
%  q   - q paramter of samplarate conversion
%  x0  - fractional delay
%
% Ouptut parameters
%  y   - Resampled signal
% 
% Author: Sergey Bakhurin (dsplib.org)
%


if(p>1 && q==1)
  y = zeros(1, floor((length(s)-1)*p/q)+1);
else
  y = zeros(1, floor(length(s)*p/q));
end

t = zeros(size(y));

s = [0, 0, s, 0, 0];

for k = 0:length(y)-1
  x = k*q/p - x0;
  
  t(k+1) = x;
  
  n = floor(x) + 5;
  d = floor(x)+1-x;

  
  a0 = s(n-1);
  a3 = 1/6 * (s(n) - s(n-3)) + 0.5*(s(n-2) - s(n-1));
  a1 = 0.5 * (s(n) - s(n-2)) - a3;
  a2 = s(n) - s(n-1) - a3 - a1;
  
  y(k+1) = a0 - a1 * d + a2*d^2 - a3*d^3;  
  
end

Использование фильтра Фарроу в качестве цифрового интерполятора сигналов
Рассмотрим пример использования фильтра цифровой передискретизации для цифровой интерполяции сигнала.
Пусть входной сигнал $s(n),$ $n = 0 \ldots N-1,$ содержит $N=8$ отсчетов сигнала, как это показано на рисунке 8a. Мы уже использовали данный сигнал в качестве входного сигнала цифрового фильтра коррекции дробной задержки.
Необходимо произвести цифровую интерполяцию сигнала и увеличить частоту дискретизации сигнала $y(k),$ $k = 0\ldots K-1$ на выходе фильтра Фарроу в $P=10$ раз.
При цифровой интерполяции количество отсчетов сигнала на выходе будет равно (смотри рисунок 8b):
\begin{equation*} K = (N-1) \cdot P +1 = (8-1)\cdot 10 +1 = 71. \end{equation*}
(8)
Пересчет значений $x_k,$ $n$ и $\Delta_k$ также производится согласно (1), (2) и (3) соотвественно.
Результат цифровой интерполяции при использовании фильтра Фарроу показан на рисунке 8b. Черным цветом отмечены узлы интерполяции, соответствующие входному сигналу.
Рисунок 8. Входной сигнал и результат цифровой интерполяции при использовании фильтра Фарроу.
Исходный код программы, реализующей расчет данных для построения рисунка 8, также приведен в соответствующем разделе справки библиотеки DSPL.
Пользователи Matlab и GNU Octave могут произвести цифровую интерполяцию с использование фильтра Фарроу выполнив скрипт resample_lagrange_interp.m.
						
clear all;
close all;
clc;

addpath('functions');

p = 10;
q = 1;
x0 = 0;

s = [1 2 2 1 -0.5 -1 -2 -0.5];

y = resample_lagrange(s, p, q, x0);

figure(1); 
stem(0:length(y)-1, y,'r'); 
hold on; 
stem((0:length(s)-1)*p, s, 'k');

Фильтр цифровой интерполяции можно трактовать как фильтр нижних частот. Импульсная характеристика $h(n)$ и амплитудно-частотная характеристика $\left|H \left( \operatorname{e}^{j\cdot \omega}\right) \right|^2$ фильтра Фарроу цифровой интерполяции сигнала показаны на рисунке 9 для $P = 10$.

Рисунок 9. Импульсная характеристика и АЧХ фильтра Фарроу.
Пунктирной линией показана импульсная характеристика при $P \rightarrow\infty$.
На рисунке 9 можно заметить, что импульсная характеристика не имеет непрерывной производной в узлах интерполяции (имеет точки излома) ввиду того, что построение интеполяционного полинома Лагранжа ведется только по отсчетам сигнала без ограничений на непрерывность производной. В результате чего уровень подавления в полосе заграждения АЧХ фильтра составляет всего 28 дБ.
Исходный код программы расчета импульсной характеристики и АЧХ фильтра Фарроу цифровой интерполяции сигнала (рисунка 9) также приведен в справке библиотеки DSPL.
Пользователи Matlab и GNU Octave могут произвести расчет характеристик цифрового фильтра интерполятора выполнив скрипт resample_lagrange_interp_filter.m.
						

clear all;
close all;
clc;

addpath('functions');

p = 10;
q = 1;
x0 = 0;

s = [0 0 1 0 0];

h = resample_lagrange(s, p, q, x0);


w = 0:0.01:pi;
H = freqz(h,1,w);

figure(1); subplot(1,2,1); 
plot(0:length(h)-1, h, 'r'); hold on; 
stem(0:length(h)-1, h); 

 
xlabel('n'), ylabel('h(n)');
grid on;

figure(1); subplot(1,2,2);  
plot(w/pi, 20*log10(abs(H))); 
axis([0, 1, -60, 30]);
xlabel('w/pi');
ylabel('|H(e^{jw})|^2, dB'); grid on;
Использование фильтра Фарроу для дробного изменения частоты дискретизации сигнала
Пусть синусоидальный входной сигнал
\begin{equation*} s(n) = \sin \left( 2\pi \cdot n \cdot \frac{f_0}{F_{\textrm{s}}} \right), ~ n =0 \ldots N-1 \end{equation*}
(9)
с частотой $f_0 = 6$ кГц содержит $N = 54$ отсчета взятых с частотой дискретизации $F_{\textrm{s}}=26.4 $ кГц.
Отношение частоты дискретизации $F_{\textrm{s}}$ синусоидального сигнала $s(n)$ к частоте данного сигнала $f_0$ не является целым числом $\frac{F_{\textrm{s}}}{f_0} = \frac{26.4}{6}$. В результате на одном периоде синусоидального сигнала не укладывается целое число отсчетов, и исходный сигнал показан на верхнем графике рисунка 10.
Произведем изменение частоты дискретизации входного сигнала $s(n)$ в $\frac{P}{Q}$ раз, где $P = 20$, $Q = 11$ и получим передискретизированный сигнал $y(k)$, отсчеты которого получены с частотой дискретизации $F_{\textrm{y}} = \frac{P}{Q} \cdot F_{\textrm{s}} = 48$ кГц. В результате на одном периоде сигнала $y(k)$ будет укладываться ровно 8 отсчетов. Дробное изменение частоты дискретизации будем производить при использовании фильтра Фарроу на основе полиномиальной интерполяции. Пересчет значений $x_k,$ $n$ и $\Delta_k$ производится согласно (1), (2) и (3) соотвественно.
На нижнем графике рисунка 10 показан сигнал $y(k)$ после цифровой передискретизации при использовании фильтра Фарроу.
Рисунок 10. Дробное изменение частоты дискретизации сигнала.
Непрерывными линиями на рисунке 10 показан исходный непрерывный синусоидальный сигнал. Из рисунка 10 следует, что после передискретизации на одном периоде исходного сигнала укладывается ровно восемь отсчетов сигнала $y(k)$.
Исходный код программы расчета данных примера цифровой передискретизации вы также можете найти в соответствующем разделе справки DSPL.
Пользователи Matlab и GNU Octave могут произвести цифровую дискретизацию выполнив скрипт resample_lagrange_pq.m.
						
clear all;
close all;
clc;

addpath('functions');

Fs = 26.4;
f0 = 6;
N = 54;


p = 20;
q = 11;
x0 = 0;

t = (0:N-1)/Fs;
s  = sin(2*pi*t*f0);
t0 = 0:0.01:2;
s0 = sin(2*pi*t0*f0);

y = resample_lagrange(s, p, q, x0);
tr = (0:length(y)-1)/Fs/p*q;

figure(1);
subplot(2,1,1); hold on;
stem(t, s);
%plot(t0 ,s0, 'r');
axis([0,2,-1.5,1.5]); grid on;
 
 
subplot(2,1,2);  hold on;
stem(tr, y);
%plot(t0 ,s0, 'r');
axis([0,2,-1.5,1.5]); 
grid on;
Выводы
В данном разделе мы рассмотрели вопрос пересчета индексов отсчетов входного сигнала при кусочно-кубической интерполяции. Были получены выражения для пересчета моментов дискретизации выходного сигнала по шкале $n$ входного сигнала, а также расчет значения $\Delta_k $ по шкале $t$.
Полученные выражения были использованы в примерах использования фильтра Фарроу для компенсации дробной задержки, цифровой интерполяции и дробного изменения частоты дискретизации сигнала.
Были приведены семейства импульсных характеристик, АЧХ и группового времени запаздывания фильтров Фарроу для различного значения дробной задержки.
Также приведены импульсная характеристика и АЧХ фильтра Фарроу цифровой интерполяции сигнала.
Даны ссылки на исходные коды программ, реализующих расчет всех показанных примеров с использованием библиотеки DSPL и пакетов Matlab и GNU Octave.

Ваши комментарии, вопросы и пожелания вы можете оставить на форуме или в гостевой книге.
Список литературы
[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бнаружили ошибку в тексте? Выделите ее мышкой и нажмите
Описание (необязательно)
Закрыть Х