Задержка как система свертки

Целые числа целочисленных задержек

Рассмотрим задержку цифрового сигнала, $y[n]=x[n-D]$ где D - целое число. Эта операция может быть представлена как фильтр свертки $$y=h*x$$, с конечной импульсной характеристикой $\mathit{h}\left[\mathit{n}\right]=\delta \left[\mathit{n}-\mathit{D}\right]$. Соответствующая передаточная функция $$H(z)=z^{-D}$$, и частая реакция является $\mathit{H}\left(\omega \right)={\mathit{e}}^{-\mathit{i}\omega \mathit{D}}$. Программно можно реализовать такой целочисленный фильтр задержки с помощью следующего кода MATLAB ®.

% Create the FIR
D = 3; % Delay value
h = [zeros(1,D) 1]

h = 1×4

     0     0     0     1

Сдвиньте последовательность путем фильтрации ее через конечную импульсную характеристику h. Обратите внимание на начальные нули в начале выхода, которые означают начальное условие, которое присуще таким фильтрам.

x = (1:10)';
dfir = dsp.FIRFilter(h);
y = dfir(x)'

y = 1×10

     0     0     0     1     2     3     4     5     6     7

Нецелочисленные задержки через D/A интерполяцию

Задержка последовательности $\mathit{x}\left[\mathit{n}-\mathit{D}\right]$ не определяется всякий раз, когда D не является целым числом. Чтобы сделать такие дробные задержки разумными, необходимо добавить промежуточный D/A каскад интерполяции, чтобы дискретизировать выход на contniuum. То есть, $\mathit{y}\left[\mathit{n}\right]=\widehat{{\mathit{x}}_{\mathit{n}}}\left(\mathit{n}-\mathit{D}\right)$ где $\widehat{{\mathit{x}}_{\mathit{n}}}$ обозначает некоторую D/A интерполяцию входного секвенирования $\mathit{x}$. Интерполяционная функция D/A$\widehat{{\mathit{x}}_{\mathit{n}}}$ может зависеть от n и может рассматриваться как представление базовой модели аналогового сигнала, от которой зависит последовательность $\mathit{x}$ был взят образец. Эта стратегия используется в других задачах повторной дискретизации, таких как преобразование скорости.

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

Полосно-ограниченные фильтры дробной задержки

Формула интерполяции Шеннона-Уиттакера $\hat{\mathit{x}}\left(\mathit{t}\right)=\sum _{\mathit{k}}\mathit{x}\left[\mathit{k}\right]\mathrm{sinc}\left(\mathit{t}-\mathit{k}\right)$ моделирует полосно-ограниченные сигналы. То есть промежуточное D/A преобразование $\hat{\mathit{x}}$ является полосно-ограниченной реконструкцией входа последовательности. Для значения задержки D, дробной задержки $\mathit{y}\left[\mathit{n}\right]=\hat{\mathit{x}}\left(\mathit{n}-\mathit{D}\right)$, который использует то же самое $\hat{\mathit{x}}$ для каждого n может быть представлен как фильтр свертки. Этот фильтр называется идеальным полосно-ограниченным дробным фильтром задержки, и его импульсная характеристика

$$h_D[k]=\text{sinc}(k-D)$$.

Соответствующая частотная характеристика (то есть DTFT) задается как $$H_D(\omega )=e^{-i\omega D}$$.

Причинно-следственная конечная импульсная характеристика Приближения идеального полосно-ограниченного сдвигового фильтра

Идеальный фильтр сдвига синуса, описанный в предыдущем разделе, является полнопроходным фильтром (т.е. $\left|{\mathit{H}}_{\mathit{d}}\left(\omega \right)\right|=1$), но он имеет бесконечную и не причинно-следственную импульсную характеристику ${\mathit{h}}_{\mathit{D}}$. В MATLAB он не может быть представлен как вектор, а скорее как функция от индекса k.

% Ideal Filter sequence
D = 0.4;
hIdeal = @(k) sinc(k-D); 

