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

Содержание

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

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

Страница проекта на SourceForge libdspl-2.0 on SourceForge

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

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

Структурная схема оптимизированного фильтра Фарроу
на основе кубических сплайнов Эрмита
Рисунок 1. Структурная схема оптимизированного фильтра Фарроу на основе кубических сплайнов Эрмита

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

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

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

Пересчет индексов отсчетов входного сигнала при кусочно-кубической сплайн-интерполяции

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

Высокая эффективность цифрового передискретизатора достигается за счет использования кусочно-кубической сплайн-интерполяции и пересчета индексов n отсчетов входного сигнала s(n), n =0,1,2, \ldots в интервал [-2, \ldots 1], как это показано на рисунке 2.

Пересчет индексов отсчетов входного сигнала
в интервал 
Рисунок 2. Пересчет индексов отсчетов входного сигнала в интервал t\in[-2, \,\, 1]

Пусть частота дискретизации входного сигнала s(n), n =0,1,2, \ldots равна F_\textrm{s}. Тогда сигнал y(k), k =0,1,2, \ldots на выходе фильтра передискретизатора будет иметь частоту дискретизации F_\textrm{y} = \frac{P F_\textrm{s}}{Q}. Кроме того, первый отсчет выходного сигнала 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 (рисунок 2а), а затем перейти к шкале t в интервале [-2 \ldots 1] для расчета коэффициентов кубического сплайна (рисунок 2б). При этом, момент взятия k-го выходного отсчета по шкале n должен быть пересчитан в значение 0 \leq \Delta_k < 1 по шкале t.

Момент x_k взятия k-го выходного отсчета после передискретизации по шкале n входного сигнала равен:

equation 1
(1)
Для расчета кубического сплайна нам требуется четыре отсчета входного сигнала s(n)s(n-1), s(n-2) и s(n-3), как это показано на рисунке 2а. Причем два отсчета s(n) и s(n-1) должны располагаться правее x_k, а два других отсчета s(n-2) и s(n-3) — левее.

Тогда индекс n отсчета s(n), соответствующего моменту времени t = 1 по шкале t от -2 до 1, равен:

equation 2
(2)
где оператор \lfloor x_k \rfloor означает округление к ближайшему меньшему (оператор floor).

Значение \Delta_k по шкале t равно:

equation 3
(3)
На рисунке 3 показан пример пересчета индексов входного и выходного сигналов для 5-го выходного отсчета y(5) при P = 4Q  = 3x_0 = 0.2.

Пример пересчета индексов при
, , 
Рисунок 3. Пример пересчета индексов при P = 4Q  = 3x_0 = 0.2

Момент x_5 по шкале n, согласно (1), равен:

