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

Содержание

DSPL-2.0 is free digital signal processing algorithms library

Distributed under v3 LGPL v3 license

GitHub project page.

Found a mistake? Select it with the mouse and press ctrl+enter

Введение

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

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

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

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

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

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

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

Figure 2
Figure 2. Пересчет индексов отсчетов входного сигнала в интервал 

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

Везде далее в данном разделе, так же как и в предыдущем разделе, переменная  индексирует отсчеты входного сигнала  , а переменная  индексирует отсчеты выходного сигнала  после ресэмплера.

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

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

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

(1)

Для расчета кубического полинома нам требуется четыре отсчета входного сигнала  и , как это показано на рисунке 2а. Причем два отсчета  и  должны располагаться правее  , а два других отсчета  и  - левее.

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

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

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

(3)

На рисунке 3 показан пример пересчета индексов входного и выходного сигналов для -го выходного отсчета  при .

Figure 3
Figure 3. Пример пересчета индексов при

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

(4)

Тогда индекс , соответственно (2), равен:

(5)
и параметр  по шкале  равен:
(6)

Таким образом, необходимо по отсчетам , , и рассчитать коэффициенты кубического полинома  согласно формулам:

(7)
после чего значение  можно получить как значение интерполяционного полинома   для , т.е. .

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

Пусть входной сигнал  содержит  отсчетов, как это показано на рисунке 4.

Figure 4
Figure 4. Входной сигнал фильтра Фарроу

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

В таблице 1 приведен процесс компенсации дробной задержки сигнала , , показанного на рисунке 4.

Таблица 1
Table 1. Компенсация дробной задержки сигнала 

Значения  и  рассчитываются согласно (1), (2) и (3) соотвественно. Далее приведены значения  для каждого значения . Рамкой обведены нулевые значения  для  и , выходящие за пределы индексации сигнала . В следующих столбцах таблицы 1 приведены значения коэффициентов кубического полинома , рассчитанные согласно (7) для текущего момента . Наконец, в последнем столбце значение сигнала на выходе фильтра Фарроу.

На рисунке 5 показан исходный сигнал , а также сигнал  после компенсации дробной задержки .

Figure 5
Figure 5. Входной сигнал  фильтра Фарроу и сигнал после компенсации дробной задержки 

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

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

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

Амплитудно-частотная характеристика (АЧХ) и групповая задержка  полученных фильтров показаны на рисунке 7.

Figure 7
Figure 7. АЧХ и групповая задержка фильтров Фарроу для различного значения дробной задержки 

Из графиков АЧХ и групповой задержки можно сделать вывод, что коррекция задержки возможна только в нормированной полосе до  (т.е. только до Гц). Для более широкой полосы коррекции фильтр Фарроу на основе полиномиальной кусочно-кубической интерполяции не подходит ввиду недопустимо высокого искажения характеристик фильтра. Если требуется использовать компенсацию дробной задержки в более широкой полосе, то необходимо применить более сложные фильтры, использующие полиномиальную интерполяцию более высокого порядка, или альтернативные КИХ-фильтры, или БИХ-фильтры коррекции групповой задержки сигналов.

Скрипт resampling_lagrange_fd.py реализующий расчет таблицы 1 и данных для построения рисунка 5, написан на языке Python. Исходный код скрипта приведен в приложении.

В библиотеке DSPL-2.0 реализована функция farrow_lagrange, для дробной передискретизации сигналов.

Расчет семейства импульсных характеристик, АЧХ и групповой задержки фильтров Фарроу для различного значения дробной задержки , показанная на рисунках 6 и 7, производился при помощи библиотеки DSPL-2.0. Исходный код программы resampling_lagrange_fd_filter.c расчета данных для построения семейства характеристик фильтров Фарроу дробной задержки приведен в приложении.

Также в приложении можно найти аналогичный скрипт resampling_lagrange_fd_filter.py на языке Python.

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

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

