Проект пиковых и узкополосных фильтров

Открыть Live Script

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

Узкие фильтры

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

F0 = 60;   % interference is at 60 Hz
Fs = 3000; % sampling frequency is 3000 Hz
notchspec = fdesign.notch('N,F0,Q',2,F0,10,Fs);
notchfilt = design(notchspec,'SystemObject',true);
fvtool(notchfilt,'Color','white');

Figure Filter Visualization Tool - Magnitude Response (dB) contains an axes and other objects of type uitoolbar, uimenu. The axes with title Magnitude Response (dB) contains an object of type line.

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

notchspec.Q = 100;
notchfilt1 = design(notchspec,'SystemObject',true);
fvt= fvtool(notchfilt, notchfilt1, 'Color','white');
legend(fvt,'Q = 10','Q = 100');

Figure Filter Visualization Tool - Magnitude Response (dB) contains an axes and other objects of type uitoolbar, uimenu. The axes with title Magnitude Response (dB) contains 2 objects of type line. These objects represent Q = 10, Q = 100.

Эквивалентный способ определения фактора качества для определения 3-dB полосы пропускания, BW. Они связаны Q = F0/BW. Определение полосы пропускания может быть более удобным способом достижения точно желаемой формы для проектируемого фильтра.

notchspec = fdesign.notch('N,F0,BW',2,60,5,3000);
notchfilt2 = design(notchspec,'SystemObject',true);
fvt= fvtool(notchfilt, notchfilt1, notchfilt2, 'Color','white');
legend(fvt,'Q = 10','Q = 100','BW = 5 Hz');

Figure Filter Visualization Tool - Magnitude Response (dB) contains an axes and other objects of type uitoolbar, uimenu. The axes with title Magnitude Response (dB) contains 3 objects of type line. These objects represent Q = 10, Q = 100, BW = 5 Hz.

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

notchspec = fdesign.notch('N,F0,Q',2,.4,100);
notchfilt = design(notchspec,'SystemObject',true);
notchspec.FilterOrder = 6;
notchfilt1 = design(notchspec,'SystemObject',true);
fvt= fvtool(notchfilt, notchfilt1, 'Color','white');
legend(fvt,'2nd-Order Filter','6th-Order Filter');

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

N = 8; F0 = 0.4; BW = 0.1;
notchspec = fdesign.notch('N,F0,BW',N,F0,BW);
notchfilt = design(notchspec,'SystemObject',true);
notchspec1 = fdesign.notch('N,F0,BW,Ap,Ast',N,F0,BW,0.5,60);
notchfilt1 = design(notchspec1,'SystemObject',true);
fvt= fvtool(notchfilt, notchfilt1, 'Color','white');
legend(fvt,'Maximally Flat 8th-Order Filter',...
    '8th-Order Filter With Passband/Stopband Ripples', ...
    'Location','NorthWest');
axis([0 1 -90 0.5]);

Пиковые фильтры

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

N = 6; F0 = 0.7; BW = 0.001;
peakspec = fdesign.peak('N,F0,BW',N,F0,BW);
peakfilt = design(peakspec,'SystemObject',true);
peakspec1 = fdesign.peak('N,F0,BW,Ast',N,F0,BW,80);
peakfilt1 = design(peakspec1,'SystemObject',true);
fvt= fvtool(peakfilt, peakfilt1, 'Color','white');
legend(fvt,'Maximally Flat 6th-Order Filter',...
    '6th-Order Filter With 80 dB Stopband Attenuation','Location','North');

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

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

Полезной начальной точкой является статический фильтр, вызываемый из динамической (потоковой) симуляции. В этом случае непосредственно создается узкополосный фильтр 2-го порядка и его коэффициенты вычисляются с помощью iirnotch. Эти расчётные параметры являются центральной частотой 1 кГц и полосой пропускания 50 Гц на -3 дБ с частотой дискретизации 8 кГц.

Fs = 8e3;           % 8 kHz sampling frequency
F0 = 1e3/(Fs/2);    % Notch at 2 kHz
BW = 500/(Fs/2);    % 500 Hz bandwidth
[b, a] = iirnotch(F0, BW)

b = 1×3

    0.8341   -1.1796    0.8341

a = 1×3

    1.0000   -1.1796    0.6682

biquad = dsp.BiquadFilter('SOSMatrix', [b, a]);

scope = dsp.SpectrumAnalyzer('PlotAsTwoSidedSpectrum', false, ...
    'SampleRate', Fs, ...
    'SpectralAverages', 16);
samplesPerFrame = 256;

nFrames = 4096;
for k = 1:nFrames
   x = randn(samplesPerFrame, 1);
   y = biquad(x);
   scope(y);
end