equation 4
(4)
Тогда индекс n, соответственно (2), равен:
equation 5
(5)
и параметр \Delta_5 по шкале t равен:
equation 6
(6)
Таким образом, необходимо рассчитать коэффициенты кубического сплайна P(t) = a_0 + a_1  t +a_2  t^2 + a_3  t^3 согласно формулам:
equation 7
(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_S(k) — значение сигнала y_S(k) = P(-\Delta_k) на выходе фильтра Фарроу с использованием сплайн-интерполяции. В последнем столбце y_L(k) для сравнения приведены значения на выходе фильтра Фарроу с использованием полиномиальной интерполяции Лагранжа.

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

На рисунке 5 показан исходный сигнал s(n), n = 0 \ldots 7 и сигнал y_S(k), n = 0 \ldots 7 после компенсации дробной задержки x_0.

Входной сигнал  фильтра Фарроу и
сигнал после компенсации дробной задержки
Рисунок 5. Входной сигнал s(n) фильтра Фарроу и сигнал после компенсации дробной задержки y_S(k)

Из таблицы 1 следует, что при компенсации дробной задержки параметры \Delta_k = x_0 для всех k = 0 \ldots N-1. В этом случае фильтр Фарроу можно трактовать как КИХ-фильтр 3-го порядка с изменяющимися АЧХ \left| H \left(\textrm{e}^{j \omega} \right) \right|^2 и групповой задержкой \tau(\omega) путем задания требуемого значения x_0.

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

АЧХ и групповая задержка фильтров Фарроу для
различного значения дробной задержки 
при использовании сплайн-интерполяции (сплошные)
и интерполяционного полинома Лагранжа (пунктирные)
Рисунок 6. АЧХ и групповая задержка фильтров Фарроу для различного значения дробной задержки x_0 при использовании сплайн-интерполяции (сплошные) и интерполяционного полинома Лагранжа (пунктирные)

Из графиков АЧХ и групповой задержки можно сделать вывод, что фильтр Фарроу при использовании сплайн-интерполяции обладает меньшей неравномерностью АЧХ, но чуть большей неравномерностью групповой задержки по сравнения с фильтром Фарроу на основе полиномиальной интерполяции Лагранжа. При этом коррекция задержки возможна только в нормированной полосе до 0.4\pi рад/с (до 0.4  F_{\textrm{s}} Гц).

Расчет данных для построения таблицы 1 производился с использованием библиотеки DSPL-2.0. Исходный код программы resampling_spline_fd.c расчета данных таблицы 1 приведен в приложении.

Расчет АЧХ и групповой задержки фильтра Фарроу на основе сплайн-интерполяции производился программой resampling_spline_filter_fd.c .

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

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

Пусть входной сигнал s(n), n = 0 \ldots N-1, содержит N=8 отсчетов сигнала, как это показано на рисунке 7а. Мы уже использовали данный сигнал в качестве входного сигнала цифрового фильтра коррекции дробной задержки.

Необходимо произвести цифровую интерполяцию сигнала и увеличить частоту дискретизации сигнала y(k), k = 0\ldots K-1 на выходе фильтра Фарроу в P=10 раз при использовании сплайн-интерполяции.

При цифровой интерполяции количество отсчетов сигнала на выходе будет равно (смотри рисунок 7б):

equation 8
(8)
Пересчет значений x_k, n и \Delta_k также производится

согласно (1), (2) и (3) соответственно.

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

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

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

Импульсная характеристика и АЧХ фильтра
Фарроу на основе сплайн-интерполяции
и полиномов Лагранжа
Рисунок 8. Импульсная характеристика и АЧХ фильтра Фарроу на основе сплайн-интерполяции и полиномов Лагранжа

На рисунке 8 можно заметить, что при использовании полиномов Лагранжа, импульсная характеристика не имеет непрерывной производной в узлах интерполяции (имеет точки излома). В результате, уровень подавления в полосе заграждения АЧХ фильтра составляет всего 29 дБ. С другой стороны, использование сплайнов Эрмита обеспечивает непрерывную производную импульсной характеристики фильтра Фарроу в узлах интерполяции. Уровень подавления в полосе заграждения АЧХ фильтра Фарроу на основе сплайн интерполяции на 12 дБ выше фильтра на основе полиномиальной интерполяции Лагранжа.

Расчет данных для построения рисунков 7 и 8 производился производился программой resampling_spline_interp.c , исходный код которой можно найти в приложении.

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

Пусть синусоидальный входной сигнал

equation 9
(9)
с частотой f_0 = 6 кГц содержит N = 54 отсчета, взятых с частотой дискретизации F_{\textrm{s}}=26.4 кГц.

Мы уже использовали данный сигнал для примера дробной передискретизации на основе полиномиальной интерполяции Лагранжа. Отношение частоты дискретизации F_{\textrm{s}} синусоидального сигнала s(n) к частоте данного сигнала \frac{F_{\textrm{s}}}{f_0} = \frac{26.4}{6} не является целым числом. В результате на одном периоде синусоидального сигнала не укладывается целое число отсчетов. Исходный сигнал показан на верхнем графике рисунка 9.

Произведем изменение частоты дискретизации входного сигнала s(n) в \frac{P}{Q} раз, где P = 20, Q = 11, и получим передискретизированные сигналы y_L(k) и y_S(k), на основе полиномиальной и сплайн-интерполяции соотвественно. Отсчеты сигналов y_L(k) и y_S(k) получены с частотой дискретизации F_{\textrm{y}} = \frac{P}{Q}  F_{\textrm{s}} = 48 кГц. В результате на одном периоде сигнала y(k) будет укладываться ровно 8 отсчетов. Пересчет значений x_k, n и \Delta_k производится

согласно (1), (2) и (3) соответственно.

На втором графике рисунка 9 показан сигнал y_L(k) после цифровой передискретизации при использовании фильтра Фарроу на основе полиномиальной интерполяции Лагранжа. Также на третьем графике приведен сигнал y_S(k) при использовании фильтра Фарроу на основе сплайн-интерполяции. На последнем графике приведены сигналы ошибок передискретизации e_L(k) и e_S(k) на выходе фильтров Фарроу при использовании полиномиальной интерполяции Лагранжа и сплайн-интерполяции соотвественно.

Также на нижних графиках приведены оценки нормированного среднего квадрата ошибок передискретизации \textrm{NMSE} \big[ e_L(k)\big] и \textrm{NMSE} \big[ e_S(k)\big].

Дробное изменение частоты дискретизации сигнала
Рисунок 9. Дробное изменение частоты дискретизации сигнала

Из рисунка 9 следует, что после передискретизации на одном периоде исходного сигнала укладывается ровно восемь отсчетов сигнала как при использовании интерполяции Лагранжа, так и при сплайн-интерполяции. При этом нормированная средне-квадратическая ошибка \textrm{NMSE} \big[ e_S(k)\big] на 0.001 (0.1\%) выше при использовании сплайн-интерполяции.

Расчет данных для построения рисунка 9 производился программой resampling_spline_pq.c , исходный код которой также приведен в приложении.

Выводы

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

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

В данном разделе мы рассмотрели вопрос пересчета индексов отсчетов входного сигнала при кусочно-кубической сплайн-интерполяции. пыли приведены выражения для пересчета моментов дискретизации выходного сигнала по шкале n входного сигнала, а также расчет значения \Delta_k по шкале t.

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

пыли приведены сравнения семейств АЧХ и группового времени запаздывания фильтров Фарроу при использовании полиномиальной интерполяции Лагранжа и сплайн-интерполяции для различного значения дробной задержки.

Лакже приведены сравнения АЧХ фильтра Фарроу цифровой интерполяции сигнала при использовании полиномиальной интерполяции Лагранжа и сплайн-интерполяции. показано, что фильтр передискретизации на основе сплайн-интерполяции обеспечивает непрерывную производную интерполированного сигнала, за счет чего достигается на 12 дБ лучшее подавление в полосе заграждения фильтра.

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

Даны ссылки на исходные коды программ, реализующих расчет всех показанных примеров.

приложения

ƒанные для построения рисунков данного раздела были просчитаны при использовании библиотеки DSPL-2.0.

Ќиже приведЄн исходный код программы на языке —и, расчета данных для реализации фильтра дробной задержки (рисунок 5 и таблица 1)


#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "dspl.h"

int main(int argc, char* argv[])
{
  double s0[8] = {1.0, 2.0, 2.0, 1.0, -0.5, -1.0, -2.0, -0.5}; /* исходный сигнал   */
  double *y = NULL;
  int    i, ny;
  
  void* hdspl; /* DSPL handle        */

  /* Загружаем libdspl.dll */
  hdspl = dspl_load();
  if(!hdspl)
  {
    printf("Cannot to load libdspl!\n");
    return 0;
  }
  
  farrow_spline(s0, 8, 1.0, 1.0, 0.25, &y, &ny);
  for(i = 0; i < ny; i++)
    printf("y[%d] = %+-6.3f\n", i, y[i]);

  /* Выгружаем libdspl.dll из памяти */
  dspl_free(hdspl);

  return 0;
}

јналогичный скрипт на языке Python:


# -*- coding: utf-8 -*-
"""
Created on Sun Aug 27 23:40:00 2017

@author: sergey
"""

import numpy as np

p = 1
q = 1
x0 = 0.25



s0 = np.array([1., 2., 2., 1., -0.5, -1., -2., -0.5])

tin = np.linspace(0.0, len(s0), len(s0), endpoint=False)

y = np.zeros(int(np.floor( float(len(s0)*q)/float(p))))

s = np.concatenate((np.array([0., 0.]), s0, np.array([0., 0.])))


print('---------------------------------------------------------', end='');
print('-------------------------------------------------------\n', end='');
print('k     x_k    n        d    s(n-3)    s(n-2)   s(n-1)    s(n)     ', end='');
print('a_0       a_1       a_2       a_3       y(k)\n', end='');
print('---------------------------------------------------------', end='');
print('-------------------------------------------------------\n', end='');

tout = np.zeros(int(np.floor( float(len(s0)*q)/float(p))))
for k  in range(len(y)):
    x =  k*q/p - x0       
    n =  int(x) + 1 + 3 
    d =  int(x) + 1 - x 

    tout[k] = x;
    
    a0 = s[n-1];
    a1 = 0.5*(s[n] - s[n-2]);
    a3 = 2.0*(s[n-2] - s[n-1]) + a1 + 0.5*(s[n-1] - s[n-3]);
    a2 = s[n-2] - s[n-1] +a3 + a1;
    
    y[k] = a0 - a1 * d + a2*d**2 - a3*d**3
    
    print('%(d0)d  %(f0)7.2f   %(d1)2d   %(f1)7.2f'  % 
        {'d0': k, 'f0': x, 'd1': n,  'f1': d}, end='')
    print('%(s0)8.1f %(s1)8.1f %(s2)8.1f %(s3)8.1f ' % 
        {'s0': s[n-3],'s1': s[n-2], 's2': s[n-1], 's3':s[n]}, end='')
  
    print('%(a0)9.3f %(a1)9.3f %(a2)9.3f %(a3)9.3f '% 
        {'a0': a0, 'a1': a1, 'a2': a2, 'a3': a3}, end='')
    print('%(y)9.4f\n' % {'y': y[k]}, end='')
 


print('---------------------------------------------------------', end='');
print('-------------------------------------------------------\n', end='');

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

"""

–асчет характеристик фильтра дробной задержки производился программой resampling_spline_fd_filter.c или скриптом resampling_spline_fd_filter.py :


#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "dspl.h"

int main(int argc, char* argv[])
{
  double s0[8] = {1.0, 2.0, 2.0, 1.0, -0.5, -1.0, -2.0, -0.5}; /* исходный сигнал   */
  double *y = NULL;
  int    i, ny;
  
  void* hdspl; /* DSPL handle        */

  /* Загружаем libdspl.dll */
  hdspl = dspl_load();
  if(!hdspl)
  {
    printf("Cannot to load libdspl!\n");
    return 0;
  }
  
  farrow_spline(s0, 8, 1.0, 1.0, 0.25, &y, &ny);
  for(i = 0; i < ny; i++)
    printf("y[%d] = %+-6.3f\n", i, y[i]);

  /* Выгружаем libdspl.dll из памяти */
  dspl_free(hdspl);

  return 0;
}


# -*- coding: utf-8 -*-
"""
Created on Sun Aug 27 23:40:00 2017

@author: sergey
"""

import numpy as np

p = 1
q = 1
x0 = 0.25



s0 = np.array([1., 2., 2., 1., -0.5, -1., -2., -0.5])

tin = np.linspace(0.0, len(s0), len(s0), endpoint=False)

y = np.zeros(int(np.floor( float(len(s0)*q)/float(p))))

s = np.concatenate((np.array([0., 0.]), s0, np.array([0., 0.])))


print('---------------------------------------------------------', end='');
print('-------------------------------------------------------\n', end='');
print('k     x_k    n        d    s(n-3)    s(n-2)   s(n-1)    s(n)     ', end='');
print('a_0       a_1       a_2       a_3       y(k)\n', end='');
print('---------------------------------------------------------', end='');
print('-------------------------------------------------------\n', end='');

tout = np.zeros(int(np.floor( float(len(s0)*q)/float(p))))
for k  in range(len(y)):
    x =  k*q/p - x0       
    n =  int(x) + 1 + 3 
    d =  int(x) + 1 - x 

    tout[k] = x;
    
    a0 = s[n-1];
    a1 = 0.5*(s[n] - s[n-2]);
    a3 = 2.0*(s[n-2] - s[n-1]) + a1 + 0.5*(s[n-1] - s[n-3]);
    a2 = s[n-2] - s[n-1] +a3 + a1;
    
    y[k] = a0 - a1 * d + a2*d**2 - a3*d**3
    
    print('%(d0)d  %(f0)7.2f   %(d1)2d   %(f1)7.2f'  % 
        {'d0': k, 'f0': x, 'd1': n,  'f1': d}, end='')
    print('%(s0)8.1f %(s1)8.1f %(s2)8.1f %(s3)8.1f ' % 
        {'s0': s[n-3],'s1': s[n-2], 's2': s[n-1], 's3':s[n]}, end='')
  
    print('%(a0)9.3f %(a1)9.3f %(a2)9.3f %(a3)9.3f '% 
        {'a0': a0, 'a1': a1, 'a2': a2, 'a3': a3}, end='')
    print('%(y)9.4f\n' % {'y': y[k]}, end='')
 


print('---------------------------------------------------------', end='');
print('-------------------------------------------------------\n', end='');

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

"""

программы resampling_spline_interp.c и resampling_spline_interp.py расчета характеристик полиномиального интерполятора:


#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "dspl.h"

#define N 500

#define P 10
#define Q 1




int main(int argc, char* argv[])
{
  /* время             */
  double t0[8] = {0.0, 1.0, 2.0, 3.0,  4.0,  5.0,  6.0,  7.0};
  /* исходный сигнал   */
  double s0[8] = {1.0, 2.0, 2.0, 1.0, -0.5, -1.0, -2.0, -0.5};
  double *ti = NULL;
  double *y = NULL;
  int    i, ny, ni;
  
  double *hl = NULL;
  double *hs = NULL;
  double  h[5] = {0.0, 0.0, 1.0, 0.0, 0.0};
  
  double w[N], H[N];
  
  
  void* hdspl; /* DSPL handle        */
  void* hplot; /* GNUPLOT handle     */
  
  /* Загружаем libdspl.dll */
  hdspl = dspl_load();
  if(!hdspl)
  {
    printf("Cannot to load libdspl!\n");
    return 0;
  }
  
  for(i = 0; i < 8; i++)
    t0[i] *= (double)P;
  
  writetxt(t0, s0, 8, "dat/resample_spline_interp_s.txt");
  farrow_spline(s0, 8, P, Q, 0, &y, &ny);
  
  printf("ny = %d\n", ny);
  
  ti = (double*)malloc(ny*sizeof(double));
  linspace(0,ny,ny,DSPL_PERIODIC, ti);

  array_scale_lin(ti,   ny,  0, ny-1, 0, 14,  ti);
  array_scale_lin(y,    ny, -2,  2, -2, 4,   y);
  writetxt(ti, y, ny-1, "dat/resample_spline_interp_y.cm");
  
 
  farrow_lagrange(h, 5, P, Q, 0, &hl, &ni);
  ti = (double*)realloc(ti, ni*sizeof(double));
  linspace(0,ni,ni,DSPL_PERIODIC, ti);
  writetxt(ti, hl, ni, "dat/resample_lagrange_interp_filter_time.txt");
  
  farrow_spline(h, 5, P, Q, 0, &hs, &ni);
  ti = (double*)realloc(ti, ni*sizeof(double));
  linspace(0,ni,ni,DSPL_PERIODIC, ti);
  writetxt(ti, hs, ni, "dat/resample_spline_interp_filter_time.txt");
  
  
  linspace(0, M_PI, N, DSPL_PERIODIC, w);
  filter_freq_resp(hl, NULL, ni-1, w, N, DSPL_FLAG_LOGMAG, H, NULL, NULL);
  writetxt(w, H, N, "dat/resample_lagrange_interp_filter_freq.txt");
  
  linspace(0, M_PI, N, DSPL_PERIODIC, w);
  filter_freq_resp(hs, NULL, ni-1, w, N, DSPL_FLAG_LOGMAG, H, NULL, NULL);
  writetxt(w, H, N, "dat/resample_spline_interp_filter_freq.txt");
  
  /* Выгружаем libdspl.dll из памяти */
  dspl_free(hdspl);

  return 0;
}


# -*- coding: utf-8 -*-
"""
Created on Mon Aug 28 23:34:59 2017

@author: sergey
"""

import numpy as np
import scipy.signal as signal

def resampling_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):
        if(q==1):
            y = np.zeros(int(float((len(s)-1) * p) / float(q)) + 1)
        else:
           y = np.zeros(int( float(len(s) * p ) / float(q))) 
    else:
        y = np.zeros(int( float(len(s) * p ) / float(q)))
        
        

    t = np.zeros(len(y))
    s = np.concatenate((np.array([0., 0.]), s, np.array([0., 0.])))

    for k in range(len(y)):
        x = k*q/p - x0
        t[k] = x  
        n = int(np.floor(x)) + 4
        d = np.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] = a0 - a1 * d + a2*d**2 - a3*d**3  
    return y