Для практических и вычислительных целей идеальный фильтр может быть усечен в конечном индексе окне за счет некоторой потери полосы пропускания. Для целевого значения задержки $\mathit{D}$ и желаемую длину N, окно индексов $\mathit{k}$удовлетворение $|\mathit{D}-\mathit{k}|\le \frac{\mathit{N}}{2}$ симметрично около $\mathit{D}$, и захватывает основную лепестку идеального фильтра. Для $\mathit{D}={\mathit{i}}_0+\mathit{FD}$ где $\mathit{0}\le \mathit{FD}\le 1$ и целое число ${\mathit{i}}_0$, явные индексы окна$$\{i_0-\lfloor \frac{N-1}{2}\rfloor ,\dots ,i_0+\lfloor \frac{N}{2}\rfloor \}$$. Целое число ${\mathit{i}}_0$ упоминается как целочисленная задержка и может быть выбрана произвольно. Чтобы сделать конечную импульсную характеристику причинно-следственную связь, установите ${\mathit{i}}_0=\lfloor \frac{\mathit{N}-1}{2}\rfloor$, поэтому окно индекса $\left\{0,\dots ,\mathit{N}-1\right\}$. Приведенный ниже код изображает обоснование причинно-следственной конечной импульсной характеристики приближения.

% FIR approximation with causal shift
N = 6;
idxWindow = (-floor((N-1)/2):floor(N/2))';
i0 = -idxWindow(1); % Causal latency
hApprox = hIdeal(idxWindow);
plot_causal_fir("sinc",D,N,i0,hApprox,hIdeal);

Усечение фильтра sinc вызывает пульсацию в частотной характеристики, которая может быть решена путем применения весов $\left\{{\mathit{w}}_{\mathit{k}}\right\}$ (такой как Кайзер или Хемминг) к КИХ коэффициентам.

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

Вы можете спроектировать такой фильтр, используя designFracDelayFIR функции и dsp.VariableFractionalDelay Системный объект в 'FIR' режим, в обоих из которых используются веса окон Кайзера.

Фракционные фильтры задержки на основе лагранжа

Фракционные фильтры задержки на основе Лагранжа используют полиномиальную аппроксимацию на движущемся окне входа отсчетов. То есть, $\widehat{{\mathit{x}}_{\mathit{n}}}\left(\mathit{t}\right)$ является полиномом некоторой фиксированной степени K. Как и основанные на sinc фильтры задержки, основанные на Лагранже фильтры задержки могут быть сформулированы как причинная свертка конечной импульсной характеристики (т.е. $\mathit{y}=\mathit{h}*\mathit{x}$) длины N = K + 1 и поддерживается в окне индекса$$\{-\lfloor \frac{N-1}{2}\rfloor ,\dots \lfloor \frac{N}{2}\rfloor \}$$. Подобно модели на основе синуса, примените причинно-следственную задержку ${\mathit{i}}_0=\lfloor \frac{\mathit{N}-1}{2}\rfloor$. Учитывая дробную задержку $\mathit{FD}$, коэффициенты конечной импульсной характеристики $\mathit{h}\left[0\right],\dots ,\mathit{h}\left[\mathit{K}\right]$ (причинно-сдвинутый) фильтр Лагранжа задержки может быть получен путем решения системы линейных уравнений, как написано ниже. Эти уравнения описывают стандартную задачу полинома Лагранжа.

Вот, ${\mathit{t}}_0,\dots ,{\mathit{t}}_{\mathit{K}}$ являются перечисленными индексами окна выборки. Реализация проста.

% Filter parameters
FD = 0.4;
K = 7;    % Polynomial degree
N = K+1;  % FIR Length
idxWindow = (-floor((N-1)/2):floor(N/2))';

% Define and solve Lagrange interpolation equations
V = idxWindow.^(0:K); % Vandermonde structure
C = FD.^(0:K);

hLagrange = C/V;  % Solve for the coefficients
i0 = -idxWindow(1); % Causal latency
plot_causal_fir("Lagrange",FD,N,i0,hLagrange);

Эта модель может быть реализована как фильтр конечной импульсной характеристики прямой формы, если значение задержки FD фиксировано, или с помощью структуры Farrow, если значение задержки изменяется. Ниже приведен раздел, посвященный реализации интерполяции Лагранжа с использованием dsp.VariableFractionalDelay в 'Farrow' режим.

Проектируйте и реализуйте основанные на синусоидальной задержке фильтры конечной импульсной характеристики задержки

Следующий раздел посвящен разработке и реализации фракционных фильтров задержки на основе sinc.

Функциональная designFracDelayFIR в режиме Проект

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

FD = 0.32381;
N = 10;
h = designFracDelayFIR(FD,N)

h = 1×10

    0.0046   -0.0221    0.0635   -0.1664    0.8198    0.3926   -0.1314    0.0552   -0.0200    0.0042

Сама реализация может быть выполнена как стандартный конечная импульсная характеристика, такой как dsp.FIRFilter Системный объект.

