Exponenta Banner

Проект дециматоров и интерполяторов

Открыть скрипт

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

Роль Lowpass в преобразовании скорости

Преобразование скорости - это процесс изменения скорости дискретного сигнала для получения нового дискретного представления базового непрерывного сигнала. Способ включает равномерное понижение дискретизации и увеличение дискретизации. Равномерное понижение дискретизации со скоростью N относится к взятию каждой N-й выборки последовательности и отбрасыванию остальных выборок. Равномерное увеличение дискретизации в N раз относится к заполнению N-1 нулями между каждые двумя последовательными выборками.

x = 1:3
L = 3; % upsampling rate
M = 2; % downsampling rate

% Upsample and downsample
xUp = upsample(x,L)
xDown = downsample(x,M)

x =

     1     2     3


xUp =

     1     0     0     2     0     0     3     0     0


xDown =

     1     3

Обе эти основные операции вводят программные продукты: понижающая дискретизация вводит сглаживание, а повышающая дискретизация - визуализацию. Чтобы уменьшить этот эффект, используйте lowpass.

  • При понижающей дискретизации со скоростью, $N$lowpass фильтр, примененный до понижающей дискретизации, ограничивает входную полосу пропускания и, таким образом, устраняет сглаживание спектра. Это аналогично аналоговому LPF, используемому в A/D преобразователях. В идеале такой сглаживающий фильтр имеет единичное усиление и частоту отключения, $\omega_c = \frac{1}{N}\omega_N$вот$\omega_N$ частота Найквиста сигнала. Примечание: базовая частота дискретизации незначительна, мы принимаем нормированные частоты (то есть) $\omega_N=1$на протяжении всего обсуждения.

  • При увеличении дискретизации со скоростью, $N$lowpass фильтр, примененный после увеличения дискретизации, известен как антивизионный фильтр. Фильтр удаляет спектральные изображения низкоскоростного сигнала. В идеале частота отключения этого антивизуального фильтра является ($\omega_c = \frac{1}{N}$как и его антиалифицирующий аналог), в то время как его коэффициент усиления равен.$N$

Для операций увеличения дискретизации и понижающей дискретизации требуется$N$ lowpass с нормализованной частотой среза. $\frac{1}{N}$Единственное различие заключается в необходимом усилении и размещении фильтра (до или после преобразования скорости).

Комбинация увеличения дискретизации сигнала в множитель, $L$с последующей фильтрацией, а затем понижающей дискретизации в множитель$M$ преобразует частоту выборки последовательности в рациональный коэффициент. $\frac{L}{M}$Это достигается повышением дискретизации по скорости$L$ с последующей фильтрацией, а затем понижением дискретизации по скорости. Операция $M$преобразования порядка скорости не может быть заменена. Один фильтр, который объединяет сглаживание и сглаживание изображения, помещается между стадиями повышающей дискретизации и понижающей дискретизации. Этот фильтр является lowpass с нормированной частотой отключения и$\omega_c=\min(1/L,1/M)$ усилением.$L$

В то время как любой lowpass конечной импульсной характеристики проекта функцию (например fir1, firpm, или fdesign) мог бы спроектировать соответствующий фильтр сглаживания и анти-визуализации, функцию designMultirateFIR дает удобный и упрощенный интерфейс. В следующих нескольких разделах показано использование этих функций для разработки фильтра и продемонстрировать, почему designMultirateFIR является предпочтительным способом.

Фильтрованное преобразование скорости: дециматоры, интерполяторы и преобразователи рациональной скорости

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

Фильтрация преобразования скорости с использованием filter, upsample, и downsample функции

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

  1. Создайте сглаживающий lowpass фильтр h

  2. Фильтруйте вход, хотя h

  3. Понижайте фильтрованную последовательность в множитель M

% Define an input sequence
x = rand(60,1);

% Implement an FIR decimator
h = fir1(L*12*2,1/M); % an arbitrary filter
xDecim = downsample(filter(h,1,x), M);

Интерполяция относится к усилению дискретизации с последующей фильтрацией. Реализация очень похожа на десятикратное уменьшение.

xInterp = filter(h,1,upsample(x,L));

Наконец, рациональное преобразование скорости состоит из интерполятора, за которым следует дециматор (в этом конкретном порядке).

xRC = downsample(filter(h,1,upsample(x,L) ), M);

Фильтрация преобразования скорости с использованием системных объектов