def resampling_spline(s, p, q, x0):
    """
    """
    if(p>1):
        if(q==1):
            y = np.zeros(int(float((len(s)-1) * p) / float(q)) + 1)
        else:
           y = np.zeros(int( float(len(s) * p ) / float(q))) 
    else:
        y = np.zeros(int( float(len(s) * p ) / float(q)))
        
        

    t = np.zeros(len(y))
    s = np.concatenate((np.array([0., 0.]), s, np.array([0., 0.])))

    for k in range(len(y)):
        x = k*q/p - x0
        t[k] = x  
        n = int(np.floor(x)) + 4
        d = np.floor(x) + 1 - x  
        a0 = s[n-1];
        a1 = 0.5*(s[n] - s[n-2]);
        a3 = 2.0*(s[n-2] - s[n-1]) + a1 + 0.5*(s[n-1] - s[n-3]);
        a2 = s[n-2] - s[n-1] +a3 + a1;
        y[k] = a0 - a1 * d + a2*d**2 - a3*d**3  
    return y


p = 10
q = 1
x0 = 0

"""
Пример интерполяции сигнала с использованием фильтра Фарроу
"""
s = np.array([1., 2., 2., 1., -0.5, -1., -2., -0.5])

t = np.linspace(0, len(s)*p, len(s), endpoint = False, dtype = 'float64')
np.savetxt("dat/resample_spline_interp_s.txt", np.transpose([t, s]), fmt="%+.9e")