Пусть входной сигнал , содержит  отсчетов сигнала, как это показано на рисунке 8а. Мы уже использовали данный сигнал в качестве входного сигнала цифрового фильтра коррекции дробной задержки.

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

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

(8)

Пересчет значений  и  также производится согласно (1), (2) и (3) соотвественно.

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

Figure 8
Figure 8. Цифровая интерполяция при использовании фильтра Фарроу:
a — входной сигнал ; б — результат интерполяции

Фильтр цифровой интерполяции можно трактовать как фильтр нижних частот.

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

Figure 9
Figure 9. Импульсная характеристика и АЧХ интерполятора Фарроу,

Пунктирной линией показана импульсная характеристика при .

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

Расчет характеристик интерполятора на основе фильтра Фарроу также производился с помощью библиотеки DSPL-2.0.

Исходный код программы resampling_lagrange_interp.c приведен в приложении. Также в приложении приведен аналогичный скрипт resampling_lagrange_interp.py на языке Python.

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

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

(9)
с частотой  кГц содержит  отсчета взятых с частотой дискретизации  кГц.

Отношение  частоты дискретизации  к частоте сигнала  не является целым числом. В результате на одном периоде синусоидального сигнала не укладывается целое число отсчетов, как это показано на верхнем графике рисунка 10.

Произведем изменение частоты дискретизации входного сигнала в  раз, где  и получим сигнал , с частотой дискретизации  кГц. В результате на одном периоде синусоидального сигнала (9) будет укладываться ровно 8 отсчетов . Дробное изменение частоты дискретизации будем производить при использовании фильтра Фарроу на основе полиномиальной интерполяции. Пересчет значений  и  производится согласно (1), (2) и (3) соотвественно.

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

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

Тонкими пунктирными линиями на рисунке 10 показаны исходный непрерывный синусоидальный сигнал (9).

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

Выводы

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

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

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

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

Смотри также

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

Приложения

Данные для построения рисунков данного раздела были просчитаны при использовании библиотеки 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 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 *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_lagrange(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]
    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
    
    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_lagrange_fd_filter.c или скриптом resampling_lagrange_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 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 *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_lagrange(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]
    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
    
    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_lagrange_interp.c и resampling_lagrange_interp.py расчета характеристик полиномиального интерполятора:


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

#define N 200

#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 *hi = 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_lagrange_interp_s.txt");
  farrow_lagrange(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_lagrange_interp_y.cm");
  
 
  farrow_lagrange(h, 5, P, Q, 0, &hi, &ni);
  ti = (double*)realloc(ti, ni*sizeof(double));
  linspace(0,ni,ni,DSPL_PERIODIC, ti);
  writetxt(ti, hi, ni, "dat/resample_lagrange_interp_filter_time.txt");
  
  
  linspace(0, M_PI, N, DSPL_PERIODIC, w);
  filter_freq_resp(hi, NULL, ni-1, w, N, DSPL_FLAG_LOGMAG, H, NULL, NULL);
  writetxt(w, H, N, "dat/resample_lagrange_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



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_lagrange_interp_s.txt", np.transpose([t, s]), fmt="%+.9e")

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

t = np.linspace(0, len(y), len(y), endpoint = False, dtype = 'float64')
np.savetxt("dat/resample_lagrange_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")









    
    
    

Программы resampling_lagrange_pq.c и resampling_lagrange_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;
  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));
  for(i = 0; i < ny; i++)
    ty[i] = (double)(Q * i) / FS / (double)P;
  
  
  array_scale_lin(ts,   N,  0, 2,  0, 16,  ts);
  array_scale_lin( s,   N, -1, 1, -2,  4,   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, -2,  4,   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, -2,  4,   y);
  writetxt(ty, y,      ny, "dat/resample_lagrange_ex_fs_y.cm");
  
  
  if(y)
    free(y);
  if(ty)
    free(ty);
  
  /* Выгружаем 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



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)

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

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, y]), fmt="%+.9e")


    
    
    

Reference

[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.


Page update: 19.06.2020 (13:45:35)