Для потоковых данных системные объекты dsp.FIRInterpolator, dsp.FIRDecimator, и dsp.FIRRateConverter инкапсулируйте изменение скорости и фильтрацию в один объект. Для примера конструкции интерполятора выполняется следующим образом.

firInterp = dsp.FIRInterpolator(L,h);

Затем передайте последовательность к вновь созданному объекту шаговым вызовом.

xInterp = firInterp(x);

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

firDecim = dsp.FIRDecimator(M,h); % Construct
xDecim = firDecim(x); % Decimate (step call)

firRC = dsp.FIRRateConverter(L,M,h); % Construct
xRC = firRC(x); % Convert rate (step call)

Использование системных объектов обычно является предпочтительным, так как они:

  • Разрешает более чистый синтаксис.

  • Сохраните состояние в качестве фильтра начального конидитона для последующих вызовов шага.

  • Самое главное, они используют очень эффективный алгоритм полифазы.

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

Создайте фильтр преобразования скорости с помощью designMultirateFIR

Функция designMultirateFIR(L,M) автоматически находит подходящую частоту масштабирования и среза для заданного коэффициента преобразования скорости. $L/M$Используйте коэффициенты конечной импульсной характеристики, возвращенные designMultirateFIR с dsp.FIRDecimator (если)$L=1$, dsp.FIRInterpolator (если)$M=1$, или dsp.FIRRateConverter (общий случай).

Давайте разработаем интерполяционный фильтр:

L = 3;
bInterp = designMultirateFIR(L,1); % Pure upsampling filter
firInterp = dsp.FIRInterpolator(L,bInterp);

Затем примените интерполятор к последовательности.

% Create a sequence
n = (0:89)';
f = @(t) cos(0.1*2*pi*t).*exp(-0.01*(t-25).^2)+0.2;
x = f(n);

% Apply interpolator
xUp = firInterp(x);
release(firInterp);

Давайте сначала рассмотрим необработанный выход интерполятора и сравним с исходной последовательностью.

plot_raw_sequences(x,xUp);

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

  • Растягивается временной интервал (как и ожидалось).

  • Сигнал имеет задержку в половине длины конечной импульсной характеристики length(h)/2 (обозначается$i_0$ далее).

  • В начале есть переходный процесс.

Чтобы сравнить, выровнять и масштабировать временные интервалы двух последовательностей. Интерполированная выборка xUp[k] соответствует входному времени.$t[k]=\frac{1}{L}(k-i_0)$

nUp = (0:length(xUp)-1);
i0 = length(bInterp)/2;
plot_scaled_sequences(n,x,(1/L)*(nUp-i0),xUp,["Original Sequence",...
            "Interpolator Output Sequence (Time Adjusted)"],[0,60]);

Эта же идея работает для понижающей дискретизации, где преобразование времени является:$t[k] = Mk-i_0$

M = 3;
bDecim = designMultirateFIR(1,M); % Pure downsampling filter
firDecim = dsp.FIRDecimator(M,bDecim);
xDown = firDecim(x);

Постройте их в одной шкале и скорректируйте на задержку. Обратите внимание, что они отлично перекрываются.

i0 = length(bDecim)/2;
nDown = (0:length(xDown)-1);
plot_scaled_sequences(n,x,M*nDown-i0,xDown,["Original Sequence",...
            "Decimator Output Sequence (Time Adjusted)"],[-10,80]);

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

hfv = fvtool(firInterp,firDecim); % Notice the gains in the passband
legend(hfv,"Interpolation Filter L="+num2str(L), ...
            "Decimation Filter M="+num2str(M));

Общие рациональные преобразования могут быть обработаны так же, как увеличение дискретизации и понижение дискретизации. Срез есть$~\omega_c=\min(\frac{1}{L},~\frac{1}{M})~$, и коэффициент усиления есть. $L$Функция designMultirateFIR вычисляет это автоматически.

L = 5;
M = 2;
b = designMultirateFIR(L,M);
firRC = dsp.FIRRateConverter(L,M,b);

Теперь сравним комбинированный фильтр с отдельными компонентами интерполяции/децимации.

firDecim = dsp.FIRDecimator(M,designMultirateFIR(1,M));
firInterp = dsp.FIRInterpolator(L,designMultirateFIR(L,1));