y = resampling_spline(s, p, q, x0)

t = np.linspace(0, len(y), len(y), endpoint = False, dtype = 'float64')
np.savetxt("dat/resample_spline_interp_y.txt", np.transpose([t, y]), fmt="%+.9e") 







"""
Расчет импульсной характеристики h[k] и частотной характеристики H(w)
"""

s = np.array([0., 0., 1., 0., 0.])

h = resampling_lagrange(s, p, q, x0)

t = np.linspace(0, len(h), len(h), endpoint = False, dtype = 'float64')
np.savetxt("dat/resample_lagrange_interp_filter_time.txt", np.transpose([t, h]), fmt="%+.9e")

w, H = signal.freqz(h)
H = 20.0 * np.log10(np.abs(H))
np.savetxt("dat/resample_lagrange_interp_filter_freq.txt", np.transpose([w, H]), fmt="%+.9e")


h = resampling_spline(s, p, q, x0)

t = np.linspace(0, len(h), len(h), endpoint = False, dtype = 'float64')
np.savetxt("dat/resample_spline_interp_filter_time.txt", np.transpose([t, h]), fmt="%+.9e")

w, H = signal.freqz(h)
H = 20.0 * np.log10(np.abs(H))
np.savetxt("dat/resample_spline_interp_filter_freq.txt", np.transpose([w, H]), fmt="%+.9e")









    
    
    

