Exponenta Banner
exponenta event banner

Проектирование прореживателей и интерполяторов

Открыть сценарий

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

Роль фильтрации нижних частот в преобразовании скорости

Преобразование скорости - это процесс изменения скорости дискретного сигнала для получения нового дискретного представления основного непрерывного сигнала. Процесс включает равномерную понижающую и повышающую дискретизацию. Равномерное понижение дискретизации на скорость 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

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

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

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

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

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

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

Отфильтрованное преобразование скорости: прореживатели, интерполяторы и рациональные преобразователи скорости

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

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

Прореживание относится к фильтрации LTI с последующим равномерным понижением дискретизации. Дециматор КИХ может быть реализован следующим образом.

  1. Проектирование сглаживающего фильтра нижних частот h

  2. Фильтрация входного сигнала через h

  3. Понизить отфильтрованную последовательность на коэффициент М

% 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)

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

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

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

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

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

Проектирование фильтра преобразования скорости с использованием designMultirateFIR

Функция designMultirateFIR(L,M) автоматически находит соответствующую частоту масштабирования и отсечки для данного коэффициента преобразования скорости. $L/M$Использовать коэффициенты FIR, возвращенные 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, есть несколько ключевых различий. В интерполированном сигнале

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

  • Сигнал имеет задержку в половину длины FIR 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]);

Настройка параметров проектирования КИХ нижних частот

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

Регулировка длины КИХ

Длина КИХ может управляться через 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. 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)

Дизайн Equiripple

Функция designMultirateFIR использует оконную конструкцию нижних частот FIR. Также могут быть применены другие способы проектирования нижних частот, такие как 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$нижним полосой -го диапазона, $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 эффективны для реализации, поскольку 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]$$

Рассмотрим эффект использования фильтра Найквиста для интерполяции. 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");

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

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

Поддержка