hfv = fvtool(firInterp,firDecim, firRC); % Notice the gains in the passband
legend(hfv,"Interpolation Filter L="+num2str(L),...
    "Decimation Filter M="+num2str(M), ...
    "Rate Conversion Filter L/M="+num2str(L)+"/"+num2str(M));

Один раз в FIRRateConverter настраивается, выполняет преобразование скорости шаговым вызовом.

xRC = firRC(x);

Постройте график входа и выхода фильтра с регулировкой времени, заданной.$t[k] = \frac{1}{L}(Mk-i_0)$

nRC = (0:length(xRC)-1)';
i0 = length(b)/2;
plot_scaled_sequences(n,x,(1/L)*(M*nRC-i0),xRC,["Original Sequence",...
        "Rate Converter Output Sequence (time adjusted)"],[0,80]);

Настройка Lowpass конечной импульсной характеристики Расчётных параметров

Использование designMultirateFIR можно также настроить длину конечной импульсной характеристики, ширину перехода и затухание в полосе задерживания.

Настройка длины конечной импульсной характеристики

Длиной конечной импульсной характеристики можно управлять через L, M и третий параметр P, называемый половинно-полифазной длиной, значение по умолчанию которого 12 (см. Выходные аргументы для получения дополнительной информации). Давайте рассмотрим две проектные точки.

% Unspecified half-length defaults to 12
b24 = designMultirateFIR(3,1);

halfPhaseLength = 20;
b40 = designMultirateFIR(3,1,halfPhaseLength);

Обычно большая половина длины полифазы приводит к более крутым переходам.

hfv = fvtool(b24,1,b40,1);
legend(hfv, 'Polyphase length = 24 (Default)','Polyphase length = 40');

Настройка ширины перехода

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

TW = 0.02;
bTW = designMultirateFIR(3,1,TW);

hfv = fvtool(b24,1,bTW,1);
legend(hfv, 'Default Design (FIR Length = 72)',"Design with TW="...
        +num2str(TW)+" (FIR Length="+num2str(length(bTW))+")");

Особый случай преобразования скорости интерполяторами и дециматорами 2: Полуполоса

Используя полупериодический фильтр (т.е.)$\omega_c=\frac{1}{2}$, можно выполнить преобразование частоты дискретизации в 2 раза. The dsp.FIRHalfbandInterpolator и dsp.FIRHalfbandDecimator объекты выполняют интерполяцию и десятикратное уменьшение в 2 раза, используя полуполосные фильтры. Этот системный объект реализован с использованием эффективной полифазной структуры, специфичной для этого преобразования скорости. Аналоги БИХ dsp.IIRHalfbandInterpolator и dsp.IIRHalfbandDecimator может быть еще более эффективным. Эти системные объекты могут также работать с пользовательскими частотой дискретизации.

Визуализируйте величину ответ с помощью fvtool. В случае интерполяции фильтр сохраняет большую часть спектра от 0 до Fs/2 при ослаблении спектральных изображений. Для десятикратного уменьшения фильтр проходит примерно половину полосы, то есть от 0 до Fs/4, и ослабляет вторую половину порядка, чтобы минимизировать сглаживание. Величина ослабления может быть установлена на любое желаемое значение как для интерполяции, так и для десятикратного уменьшения. Если не задано, значение по умолчанию равно 80 дБ.

Fs = 1e6;
hbInterp = dsp.FIRHalfbandInterpolator('TransitionWidth',Fs/10,...
    'SampleRate',Fs);
fvtool(hbInterp) % Notice gain of 2 (6 dB) in the passband
hbDecim  = dsp.FIRHalfbandDecimator('TransitionWidth',Fs/10,...
    'SampleRate',Fs);
fvtool(hbDecim)

Эквириппл- Проект

Функция designMultirateFIR использует оконные проекты конечной импульсной характеристики lowpass. Также могут быть применены другие методы проекта lowpass, такие как equiripple. Для дополнительного контроля над процессом проекта используйте fdesign создания фильтра. Следующий пример проектирует дециматор, используя fdesign.decimator функция.

M   = 4;   % Decimation factor
Fp  = 80;  % Passband-edge frequency
Fst = 100; % Stopband-edge frequency
Ap  = 0.1; % Passband peak-to-peak ripple
Ast = 80;  % Minimum stopband attenuation
Fs  = 800; % Sampling frequency
fdDecim = fdesign.decimator(M,'lowpass',Fp,Fst,Ap,Ast,Fs) %#ok