Figure Spectrum Analyzer contains an axes and other objects of type uiflowcontainer, uimenu, uitoolbar. The axes contains an object of type line. This object represents Channel 1.

Коэффициенты изменяющихся во времени фильтров должны изменяться с течением времени из-за изменений во время выполнения расчётных параметров (например, центральная частота для узкополосного фильтра). Для каждого изменения параметров расчётных параметров коэффициентов b и a должны быть пересчитаны и биквад. Для SOSMatrix задано новое значение. Это может быть дорогим в вычислительном отношении. Для этого типа приложения dsp. CoupledAllpassFilter может предложить более удобные структуры фильтра. Преимущества включают - Внутренняя стабильность - Коэффициенты, разделенные по отношению к расчётным параметрам

Создайте связанный решетчатый фильтр allpass, эквивалентный biquad

[k1, k2] = tf2cl(b, a) %#ok

k1 =

     []

apnotch = dsp.CoupledAllpassFilter('Lattice', k1, k2);
fvtool(apnotch, 'Fs', Fs)

Figure Filter Visualization Tool - Magnitude Response (dB) contains an axes and other objects of type uitoolbar, uimenu. The axes with title Magnitude Response (dB) contains an object of type line.

Одним из преимуществ этой структуры, основанной на альпасе, является то, что новые коэффициенты разделены относительно расчётных параметров F0 и BW. Так, например, изменение F0 на 3 кГц даст

F0 = 3e3/(Fs/2);
[b, a] = iirnotch(F0, BW)

b = 1×3

    0.8341    1.1796    0.8341

a = 1×3

    1.0000    1.1796    0.6682

[k1, k2] = tf2cl(b, a)%#ok

k1 =

     []

ПРИМЕЧАНИЕ: в то время как a и b оба изменились, в форме allpass решетки изменение проекта затронуло только k2 (1).

Теперь измените пропускную способность на 1 кГц

BW = 1e3/(Fs/2);
[b, a] = iirnotch(F0, BW)

b = 1×3

    0.7071    1.0000    0.7071

a = 1×3

    1.0000    1.0000    0.4142

[k1, k2] = tf2cl(b, a)%#ok

k1 =

     []

Изменение проекта теперь затронуло только k2 (2).

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

Следующий принцип применяется к изменению расчётных параметров во время динамической симуляции, включая живую визуализацию эффектов на оценку передаточной функции фильтра. В реальных приложениях обычно можно идентифицировать отдельные выражения, которые связывают каждые расчётные параметры с соответствующим коэффициентом allpass решетки. Использование функций вместо функций всей фильтрации, таких как iirnotch и tf2cl в основном цикле симуляции, повысило бы эффективность.

% Notch filter parameters - how they vary over time
Fs = 8e3;       % 8 kHz sampling frequency
f0 = 1e3*[0.5, 1.5, 3, 2]/(Fs/2);   % in [0, 1] range
bw = 500/(Fs/2) * ones(1,4);        % in [0, 1] range

myChangingParams = struct('F0', num2cell(f0), 'BW', num2cell(bw));
paramsChangeTimes = [0, 7, 15, 20]; % in seconds

% Simulation time management
nSamplesPerFrame = 256;
tEnd = 30;
nSamples = ceil(tEnd * Fs);
nFrames = floor(nSamples / nSamplesPerFrame);

% Object creation
apnotch = dsp.CoupledAllpassFilter('Lattice');
scope = dsp.SpectrumAnalyzer('PlotAsTwoSidedSpectrum', false, ...
    'SampleRate', Fs, ...
    'SpectralAverages', 4, ...
    'RBWSource', 'Auto');
paramtbl = ParameterTimeTable('Time', paramsChangeTimes, ...
    'Values', myChangingParams, ...
    'SampleRate', Fs/nSamplesPerFrame);

% Actual simulation loop
for frameIdx = 1:nFrames
    % Get current F0 and BW
    [params, update] = paramtbl();
    
    if(update)
        % Recompute filter coefficients if parameters changed
        [b, a] = iirnotch(params.F0, params.BW);
        [k1, k2] = tf2cl(b, a);
        % Set filter coefficients to new values
        apnotch.LatticeCoefficients1{1} = k1;
        apnotch.LatticeCoefficients2{1} = k2;
    end
    % Generate vector of white noise samples
    x = randn(nSamplesPerFrame, 1);
    % Filter noise
    y = apnotch(x);
    % Visualize spectrum
    scope(y);
end

Figure Spectrum Analyzer contains an axes and other objects of type uiflowcontainer, uimenu, uitoolbar. The axes contains an object of type line. This object represents Channel 1.

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

Проект пиковых и узкополосных фильтров

Узкие фильтры

Пиковые фильтры

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

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

Поддержка