программы resampling_spline_pq.c и resampling_spline_pq.py расчета данных примера цифровой передискретизации:


#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "dspl.h"

#define N 54
#define P 20
#define Q 11
#define FS 26.4
#define L 10
#define F0 6.0

int main(int argc, char* argv[])
{
  double tc[N*L], c[N*L], ts[N], s[N];
  double *y=NULL;
  double *ty=NULL;
  double *err = NULL;
  double e0, es, el;
  int i, ny;
  void* hdspl; /* DSPL handle        */
  
  /* Загружаем libdspl.dll */
  hdspl = dspl_load();
  if(!hdspl)
  {
    printf("Cannot to load libdspl!\n");
    return 0;
  }
  
  linspace(0, (double)N/FS, N*L, DSPL_PERIODIC, tc);
  linspace(0, (double)N/FS, N,   DSPL_PERIODIC, ts);
  
  for(i = 0; i < N; i++)
    s[i] = sin(M_2PI * F0 * ts[i]);

  for(i = 0; i < N*L; i++)
    c[i] = sin(M_2PI * F0 * tc[i]);
  
  farrow_lagrange(s, N, P, Q, 0, &y, &ny);
  
  ty  = (double*)malloc(ny*sizeof(double));
  err = (double*)malloc(ny*sizeof(double));
  e0 = 0.0; 
  el = 0.0;
  for(i = 0; i < ny; i++)
  {
    ty[i]  = (double)(Q * i) / FS / (double)P;
    err[i] = sin(M_2PI * F0 * ty[i]) - y[i];
    e0 += sin(M_2PI * F0 * ty[i]) * sin(M_2PI * F0 * ty[i]);
    el += err[i] * err[i];
  }
  printf("E0 = %8.3f, el = %8.3f, NMSE lagrange = %9.6f\n", e0, el, el/e0);
  
  
  array_scale_lin(ts,   N,  0, 2,  0, 16,  ts);
  array_scale_lin( s,   N, -1, 1, -1,  2,   s); 
  writetxt(ts, s,       N, "dat/resample_lagrange_ex_fs_s.cm"); 

  array_scale_lin(tc,   N*L,  0, 2,  0, 16,  tc);
  array_scale_lin(c,    N*L, -1, 1, -1,  2,   c);
  writetxt(tc, c, (N-1)*L, "dat/resample_lagrange_ex_fs_c.cm");
  
  array_scale_lin(ty,   ny,  0, 2,  0, 16,  ty);
  array_scale_lin(y,    ny, -1, 1, -1,  2,   y);
  writetxt(ty, y,       ny, "dat/resample_lagrange_ex_fs_y.cm");
  
  array_scale_lin(err,  ny, -0.15, 0.15, -1,  2,   err);
  writetxt(ty, err,     ny, "dat/resample_lagrange_ex_fs_e.cm");
  
  linspace(0, (double)N/FS, N*L, DSPL_PERIODIC, tc);
  linspace(0, (double)N/FS, N,   DSPL_PERIODIC, ts);
  for(i = 0; i < N; i++)
    s[i] = sin(M_2PI * F0 * ts[i]);

  for(i = 0; i < N*L; i++)
    c[i] = sin(M_2PI * F0 * tc[i]);
  
  
  farrow_spline(s, N, P, Q, 0, &y, &ny);
  
  e0 = 0.0;
  es = 0.0;
  for(i = 0; i < ny; i++)
  {
    ty[i]  = (double)(Q * i) / FS / (double)P;
    err[i] = sin(M_2PI * F0 * ty[i]) - y[i];
    e0 += sin(M_2PI * F0 * ty[i]) * sin(M_2PI * F0 * ty[i]);
    es += err[i] * err[i];
  }
  printf("E0 = %8.3f, es = %8.3f, NMSE spline =   %9.6f\n", e0, es, es/e0);

  
  
  array_scale_lin(ts,   N,  0, 2,  0, 16,  ts);
  array_scale_lin( s,   N, -1, 1, -1,  2,   s); 
  writetxt(ts, s,       N, "dat/resample_spline_ex_fs_s.cm"); 

  array_scale_lin(tc,   N*L,  0, 2,  0, 16,  tc);
  array_scale_lin(c,    N*L, -1, 1, -1,  2,   c);
  writetxt(tc, c, (N-1)*L, "dat/resample_spline_ex_fs_c.cm");
  
  array_scale_lin(ty,   ny,  0, 2,  0, 16,  ty);
  array_scale_lin(y,    ny, -1, 1, -1,  2,   y);
  writetxt(ty, y,       ny, "dat/resample_spline_ex_fs_y.cm");
  
  array_scale_lin(err,  ny, -0.15, 0.15, -1,  2,   err);
  writetxt(ty, err,     ny, "dat/resample_spline_ex_fs_e.cm");
  
  
  if(y)
    free(y);
  if(ty)
    free(ty);
  if(err)
    free(err);
  
  /* Выгружаем libdspl.dll из памяти */
  dspl_free(hdspl);

  return 0;
}