fdDecim = 

  decimator with properties:

          MultirateType: 'Decimator'
               Response: 'Lowpass'
       DecimationFactor: 4
          Specification: 'Fp,Fst,Ap,Ast'
            Description: {4x1 cell}
    NormalizedFrequency: 0
                     Fs: 800
                  Fs_in: 800
                 Fs_out: 200
                  Fpass: 80
                  Fstop: 100
                  Apass: 0.1000
                  Astop: 80

Спецификации для фильтра определяют, что диапазон переходов 20 Гц является приемлемым между 80 и 100 Гц и что минимальное ослабление для компонентов вне полосы составляет 80 дБ. Кроме того, максимальное искажение для интересующих компонентов составляет 0,05 дБ (половина неравномерности в полосе пропускания от пика до пика). Фильтр equiripple, который соответствует этим спецификациям, может быть легко получен fdesign интерфейс.

eqrDecim = design(fdDecim,'equiripple', 'SystemObject', true);
measure(eqrDecim)

ans = 

Sample Rate      : 800 Hz     
Passband Edge    : 80 Hz      
3-dB Point       : 85.621 Hz  
6-dB Point       : 87.8492 Hz 
Stopband Edge    : 100 Hz     
Passband Ripple  : 0.092414 dB
Stopband Atten.  : 80.3135 dB 
Transition Width : 20 Hz

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

fvtool(eqrDecim)

Фильтры Nyquist и согласованность интерполяции

Цифровой фильтр свертки$h$ называется L-м фильтром Найквиста, если он периодически исчезает с каждой$L$ выборки, кроме центрального индекса. Другими словами, дискретизация в$h$ множителе$L$ приводит к импульсу:

$$h[kL] = \delta_k$$

Идеальным $L$lowpass -го полосы, $h[m] = sinc(\frac{m}{L})$например, является -$L$й фильтр Найквиста. Другой пример - треугольное окно.

L=3;
t = linspace(-3*L,3*L,1024);
n = (-3*L:3*L);
hLP = @(t) sinc(t/L);
hTri = @(t) (1-abs(t/L)).*(abs(t/L)<=1);

plot_nyquist_filter(t,n,hLP,hTri,L);

Функция designMultirateFIR приводит к фильтрам Nyquist, поскольку он основан на взвешенных и усеченных версиях идеальных фильтров Nyquist.

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

Согласованность интерполяции

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

Интерполяционная согласованность сохранена в фильтре Найквиста, поскольку коэффициенты равны нулю на каждых L выборках (кроме в центре). Доказательство просто. Предположим, что$x_L$ это увеличенная версия ($x$с нулями, вставленными между выборками), так что, $x_L[Ln]=x[n]$и это интерполированный$y=h*x_L$ сигнал. Дискретизируйте$y$ равномерно и получите следующее уравнение.

$$ y[nL] = (h*x_L)[Ln] = \sum_k h[Ln-k] x_L[k] = \sum_m \underbrace{h[Ln-Lm]}_{\delta[L(n-m)]} x_L[Lm] = x_L[Ln] = x[n]$$

Давайте рассмотрим эффект использования фильтра Nyquist для интерполяции. The designMultirateFIR функция создает фильтры Nyquist. Как видно на изображении ниже, входные значения совпадают с интерполированными значениями.

% Generate input
n = (0:20)';
xInput = (n<=10).*cos(pi*0.05*n).*(-1).^n;

L = 4;
hNyq = designMultirateFIR(L,1);
firNyq = dsp.FIRInterpolator(L,hNyq);
xIntrNyq = firNyq(xInput);
release(firNyq);
plot_shape_and_response(hNyq,xIntrNyq,xInput,L,num2str(L)+"-Nyuist");

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

hNotNyq = firpm(length(hNyq)-1,[0 1/L 1.5/L 1],[1 1 0 0]);
hNotNyq = hNotNyq/max(hNotNyq); % Adjust gain
firIntrNotNyq = dsp.FIRInterpolator(L,hNotNyq);
xIntrNotNyq= firIntrNotNyq(xInput);
release(firIntrNotNyq);

plot_shape_and_response(hNotNyq,xIntrNotNyq,xInput,L,"equiripple, not Nyquist");

См. также

Функции

Объекты

Похожие темы

Документация DSP System Toolbox

Поддержка