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.Butterworth — Function

Butterworth(n)

-полюсный фильтр Баттерворта.

# 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 (по умолчанию), синтезируемый КИХ-фильтр масштабируется так, чтобы было верно следующее.

Для фильтров Lowpass и Bandstop частотная характеристика равна единице при 0 (прямой ток).
Для фильтров Highpass частотная характеристика равна единице при частоте Найквиста.
Для фильтров Bandpass частотная характеристика равна единице в центре полосы пропускания.

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.iirnotch — Function

iirnotch(Wn, bandwidth[; fs])

Цифровой запирающий БИХ-фильтр второго порядка ^[4] на частоте Wn с шириной полосы bandwidth. Если аргумент fs не указан, Wn интерпретируется как нормализованная частота в полуциклах на отсчет.

Характеристика фильтра

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

1. О. Херрманн (Herrmann, O.), Лоренс Р. Рабинер (Lawrence R. Rabiner) и Д. С. К. Чан (D. S. K. Chan), Practical design rules for optimum finite impulse response lowpass digital filters. Bell System Technical Journal 52.6 (1973): 769—799.

2. Дж. Х. Макклеллан (J. H. McClellan) и Т. В. Паркс (T. W. Parks), A unified approach to the design of optimum FIR linear phase digital filters, IEEE Trans. Circuit Theory, том CT-20, стр. 697—701, 1973 г.

3. Дж. Х. Макклеллан (J. H. McClellan), Т. В. Паркс (T. W. Parks) и Л. Р. Рабинер (L. R. Rabiner), A Computer Program for Designing Optimum FIR Linear Phase Digital Filters, IEEE Trans. Audio Electroacoust., том AU-21, стр. 506—525, 1973 г.

4. С. Дж. Орфанидис (Orfanidis, S. J.), 1996 г., Introduction to signal processing, Englewood Cliffs, N.J: Prentice Hall, стр. 370.