Проект достигания максимума и отметки фильтров

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

Отметьте фильтры

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

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 дБ, 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');

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 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]);

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 Maximally Flat 8th-Order Filter, 8th-Order Filter With Passband/Stopband Ripples.

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

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

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');

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 Maximally Flat 6th-Order Filter, 6th-Order Filter With 80 dB Stopband Attenuation.

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

Используя изменяющиеся во времени фильтры требует изменения коэффициентов фильтра, в то время как симуляция запускается. Чтобы дополнить автоматический рабочий процесс создания фильтра на основе объектов fdesign, DSP System Toolbox предусматривает другие возможности, включая функции, чтобы вычислить коэффициенты фильтра непосредственно, e.g. 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.

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

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

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

     []
k2 = 2×1

   -0.7071
    0.6682

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.

Одно преимущество этой находящейся в allpass структуры - то, что новые коэффициенты разъединяются относительно расчетных параметров 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 =

     []
k2 = 2×1

    0.7071
    0.6682

ПРИМЕЧАНИЕ: в то время как 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×1

    0.7071
    0.4142

Конструктивное изменение теперь только влияло на 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.