# -*- coding: utf-8 -*-
"""
Created on Sat Sep  2 12:18:52 2017

@author: sergey
"""

# -*- coding: utf-8 -*-
"""
Created on Mon Aug 28 23:34:59 2017

@author: sergey
"""

import numpy as np
import scipy.signal as signal

def resampling_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):
        if(q==1):
            y = np.zeros(int(float((len(s)-1) * p) / float(q)) + 1)
        else:
           y = np.zeros(int( float(len(s) * p ) / float(q))) 
    else:
        y = np.zeros(int( float(len(s) * p ) / float(q)))
        
        

    t = np.zeros(len(y))
    s = np.concatenate((np.array([0., 0.]), s, np.array([0., 0.])))

    for k in range(len(y)):
        x = k*q/p - x0
        t[k] = x  
        n = int(np.floor(x)) + 4
        d = np.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] = a0 - a1 * d + a2*d**2 - a3*d**3  
    return y
    
    

def resampling_spline(s, p, q, x0):
    """
    """
    if(p>1):
        if(q==1):
            y = np.zeros(int(float((len(s)-1) * p) / float(q)) + 1)
        else:
           y = np.zeros(int( float(len(s) * p ) / float(q))) 
    else:
        y = np.zeros(int( float(len(s) * p ) / float(q)))
        
        

    t = np.zeros(len(y))
    s = np.concatenate((np.array([0., 0.]), s, np.array([0., 0.])))

    for k in range(len(y)):
        x = k*q/p - x0
        t[k] = x  
        n = int(np.floor(x)) + 4
        d = np.floor(x) + 1 - x  
        a0 = s[n-1];
        a1 = 0.5*(s[n] - s[n-2]);
        a3 = 2.0*(s[n-2] - s[n-1]) + a1 + 0.5*(s[n-1] - s[n-3]);
        a2 = s[n-2] - s[n-1] +a3 + a1;
        y[k] = a0 - a1 * d + a2*d**2 - a3*d**3  
    return y