% Create an FIR filter object
fdfir = dsp.FIRFilter(h); 

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

% Generate some input
n = (1:100)';
x = gen_input_signal(n);

% Filter the input signal
y = fdfir(x);
plot_sequences(n,x, n,y);
legend('Filter Output','Original Sequence')
title('Raw Filter Output v.s. Input Sequence')

release(fdfir);

Заметьте, что фактическая задержка фильтра не $\mathit{FD}$, но скорее $\mathit{FD}+{\mathit{i}}_0$ из-за причинно-следственной целочисленной задержки ${\mathit{i}}_0$. Эта задержка возвращается из designFracDelayFIR функция также. Сдвиньте график последовательности входа на $\mathit{FD}+{\mathit{i}}_0$ для выравнивания выходов фильтра по ожидаемому результату.

[h,i0] = designFracDelayFIR(FD,N);
y = fdfir(x);
plot_sequences(n+i0+FD,x, n,y);
legend('Filter Output','Input Sequence (shifted by FD+i0)')
title('Filter Output v.s. Time Adjusted Input Sequence')

Обратите внимание, что сдвинутые входные маркеры, расположенные в $\left(\mathit{n}+\mathit{FD}+{\mathit{i}}_0,\mathit{x}\left[\mathit{n}\right]\right)$ обычно не совпадают с маркерами выходных выборок $\left(\mathit{n},\mathit{y}\left[\mathit{n}\right]\right)$, потому что $\mathit{n}+\mathit{FD}+{\mathit{i}}_0$ падает на нецелочисленные значения на оси X, тогда как n целого числа. Скорее сдвинутые входные выборки падают приблизительно на линию, соединяющую каждые два последовательных выхода выборки.

plot_sequences(n+i0+FD,x, n,y,'with line');
legend('Filter Output','Input Sequence (shifted by FD+i0)')
title('Output Samples v.s. Shifted Input Samples ')
xlim([20,30])

The dsp.VariableFractionalDelay Системный объект в 'FIR' Способ

Аналогично designFracDelayFIR, а dsp.VariableFractionalDelay объект может также проектировать основанные на sinc фильтры задержки при использовании с 'FIR' режим интепроляции. Начните с создания образца Системного объекта. Длина конечной импульсной характеристики всегда четная и задается как параметр половинной длины.

vfd_fir = dsp.VariableFractionalDelay('InterpolationMethod','FIR','FilterHalfLength',N/2);
i0_vfd_fir = vfd_fir.FilterHalfLength;    % Interger latency

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

y = vfd_fir(x,i0+FD);
release(vfd_fir)
plot_sequences(n+i0+FD,x, n,y);
legend('Filter Output','Input Sequence (shifted by FD+i0)')
title('dsp.VariableFractionalDelay in FIR Mode')

Сравнение designFracDelayFIR и dsp.VariableFractionalDelay в режиме ' конечной импульсной характеристике'

Оба designFracDelayFIR и dsp.VariableFractionalDelay в 'FIR' mode обеспечивают основанные на sinc фильтры дробной задержки, но их реализации различаются.

Использование designFractionalDelayFIR предпочтительнее dsp.VariableFractionalDelay в 'FIR' СПОСОБ ДЛЯ ЕГО ПРОСТОТЫ, ПОВЫШЕНИЯ ПРОИЗВОДИТЕЛЬНОСТИ И ЭФФЕКТИВНОСТИ. На рисунке ниже фильтр разработан с dsp.VariableFractionalDelay имеет более короткую полосу пропускания, и ее групповая задержка отключена на ~ 0,02 от номинального значения .

% Obtain the FIR coefficients from the dsp.VariableFractionalDelay object
h_vfd_fir = vfd_fir([1;zeros(31,1)],i0_vfd_fir+FD);
release(vfd_fir);
plot_freq_and_gd(h,i0,[],"designFracDelayFIR", h_vfd_fir,i0_vfd_fir,[],"dsp.VariableFractionalDelay FIR mode");
hold on;
yline(FD,'DisplayName','Target Fractional Delay');
ylim([-0.1,0.4])

Проектирование и реализация основанных на Лагранже фильтров задержки

Фракционный фильтр задержки на основе Лагранжа является в вычислительном отношении дешевым и может быть эффективно реализован с использованием структуры Фэрроу. Фильтр Фэрроу является специальным типом конечной импульсной характеристикой, который реализован с использованием только элементарных алгебраических операций, таких как скалярные сложения и умножения. В отличие от проектов на основе sinc, фильтры Фэрроу не требуют специализированных функций (таких как sinc или Bessel), чтобы вычислить коэффициенты конечной импульсной характеристики задержки. Это делает фильтры дробной задержки Farrow особенно простыми в реализации на основном оборудовании.

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

