В этом примере показано, как спроектировать пиковые и узкополосные фильтры. Фильтры, которые пиковые или узкополосные на определенной частоте, полезны, чтобы сохранить или исключить конкретную частотную составляющую сигнала. Расчётные параметры для фильтра - это частота, на которой требуется пик или надрез, и либо 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');
Коэффициент качества или Q-коэффициент фильтра является мерой того, насколько хорошо требуемая частота изолирована от других частот. Для фиксированного порядка фильтра более высокий Q-коэффициент достигается путем толкания полюсов ближе к нулям.
notchspec.Q = 100; notchfilt1 = design(notchspec,'SystemObject',true); fvt= fvtool(notchfilt, notchfilt1, 'Color','white'); legend(fvt,'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');
Поскольку пока можно только толкать полюса и оставаться стабильным, для улучшения приближения кирпичной стенки фильтра необходимо увеличить порядок фильтра.
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
Коэффициенты изменяющихся во времени фильтров должны изменяться с течением времени из-за изменений во время выполнения расчётных параметров (например, центральная частота для узкополосного фильтра). Для каждого изменения параметров расчётных параметров коэффициентов b и a должны быть пересчитаны и биквад. Для SOSMatrix задано новое значение. Это может быть дорогим в вычислительном отношении. Для этого типа приложения dsp. CoupledAllpassFilter может предложить более удобные структуры фильтра. Преимущества включают - Внутренняя стабильность - Коэффициенты, разделенные по отношению к расчётным параметрам
Создайте связанный решетчатый фильтр allpass, эквивалентный biquad
[k1, k2] = tf2cl(b, a) %#okk1 =
[]
k2 = 2×1
-0.7071
0.6682
apnotch = dsp.CoupledAllpassFilter('Lattice', k1, k2); fvtool(apnotch, 'Fs', Fs)
Одним из преимуществ этой структуры, основанной на альпасе, является то, что новые коэффициенты разделены относительно расчётных параметров 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)%#okk1 =
[]
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)%#okk1 =
[]
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