P = 20
Q = 11
N = 54
L = 10
Fs = 26.4
f0 = 6.0

tc = np.linspace(0, N*L, N*L, endpoint = False, dtype = 'float64')/(Fs * float(L))
c = np.sin(np.pi*2.0*f0*tc)

ts = np.linspace(0, N, N, endpoint = False, dtype = 'float64')/Fs
s = np.sin(np.pi*2.0*f0*ts)

yL = resampling_lagrange(s, P, Q, 0.0)
yS = resampling_spline(s, P, Q, 0.0)
ty = float(Q) * np.linspace(0, len(yL), len(yL), endpoint = False, dtype = 'float64') / (Fs * float(P))

eL = yL - np.sin(np.pi*2.0*f0*ty)
eS = yS - np.sin(np.pi*2.0*f0*ty)

np.savetxt("dat/resample_lagrange_ex_fs_c.txt", np.transpose([tc[0:(N-1)*L], c[0:(N-1)*L] ]), fmt="%+.9e")
np.savetxt("dat/resample_lagrange_ex_fs_s.txt", np.transpose([ts, s]), fmt="%+.9e")
np.savetxt("dat/resample_lagrange_ex_fs_y.txt", np.transpose([ty[0:len(ty)-2], yL[0:len(ty)-2]]), fmt="%+.9e")
np.savetxt("dat/resample_lagrange_ex_fs_e.txt", np.transpose([ty[0:len(ty)-2],eL[0:len(ty)-2]]), fmt="%+.9e")
np.savetxt("dat/resample_spline_ex_fs_y.txt", np.transpose([ty[0:len(ty)-2], yS[0:len(ty)-2]]), fmt="%+.9e")
np.savetxt("dat/resample_spline_ex_fs_e.txt", np.transpose([ty[0:len(ty)-2], eS[0:len(ty)-2]]), fmt="%+.9e")


nmseL = np.sum(eL**2)/np.sum(np.sin(np.pi*2.0*f0*ty)**2)
nmseS = np.sum(eS**2)/np.sum(np.sin(np.pi*2.0*f0*ty)**2)

#print("nmseL = " ,nmseL) 
#print("nmseS = ", nmseS)   
    
    

Список литературы
[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] Franck A. Efficient Algorithms for Arbitrary Sample Rate Conversion with Application to Wave Field Synthesis. PhD thesis. Universitätsverlag Ilmenau, Ilmenau, 2012.

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

Последнее изменение страницы: 20.01.2024 (19:27:58)
Страница создана Latex to HTML translator ver. 5.20.11.14