Системный объект dsp.VariableFractionalDelay в 'Farrow' mode

Используйте системный объект dsp.VariableFractionalDelay в 'Farrow' режим для создания и реализации фильтров задержки Farrow. Начните с создания образца системного объекта:

vfd = dsp.VariableFractionalDelay('InterpolationMethod','Farrow','FilterLength',8);
i0var = floor(vfd.FilterLength/2)  % Interger latency of the filter

i0var = 4

Примените созданный объект к входному сигналу и постройте график результата.

y = vfd(x,i0var+FD);
plot_sequences(n+i0var+FD,x, n,y);
legend('Farrow Fractional Delay Output','Input Sequence (shifted by FD+i0)')
title('dsp.VariableFractionalDelay in Farrow Mode')

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

release(vfd)
FDs = i0var+5*(0:0.2:0.8); % Fractional delays vector
xsource = dsp.SignalSource(x,20);
ysink = dsp.AsyncBuffer;
for FD=FDs
    xk = xsource();
    yk = vfd(xk, FD);
    write(ysink,yk);
end
y = read(ysink);

plot_sequences(n+i0var,x, n,y);
legend('Variable Fractional Delay Output','Original Sequence (shifted by i0)')
title('dsp.VariableFractionalDelay in Farrow Mode, Varying Delay')

Пропускная способность конечной импульсной характеристики фильтров дробной задержки: анализ и Проект

Более длинные фильтры дают лучшее приближение идеального фильтра задержки. Действительно, с точки зрения необработанных квадратичных норм это так. Однако нам нужна метрика, которая является более практически значимой, например, полоса пропускания. Функция designFracDelayFIR измеряет комбинированную полосу пропускания, которая определяется как частотная область значений, в котором и коэффициент усиления, и групповая задержка находятся в пределах 1% от их номинальных значений. Измеренная комбинированная полоса пропускания может быть получена как возврат значение designFracDelayFIR функция. Сравните фильтр длины 16 (синяя) с фильтром длины 256 (красная) на рисунке ниже. Как ожидалось, более длинный фильтр имеет значительно более высокую комбинированную полосу пропускания.

FD = 0.3;
N1 = 16;
N2 = 256;
[h1,i1,bw1] = designFracDelayFIR(FD, N1);
[h2,i2,bw2] = designFracDelayFIR(FD, N2);

plot_freq_and_gd(h1,i1,bw1,"N="+num2str(N1), h2,i2,bw2,"N="+num2str(N2));
ylim([-0.2,0.6])

Функциональная designFracDelayFIR в режиме проекта полосы пропускания

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

FD = 0.3;
bwLower = 0.9; % Target bandwidth lower limit
[h,i0fixed,bw] = designFracDelayFIR(FD,bwLower);
fdfir = dsp.FIRFilter(h);
info(fdfir)

ans = 6x35 char array
    'Discrete-Time FIR Filter (real)    '
    '-------------------------------    '
    'Filter Structure  : Direct-Form FIR'
    'Filter Length     : 52             '
    'Stable            : Yes            '
    'Linear Phase      : No             '

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

Искажение в сигналах с высокой пропускной способностью

В этом разделе мы сравниваем эффективность двух проектов (длинный sinc против короткого Lagrange) с высоким входом в полосу пропускания. The dsp.VariableFractionalDelay в предыдущем разделе находится 8-градусная структура Фэрроу, фактически конечная импульсная характеристика длины 9. Фильтр, полученный designFracDelayFIR(FD,0.9) имеют длину 52 выборки. Сложение двух конечных импульсных характеристик частотных характеристик вместе на том же графике демонстрирует различие полосы пропускания между ними.

release(vfd);
hvar = vfd([1;zeros(31,1)],i0var+FD);
plot_freq_and_gd(h,i0fixed,bw,"Sinc-based", hvar,i0var,[],"Farrow");
ylim([-0.2,0.6])

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

n=(1:80)';
x = high_bw_signal(n);

y1 = fdfir(x);
y2 = vfd(x,i0var+FD);

plot_signal_comparison(n,x,y1,y2,h,hvar,i0fixed,i0var,FD);

Который следует использовать: dsp.VariableFractionalDelay или designFracDelayFIR ?

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