Filters — синтезирование фильтров и фильтрация
В DSP.jl различаются понятия коэффициентов фильтров и фильтров с сохранением состояния. Объекты коэффициентов фильтра определяют характеристику фильтра в одной из стандартных форм. Объекты фильтров с сохранением состояния содержат состояние фильтра вместе с его коэффициентами в реализуемой форме (PolynomialRatio
, Biquad
или SecondOrderSections
). При вызове для объекта коэффициента фильтра filt
не сохраняет состояние.
Объекты коэффициентов фильтра
DSP.jl поддерживают стандартные представления фильтров. Коэффициенты фильтра можно преобразовывать из одного типа в другой с помощью convert
.
#
DSP.Filters.ZeroPoleGain
— Type
ZeroPoleGain(z, p, k)
Представление фильтра в виде нулей z
, полюсов p
и коэффициента усиления k
:
#
DSP.Filters.PolynomialRatio
— Type
PolynomialRatio(b, a)
Представление фильтра в виде коэффициентов числителя b
и знаменателя a
передаточной функции:
или, что то же самое:
b
и a
можно указывать в виде объектов Polynomial
или векторов, упорядоченных в порядке убывания степени.
#
DSP.Filters.Biquad
— Type
Biquad(b0, b1, b2, a1, a2)
Представление фильтра в виде передаточной функции одного сечения второго порядка, заданного следующим уравнением:
или, что то же самое:
#
DSP.Filters.SecondOrderSections
— Type
SecondOrderSections(biquads, gain)
Представление фильтра в виде передаточной функции каскада сечений второго порядка и коэффициента усиления. Аргумент biquads
должен задаваться в виде вектора Biquads
.
Объекты коэффициентов фильтра поддерживают следующие арифметические операции: инверсия (inv
), умножение (*
) для последовательного соединения и возведение в целую степень (^
) для многократного умножения на себя. Например:
julia> H = PolynomialRatio([1.0], [1.0, 0.3])
PolynomialRatio{:z, Float64}(Polynomials.LaurentPolynomial(1.0), Polynomials.LaurentPolynomial(0.3*z⁻¹ + 1.0))
julia> inv(H)
PolynomialRatio{:z, Float64}(Polynomials.LaurentPolynomial(0.3*z⁻¹ + 1.0), Polynomials.LaurentPolynomial(1.0))
julia> H * H
PolynomialRatio{:z, Float64}(Polynomials.LaurentPolynomial(1.0), Polynomials.LaurentPolynomial(0.09*z⁻² + 0.6*z⁻¹ + 1.0))
julia> H^2
PolynomialRatio{:z, Float64}(Polynomials.LaurentPolynomial(1.0), Polynomials.LaurentPolynomial(0.09*z⁻² + 0.6*z⁻¹ + 1.0))
julia> H^-2
PolynomialRatio{:z, Float64}(Polynomials.LaurentPolynomial(0.09*z⁻² + 0.6*z⁻¹ + 1.0), Polynomials.LaurentPolynomial(1.0))
Объекты фильтров с сохранением состояния
#
DSP.Filters.DF2TFilter
— Type
DF2TFilter(coef[, si])
Создает фильтр транспонированной прямой формы II с сохранением состояния с коэффициентами coef
. si
— необязательный массив, представляющий начальное состояние фильтра (по умолчанию заполнен нулями). Если f
— это объект PolynomialRatio
, Biquad
или SecondOrderSections
, фильтрация применяется напрямую. Если f
— это объект ZeroPoleGain
, сначала он преобразуется в объект SecondOrderSections
.
Тип FIRFilter
в DSP.jl сохраняет состояние между вызовами filt
, что позволяет фильтровать сигнал неопределенной длины фрагментами, помещающимися в оперативной памяти. Объект FIRFilter
содержит лишь состояние фильтра и FIRKernel
. Есть пять разных видов FIRKernel
для преобразования с постоянной частотой дискретизации, повышением частоты дискретизации, понижением частоты дискретизации, рациональной передискретизацией и произвольной частотой дискретизации. Указывать тип ядра не нужно. Конструктор FIRFilter
выбирает подходящее ядро в зависимости от входных параметров.
#
DSP.Filters.FIRFilter
— Type
FIRFilter(h[, ratio])
Создает объект FIRFilter с сохранением состояния из вектора отводов фильтра h
. ratio
— необязательное рациональное целое число, определяющее отношение частоты дискретизации между входом и выходом (например, 147//160
для преобразования аудиозаписи из 48 кГц в 44,1 кГц).
FIRFilter(h, rate[, Nϕ])
Возвращает многофазный объект FIRFilter из вектора отводов фильтра h
. rate
— это число с плавающей запятой, определяющее отношение частоты дискретизации между входом и выходом . Nϕ
— это необязательный параметр, который определяет количество фаз, создаваемых из h
. Значение Nϕ
по умолчанию — 32.
Применение фильтров
#
DSP.filt
— Function
filt(f, x[, si])
Применяет фильтр или коэффициенты фильтра f
по первому измерению массива x
. Если f
— объект коэффициентов фильтра, то si
— это необязательный массив, представляющий начальное состояние фильтра (по умолчанию заполнен нулями). Если f
— это объект PolynomialRatio
, Biquad
или SecondOrderSections
, фильтрация применяется напрямую. Если f
— это объект ZeroPoleGain
, сначала он преобразуется в объект SecondOrderSections
. Если f
— это вектор, он интерпретируется как КИХ-фильтр, и в зависимости от длины данных и фильтра выбирается наивный алгоритм или алгоритм на основе БПФ.
filt(b, a, x, [si])
Применяет фильтр, описываемый векторами a
и b
, к вектору x
с необязательным вектором начального состояния фильтра si
(по умолчанию заполнен нулями).
#
DSP.filt!
— Function
filt!(out, f, x[, si])
То же, что и filt()
, но записывает результат в аргумент out
. Выходной массив out
может не являться псевдонимом x
, то есть фильтрация может выполняться не на месте.
filt!(out, b, a, x, [si])
То же, что и filt
, но записывает результат в аргумент out
, который может являться псевдонимом входного x
для его изменения на месте.
#
DSP.Filters.filtfilt
— Function
filtfilt(coef, x)
Фильтрует x
в прямом и обратном направлениях с использованием коэффициентов фильтра coef
. Начальное состояние фильтра рассчитывается так, чтобы его отклик на ступенчатую функцию имел устойчивое состояние. Перед фильтрацией данные на обоих концах экстраполируются с использованием нечетно-симметрического расширения длиной 3*(max(length(b), length(a))-1)
.
Так как filtfilt
применяет фильтр дважды, действительный порядок фильтра вдвое больше порядка coef
. Итоговый сигнал имеет нулевое фазовое искажение.
#
DSP.Filters.fftfilt
— Function
fftfilt(h, x)
Применяет отводы КИХ-фильтра h
по первому измерению массива x
с использованием алгоритма перекрытий с сохранением на основе БПФ.
#
DSP.Filters.fftfilt!
— Function
fftfilt!(out, h, x)
Действует аналогично fftfilt
, но записывает результат в массив out.
#
DSP.Filters.tdfilt
— Function
tdfilt(h, x)
Применяет фильтр или коэффициенты фильтра h
по первому измерению массива x
с использованием наивного алгоритма временной области.
#
DSP.Filters.tdfilt!
— Function
tdfilt!(out, h, x)
Действует аналогично tdfilt
, но записывает результат в массив out
. Выходной массив out
может не являться псевдонимом x
, то есть фильтрация может выполняться не на месте.
#
DSP.Filters.resample
— Function
resample(x, rate[, coef])
Повторно дискретизирует x
с рациональной или произвольной частотой rate
. coef
— это необязательный вектор отводов КИХ-фильтра. Если аргумент coef
не указан, отводы вычисляются с использованием окна Кайзера.
На внутреннем уровне функция resample
использует многофазный объект FIRFilter
, но выполняет дополнительные операции для упрощения повторной дискретизации сигнала. Она компенсирует запаздывание (линейное нарастание) FIRFilter
и добавляет нули к x
. В результате, когда входной и выходной сигналы отрисовываются на графике с наложением, они очень хорошо коррелируют, но один из сигналов будет иметь больше отсчетов, чем другой.
resample(x::AbstractArray, rate::Real, h::Vector = resample_filter(rate); dims)
Повторно дискретизирует массив x
по измерению dims
.
Синтезирование фильтров
Большинство аналоговых и цифровых фильтров создаются путем композиции типов характеристик, которые определяют частотную характеристику фильтра, с методами синтеза, которые определяют способ построения фильтра. Возможные типы характеристик: Lowpass
, Highpass
, Bandpass
и Bandstop
. Они включают в себя границы полос. Возможные методы синтеза: Butterworth
, Chebyshev1
, Chebyshev2
, Elliptic
и FIRWindow
. Они включают в себя все необходимые параметры метода, влияющие на форму характеристики, например порядок фильтра, неравномерность характеристики и затухание. В buttord
, cheb1ord
, cheb2ord
и ellipord
доступны методы оценки порядка фильтра, если известны частоты излома для различных типов БИХ-фильтров. remezord
можно использовать для начальной оценки порядка КИХ-фильтра.
#
DSP.Filters.analogfilter
— Function
analogfilter(responsetype, designmethod)
Создает аналоговый фильтр. Возможные типы характеристик и фильтров см. ниже.
#
DSP.Filters.digitalfilter
— Function
digitalfilter(responsetype, designmethod)
Создает цифровой фильтр. Возможные типы характеристик и фильтров см. ниже.
Для некоторых фильтров метод синтеза является более общим или изначально предполагает определенный тип характеристики. К таким прямым методам синтеза относятся remez
для синтезирования КИХ-фильтров с равномерной пульсацией АЧХ любых типов и iirnotch
для синтезирования «биквадратного» запирающего БИХ-фильтра 2-го порядка.
Типы характеристик фильтров
#
DSP.Filters.Lowpass
— Type
Lowpass(Wn[; fs])
Низкочастотный фильтр с частотой среза Wn
. Если аргумент fs
не указан, Wn
интерпретируется как нормализованная частота в полуциклах на отсчет.
#
DSP.Filters.Highpass
— Type
Highpass(Wn[; fs])
Высокочастотный фильтр с частотой среза Wn
. Если аргумент fs
не указан, Wn
интерпретируется как нормализованная частота в полуциклах на отсчет.
#
DSP.Filters.Bandpass
— Type
Bandpass(Wn1, Wn2[; fs])
Полосовой пропускающий фильтр с нормализованной полосой пропускания (Wn1
, Wn2
). Если аргумент fs
не указан, Wn1
и Wn2
интерпретируются как нормализованные частоты в полуциклах на отсчет.
#
DSP.Filters.Bandstop
— Type
Bandstop(Wn1, Wn2[; fs])
Полосовой задерживающий фильтр с нормализованной полосой задержки (Wn1
, Wn2
). Если аргумент fs
не указан, Wn1
и Wn2
интерпретируются как нормализованные частоты в полуциклах на отсчет.
Методы синтеза фильтров
Методы синтеза БИХ-фильтров
#
DSP.Filters.Chebyshev1
— Function
Chebyshev1(n, ripple)
n
-полюсный фильтр Чебышева первого рода с неравномерностью характеристики ripple
дБ в полосе пропускания.
#
DSP.Filters.Chebyshev2
— Function
Chebyshev2(n, ripple)
n
-полюсный фильтр Чебышева второго рода с неравномерностью характеристики ripple
дБ в полосе задержки.
#
DSP.Filters.Elliptic
— Function
Elliptic(n, rp, rs)
n
-полюсный эллиптический фильтр (фильтр Кауэра) с неравномерностью характеристики rp
дБ в полосе пропускания и затуханием rs
дБ в полосе задержки.
Методы оценки порядка фильтров
Методы оценки порядка БИХ-фильтров
#
DSP.Filters.buttord
— Function
(N, ωn) = buttord(Wp::Tuple{Real,Real}, Ws::Tuple{Real,Real}, Rp::Real, Rs::Real; domain=:z)
Оценка порядка для фильтров Баттерворта пропускающего и задерживающего типов. Wp
и Ws
— это двухэлементные границы частот полос пропускания и задержки с неравномерностью характеристики в полосе пропускания не более Rp
дБ и затуханием в полосе задержки не более Rs
дБ. В зависимости от порядка следования границ полосы пропускания и полосы задержки выводится тип фильтра: пропускающий или задерживающий. N
— это целое число, указывающее на нижний оценочный порядок фильтра, а ωn
указывает частоту среза или частоту «-3 дБ». Если указана область :s
, частоты полосы пропускания и полосы задержки интерпретируются как выраженные в радианах в секунду при заданном порядке и собственных частотах аналогового фильтра. Область по умолчанию — :z
. При этом входные частоты интерпретируются как нормализованные в диапазоне от 0 до 1, где 1 соответствует π радиан/отсчет.
(N, ωn) = buttord(Wp::Real, Ws::Real, Rp::Real, Rs::Real; domain=:z)
Порядок ФНЧ/ФВЧ-фильтра Баттерворта и частотная аппроксимация -3 дБ. Wp
и Ws
— это частоты полосы пропускания и полосы задержки, а Rp и Rs — затухание неравномерности в полосах пропускания и задержки в дБ. Если полоса пропускания больше полосы задержки, тип фильтра выводится как предназначенный для оценки порядка высокочастотного фильтра. N
определяет наименьший возможный целочисленный порядок фильтра, а ωn
— это частота среза или частота «-3 дБ». Если указана область :s
, границы полосы пропускания и полосы задержки интерпретируются как выраженные в радианах в секунду при заданном порядке и полученной собственной частоте аналогового фильтра. Область по умолчанию — :z
. При этом входные частоты интерпретируются как нормализованные в диапазоне от 0 до 1, где 1 соответствует π радиан/отсчет.
#
DSP.Filters.cheb1ord
— Function
(N, ωn) = cheb1ord(Wp::Real, Ws::Real, Rp::Real, Rs::Real; domain::Symbol=:z)
Оценка целочисленного порядка и собственной частоты для фильтров Чебышева первого рода. Wp
и Ws
означают границы частот полосы пропускания и полосы задержки, а Rp
и Rs
— максимальные потери в полосе пропускания и минимальное затухание в полосе задержки (в дБ). В зависимости от порядка следования границ полосы пропускания и полосы задержки выводится тип фильтра: низкочастотный или высокочастотный. N
означает наименьший целочисленный порядок фильтра, при котором выполняются заданные требования, а ωn
содержит собственную частоту фильтра (в данном случае просто границу полосы пропускания). Если указана область :s
, границы полосы пропускания и полосы задержки интерпретируются как выраженные в радианах в секунду при заданном порядке и полученной собственной частоте аналогового фильтра. Область по умолчанию — :z
. При этом входные частоты интерпретируются как нормализованные в диапазоне от 0 до 1, где 1 соответствует π радиан/отсчет.
(N, ωn) = cheb1ord(Wp::Tuple{Real, Real}, Ws::Tuple{Real, Real}, Rp::Real, Rs::Real; domain::Symbol=:z)
Оценка целочисленного порядка и собственной частоты для фильтров Чебышева первого рода. Wp
и Ws
— это двухэлементные границы частот для пропускающего и задерживающего вариантов, а Rp
и Rs
представляют максимальные потери вследствие неравномерности характеристики в полосе пропускания и минимальное затухание вследствие неравномерности характеристики в полосе задержки в дБ. В зависимости от порядка следования границ полосы пропускания и полосы задержки выводится тип фильтра: пропускающий или задерживающий. N
— это целое число, указывающее на нижний оценочный порядок фильтра, а ωn
указывает частоту среза или частоту «-3 дБ». Если указана область :s
, частоты полосы пропускания и полосы задержки интерпретируются как выраженные в радианах в секунду при заданном порядке и собственных частотах аналогового фильтра. Область по умолчанию — :z
. При этом входные частоты интерпретируются как нормализованные в диапазоне от 0 до 1, где 1 соответствует π радиан/отсчет.
#
DSP.Filters.cheb2ord
— Function
(N, ωn) = cheb2ord(Wp::Real, Ws::Real, Rp::Real, Rs::Real; domain::Symbol=:z)
Оценка целочисленного порядка и собственной частоты для фильтров Чебышева второго рода (инверсных). Wp
и Ws
— это границы частот для пропускающего и задерживающего вариантов, а Rp
и Rs
представляют максимальные потери вследствие неравномерности характеристики в полосе пропускания и минимальное затухание вследствие неравномерности характеристики в полосе задержки в дБ. В зависимости от порядка следования границ полосы пропускания и полосы задержки выводится тип фильтра: низкочастотный или высокочастотный. N
— это целое число, указывающее на нижний порядок фильтра, а ωn
указывает частоту среза «-3 дБ». Если указана область :s
, частоты полосы пропускания и полосы задержки интерпретируются как выраженные в радианах в секунду при заданном порядке и собственных частотах аналогового фильтра. Область по умолчанию — :z
. При этом входные частоты интерпретируются как нормализованные в диапазоне от 0 до 1, где 1 соответствует π радиан/отсчет.
(N, ωn) = cheb2ord(Wp::Tuple{Real, Real}, Ws::Tuple{Real, Real}, Rp::Real, Rs::Real; domain::Symbol=:z)
Оценка целочисленного порядка и собственной частоты для фильтров Чебышева второго рода (инверсных). Wp
и Ws
— это двухэлементные границы частот для пропускающего и задерживающего вариантов, а Rp
и Rs
представляют максимальные потери вследствие неравномерности характеристики в полосе пропускания и минимальное затухание вследствие неравномерности характеристики в полосе задержки в дБ. В зависимости от порядка следования границ полосы пропускания и полосы задержки выводится тип фильтра: пропускающий или задерживающий. N
— это целое число, указывающее на нижний оценочный порядок фильтра, а ωn
указывает частоту среза или частоту «-3 дБ». Если указана область :s
, частоты полосы пропускания и полосы задержки интерпретируются как выраженные в радианах в секунду при заданном порядке и собственных частотах аналогового фильтра. Область по умолчанию — :z
. При этом входные частоты интерпретируются как нормализованные в диапазоне от 0 до 1, где 1 соответствует π радиан/отсчет.
#
DSP.Filters.ellipord
— Function
(N, ωn) = ellipord(Wp::Real, Ws::Real, Rp::Real, Rs::Real; domain::Symbol=:z)
Оценка целочисленного порядка и собственной частоты для эллиптических фильтров (фильтров Кауэра). Wp
и Ws
означают границы частот полосы пропускания и полосы задержки, а Rp
и Rs
— максимальные потери в полосе пропускания и минимальное затухание в полосе задержки (в дБ). В зависимости от порядка следования границ полосы пропускания и полосы задержки выводится тип фильтра: низкочастотный или высокочастотный. N
означает наименьший целочисленный порядок фильтра, при котором выполняются заданные требования, а ωn
содержит собственную частоту фильтра (в данном случае просто границу полосы пропускания). Если указана область :s
, границы полосы пропускания и полосы задержки интерпретируются как выраженные в радианах в секунду при заданном порядке и полученной собственной частоте аналогового фильтра. Область по умолчанию — :z
. При этом входные частоты интерпретируются как нормализованные в диапазоне от 0 до 1, где 1 соответствует π радиан/отсчет.
(N, ωn) = ellipord(Wp::Tuple{Real, Real}, Ws::Tuple{Real, Real}, Rp::Real, Rs::Real; domain::Symbol=:z)
Оценка целочисленного порядка и собственной частоты для эллиптических фильтров (фильтров Кауэра). Wp
и Ws
— это двухэлементные границы частот для пропускающего и задерживающего вариантов, а Rp
и Rs
представляют максимальные потери вследствие неравномерности характеристики в полосе пропускания и минимальное затухание вследствие неравномерности характеристики в полосе задержки в дБ. В зависимости от порядка следования границ полосы пропускания и полосы задержки выводится тип фильтра: пропускающий или задерживающий. N
— это целое число, указывающее на нижний оценочный порядок фильтра, а ωn
указывает частоту среза или частоту «-3 дБ». Если указана область :s
, частоты полосы пропускания и полосы задержки интерпретируются как выраженные в радианах в секунду при заданном порядке и собственных частотах аналогового фильтра. Область по умолчанию — :z
. При этом входные частоты интерпретируются как нормализованные в диапазоне от 0 до 1, где 1 соответствует π радиан/отсчет.
Методы оценки порядка КИХ-фильтров
#
DSP.Filters.remezord
— Function
N = remezord(Wp, Ws, Rp, Rs)
Оценка порядка для случаев низкочастотных цифровых фильтров на основе уравнений и коэффициентов из работы [1]. Оригинальное уравнение возвращало минимальную длину фильтра, в то время как данная реализация возвращает порядок (N=L-1). Wp
и Ws
— это нормализованные частоты полосы пропускания и полосы задержки, Rp
— затухание в полосе пропускания, а Rs
— неравномерность характеристики в полосе задержки (линейная).
ПРИМЕЧАНИЕ. Значение N представляет собой первое приближение. Если оценка занижена, порядок следует увеличить, пока не будут выполнены требования к синтезу.
Методы синтеза КИХ-фильтров
#
DSP.Filters.FIRWindow
— Type
FIRWindow(window; scale=true)
Синтезирование КИХ-фильтра с использованием окна window
, вектора, длина которого соответствует количеству отводов в итоговом фильтре.
Если аргумент scale
имеет значение true
(по умолчанию), синтезируемый КИХ-фильтр масштабируется так, чтобы было верно следующее.
FIRWindow(; transitionwidth, attenuation=60, scale=true)
Синтезирование КИХ-фильтра с окном Кайзера. Требуемое число отводов рассчитывается на основе transitionwidth
(в полуциклах на отсчет) и полосы задержки attenuation
(в дБ). attenuation
по умолчанию имеет значение 60 дБ.
Прямые методы синтеза фильтров
#
DSP.Filters.remez
— Function
remez(numtaps::Integer, band_defs;
Hz::Real=1.0,
neg::Bool=false,
maxiter::Integer=25,
grid_density::Integer=16)
Вычисляет оптимальный по минимаксному критерию фильтр с использованием алгоритма замены Ремеза [2] [3].
Это упрощенный API, который принимает всего два обязательных аргумента (numtaps, band_defs). Имеется также совместимая с SciPy версия с тремя аргументами (numtaps, bands, desired).
Вычисляет коэффициенты для фильтра с конечной импульсной характеристикой (КИХ), передаточная функция которого минимизирует максимальную погрешность между требуемым и реальным коэффициентами усиления в заданных полосах частот, с использованием алгоритма замены Ремеза.
Аргументы
-
numtaps::Integer
: требуемое количество отводов в фильтре. Количество отводов — это количество членов в фильтре или порядок фильтра плюс один. -
bands_defs
: последовательность определений полос. Эта последовательность определяет полосы. Каждый элемент представляет собой пару. Первый элемент пары — это кортеж границ полосы (нижняя и верхняя). Второй элемент пары определяет требуемую характеристику и вес данной полосы. Вес не является обязательным и по умолчанию имеет значение 1,0. Требуемая характеристика и вес могут быть либо скалярами, либо функциями. Функция должна принимать вещественную частоту и возвращать вещественную требуемую характеристику или вещественный вес. Примеры:-
ФНЧ с единичными весами.
[(0, 0.475) => 1, (0.5, 1.0) => 0]
-
ФНЧ с весом 2 в полосе задержки.
[(0, 0.475) => (1, 1), (0.5, 1.0) => (0, 2)]
-
ППФ с единичными весами.
[(0, 0.375) => 0, (0.4, 0.5) => 1, (0.525, 1.0) => 0]
-
Преобразователь Гильберта.
[(0.1, 0.95) => 1]; neg=true
-
Дифференциатор.
[(0.01, 0.99) => (f -> f/2, f -> 1/f)]; neg=true
-
-
Hz::Real
: частота дискретизации в Гц. Значение по умолчанию — 1. -
neg::Bool
: обладает ли фильтр отрицательной симметрией. Значение по умолчанию — false. При значении false фильтр является четно-симметричным. При значении true фильтр является нечетно-симметричным. neg=true означает, что h[n]=-h[end+1-n]; neg=false означает, что h[n]=h[end+1-n]. -
maxiter::Integer
: (необязательно) максимальное число итераций алгоритма. Значение по умолчанию — 25. -
grid_density:Integer
: (необязательно) плотность сетки. Плотность сетки, используемая в алгоритмеremez
, имеет размер(numtaps + 1) * grid_density
. Значение по умолчанию — 16.
Возвращаемые значения
-
h::Array{Float64,1}
: одноранговый массив, содержащий коэффициенты оптимального (по минимаксному критерию) фильтра.
Примеры
Создание фильтра длиной 35 с полосой пропускания 0,15—0,4 Гц (требуемая характеристика 1) и полосами задержки 0—0,1 Гц и 0,45—0,5 Гц (требуемая характеристика 0). Обратите внимание: поведение в частотных диапазонах между этими полосами, то есть в переходных полосах, не определено.
julia> bpass = remez(35, [(0, 0.1)=>0, (0.15, 0.4)=>1, (0.45, 0.5)=>0]);
За счет изменения ширины переходных полос можно добиваться меньшей максимальной погрешности. Чем шире переходные полосы, тем ниже максимальная погрешность в заданных полосах. Вот пропускающий фильтр с той же полосой пропускания, но более широкими переходными полосами:
julia> bpass2 = remez(35, [(0, 0.08)=>0, (0.15, 0.4)=>1, (0.47, 0.5)=>0]);
Здесь мы вычисляем частотные характеристики и строим их график в дБ.
julia> using PyPlot
julia> b = DSP.Filters.PolynomialRatio(bpass, [1.0])
julia> b2 = DSP.Filters.PolynomialRatio(bpass2, [1.0])
julia> f = range(0, stop=0.5, length=1000)
julia> plot(f, 20*log10.(abs.(freqresp(b,f,1.0))))
julia> plot(f, 20*log10.(abs.(freqresp(b2,f,1.0))))
julia> grid()
Примеры из модульных тестов — стандартная (четная) симметрия
ФНЧ (фильтр низких частот) длиной 151.
julia> h = remez(151, [(0, 0.475) => 1, (0.5, 1.0) => 0]; Hz=2.0);
ФНЧ длиной 152. Нестандартный входной вес.
julia> h = remez(152, [(0, 0.475) => (1, 1), (0.5, 1.0) => (0, 2)]; Hz=2.0);
ФВЧ (фильтр высоких частот) длиной 51.
julia> h = remez(51, [(0, 0.75) => 0, (0.8, 1.0) => 1]; Hz=2.0);
ППФ (полосно-пропускающий фильтр) длиной 180.
julia> h = remez(180, [(0, 0.375) => 0, (0.4, 0.5) => 1, (0.525, 1.0) => 0]; Hz=2.0, maxiter=30);
Примеры из модульных тестов — нечетно-симметричные фильтры: Гильберта и дифференциатор.
При четной длине аппроксимация гораздо лучше, так как характеристика не ограничена нулем на частоте Найквиста. Преобразователь Гильберта длиной 20.
julia> h = remez(20, [(0.1, 0.95) => 1]; neg=true, Hz=2.0);
Преобразователь Гильберта длиной 21.
julia> h = remez(21, [(0.1, 0.95) => 1]; neg=true, Hz=2.0);
Дифференциатор длиной 200.
julia> h = remez(200, [(0.01, 0.99) => (f -> f/2, f -> 1/f)]; neg=true, Hz=2.0);
Дифференциатор длиной 201.
julia> h = remez(201, [(0.05, 0.95) => (f -> f/2, f -> 1/f)]; neg=true, Hz=2.0);
Обратный sinc-фильтр — пользовательская функция характеристики
julia> L = 64; Fs = 4800*L;
julia> passband_response_function = f -> (f==0) ? 1.0 : abs.((π*f/4800) ./ sin.(π*f/4800));
julia> h = remez(201, [( 0.0, 2880.0) => (passband_response_function, 1.0),
(10000.0, Fs/2) => (0.0, 100.0)]; Hz=Fs);
remez(numtaps::Integer, bands::Vector, desired::Vector; weight::Vector=[], Hz::Real=1.0, filter_type::RemezFilterType=filter_type_bandpass, maxiter::Integer=25, grid_density::Integer=16)
Это совместимая с SciPy версия, требующая трех аргументов (numtaps, bands, desired). Имеется также упрощенная версия API с двумя аргументами (numtaps, band_defs). Синтезируемые фильтры эквивалентны, просто входные данные задаются по-разному. Ниже приведены аргументы и примеры не для упрощенной версии API.
Аргументы
-
bands::Vector
: монотонная последовательность, содержащая границы полос в Гц. Все элементы должны быть неотрицательными и меньше половины частоты дискретизации, заданной вHz
. -
desired::Vector
: последовательность длиной в два раза меньше размера полос, содержащая требуемые коэффициенты усиления в каждой из указанных полос. -
weight::Vector
: (необязательно) относительный вес области каждой полосы. Длина вектораweight
должна быть вполовину меньше длиныbands
. -
filter_type::RemezFilterType
: значение по умолчанию —filter_type_bandpass
. Тип фильтра:-
filter_type_bandpass
: плоская характеристика в полосах. Тип по умолчанию. -
filter_type_differentiator
: частотно-пропорциональная характеристика в полосах. Нечетная симметрия, как и в случае сfilter_type_hilbert
, но с линейной наклонной требуемой характеристикой. -
filter_type_hilbert
: фильтр с нечетной симметрией, то есть фильтр с линейной фазовой характеристикой 3-го рода (для четного порядка) или 4-го рода (для нечетного порядка).
-
Примеры
Сравните примеры с упрощенным API и API SciPy. В каждом из следующих блоков сначала синтезируется фильтр с помощью упрощенного (рекомендуемого) API, а затем тот же фильтр синтезируется с помощью API, совместимого с SciPy.
julia> bpass = remez(35, [(0, 0.1)=>0, (0.15, 0.4)=>1, (0.45, 0.5)=>0]);
julia> bpass = remez(35, [0, 0.1, 0.15, 0.4, 0.45, 0.5], [0, 1, 0]);
julia> bpass2 = remez(35, [(0, 0.08)=>0, (0.15, 0.4)=>1, (0.47, 0.5)=>0]);
julia> bpass2 = remez(35, [0, 0.08, 0.15, 0.4, 0.47, 0.5], [0, 1, 0]);
julia> h = remez(20, [(0.1, 0.95) => 1]; neg=true, Hz=2.0);
julia> h = remez(20, [0.1, 0.95], [1]; filter_type=filter_type_hilbert, Hz=2.0);
julia> h = remez(200, [(0.01, 0.99) => (f -> f/2, f -> 1/f)]; neg=true, Hz=2.0);
julia> h = remez(200, [0.01, 0.99], [1]; filter_type=filter_type_differentiator, Hz=2.0);
Характеристика фильтра
#
DSP.Filters.freqresp
— Function
H, w = freqresp(filter)
Частотная характеристика H
фильтра filter
при (нормализованных) частотах w
в радианах на отсчет для цифрового фильтра или радианах в секунду для аналогового фильтра, выбранная в качестве разумного значения по умолчанию.
freqresp(filter::FilterCoefficients{:z}, w)
Частотная характеристика цифрового фильтра filter
при нормализованной частоте или частотах w
в радианах на отсчет.
freqresp(filter::FilterCoefficients{:s}, w)
Частотная характеристика аналогового фильтра filter
на частоте или частотах w
в радианах в секунду.
#
DSP.Filters.phaseresp
— Function
phi, w = phaseresp(filter)
Фазовая характеристика phi
фильтра filter
при (нормализованных) частотах w
в радианах на отсчет для цифрового фильтра или радианах в секунду для аналогового фильтра, выбранная в качестве разумного значения по умолчанию.
phaseresp(filter, w)
Фазовая характеристика фильтра filter
при (нормализованной) частоте или частотах w
в радианах на отсчет для цифрового фильтра или радианах в секунду для аналогового фильтра.
#
DSP.Filters.grpdelay
— Function
tau, w = grpdelay(filter)
Групповое запаздывание tau
фильтра filter
при (нормализованных) частотах w
в радианах на отсчет для цифрового фильтра или радианах в секунду для аналогового фильтра, выбранное в качестве разумного значения по умолчанию.
grpdelay(fliter, w)
Групповое запаздывание цифрового фильтра 'filter' при нормализованной частоте или частотах 'w' в радианах на отсчет.
#
DSP.Filters.impresp
— Function
impresp(filter, n=100)
Импульсная характеристика цифрового фильтра filter
с n
точками.
#
DSP.Filters.stepresp
— Function
stepresp(filter, n=100)
Переходная характеристика цифрового фильтра filter
с n
точками.
Прочее
#
DSP.Filters.coefb
— Function
coefb(f)
Коэффициенты числителя объекта PolynomialRatio, начиная с наибольшей степени, то есть b
, переданное в filt()
.
#
DSP.Filters.coefa
— Function
coefa(f)
Коэффициенты знаменателя объекта PolynomialRatio, начиная с наибольшей степени, то есть a
, переданное в filt()
.
Примеры
Создание эллиптического низкочастотного фильтра 4-го порядка с нормализованной частотой среза 0,2 дБ, неравномерностью характеристики в полосе пропускания 0,5 дБ и затуханием в полосе задержки 30 дБ и извлечение коэффициентов числителя и знаменателя передаточной функции:
responsetype = Lowpass(0.2)
designmethod = Elliptic(4, 0.5, 30)
tf = convert(PolynomialRatio, digitalfilter(responsetype, designmethod))
numerator_coefs = coefb(tf)
denominator_coefs = coefa(tf)
Фильтрация данных в x
, дискретизированных с частотой 1000 Гц, с помощью пропускающего фильтра Баттерворта 4-го порядка в диапазоне от 10 до 40 Гц:
responsetype = Bandpass(10, 40; fs=1000)
designmethod = Butterworth(4)
filt(digitalfilter(responsetype, designmethod), x)
Фильтрация данных в x
, дискретизированных с частотой 50 Гц, с помощью 64-отводного низкочастотного КИХ-фильтра с окном Хеннинга на частоте 5 Гц:
responsetype = Lowpass(5; fs=50)
designmethod = FIRWindow(hanning(64))
filt(digitalfilter(responsetype, designmethod), x)
Оценка порядка низкочастотного эллиптического фильтра с нормализованной частотой среза полосы пропускания 0,2 дБ, частотой среза полосы задержки 0,4 дБ, неравномерностью характеристики в полосе пропускания 3 дБ и затуханием в полосе задержки 40 дБ:
(N, ωn) = ellipord(0.2, 0.4, 3, 40)