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

Применение фильтров библиотеки DSP

Открыть пример в Engee

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

Здесь мы рассмотрим основные фильтры:

  1. фильтры с конечной импульсной характеристикой,
  2. низкочастотные фильтры,
  3. высокочастотные фильтры.

Более детально со всеми принципами построения фильтров можно ознакомиться в этом примере: https://engee.com/community/ru/catalogs/projects/sravnitelnyi-analiz-filtrov

Начнём с подключения данной библиотеки. Важно отметить, что библиотека DSP различает коэффициенты фильтра и фильтры с сохранением состояния.

In [ ]: 
Pkg.add("DSP")
using DSP
using FFTW

   Resolving package versions...
  No Changes to `~/.project/Project.toml`
  No Changes to `~/.project/Manifest.toml`

Объявим входной сигнал, который впоследствии будем применять для тестирования наших фильтров.

А также объявим вспомогательную функцию.

In [ ]: 
Fs = 1000  # Частота дискретизации
t = 0:1/Fs:1-1/Fs  # Временной интервал
x = sin.(2π .* 0.1 .* t) .+ randn(length(t))  # Сигнал с шумом (0.1 Гц компонент и шум)
plot(x)
Out[0]:
In [ ]: 
# Расчёт спектра сигнала
function compute_spectrum(signal, fs)
    n = length(signal)
    spectrum = abs.(fft(signal)) / n
    freqs = (0:n-1) .* (fs / n)
    spectrum[1:Int(n/2)], freqs[1:Int(n/2)]  # Вернуть половину спектра (для удобства)
end
Out[0]:
compute_spectrum (generic function with 2 methods)
In [ ]: 
spectrum_x, freqs_x = compute_spectrum(x, Fs)
plot(freqs_x, spectrum_x, xlabel="Frequency (Hz)", ylabel="Amplitude", legend=:topright)
Out[0]:

Далее перейдём к примерам реализации различных фильтров.

FIR-фильтры (Finite Impulse Response)

FIR-фильтры имеют конечный отклик на импульс, что делает их стабильно решаемыми и легко проектируемыми. Они не имеют внутренних обратных связей.

Тип фильтра, применяемый нами, – это низкочастотный фильтр (Lowpass) с частотой среза 5 Гц. Применён к сигналу с помощью FIR-фильтра, спроектированного методом оконного преобразования (Hanning window) длиной 64.

In [ ]: 
responsetype = Lowpass(5; fs=Fs)
designmethod = FIRWindow(hanning(64))
y = filt(digitalfilter(responsetype, designmethod), x)

plot(x, label="Original Signal", title="Signal")
plot!(y, label="Filtered Signal", xlabel="Sample", ylabel="Amplitude")
Out[0]:
In [ ]: 
spectrum_y, freqs_y = compute_spectrum(y, Fs)

plot(freqs_x, spectrum_x, label="Original Signal Spectrum", xlabel="Frequency (Hz)", ylabel="Amplitude", legend=:topright)
plot!(freqs_y, spectrum_y, label="Filtered Signal Spectrum", title="Signal Spectrum")
Out[0]:

После фильтрации сигнал стал сглаженным, подавляя шумовые компоненты. Форма сигнала теперь лучше отражает низкочастотную составляющую.

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

Частоты выше 32 Гц полностью подавлены. Основные амплитуды спектра сосредоточены на частотах 0 Гц (DC-компонента) и 3 Гц.

Это говорит о том, что фильтр эффективно подавил шум выше частоты среза.

Низкочастотные фильтры (Low-pass filters)

Эти фильтры пропускают сигналы с частотами ниже заданного порога и ослабляют сигналы с более высокими частотами.

Тип фильтра, применяемый нами, – это низкочастотный фильтр (Lowpass) с частотой среза 0.2 Гц. Спроектирован методом эллиптического фильтра 4-го порядка с допустимыми отклонениями 0.5 дБ в полосе пропускания и 30 дБ в полосе заграждения. Реализован в виде передаточной функции с использованием полиномиального представления (PolynomialRatio).

In [ ]: 
responsetype = Lowpass(0.2)
designmethod = Elliptic(4, 0.5, 30)
tf = convert(PolynomialRatio, digitalfilter(responsetype, designmethod))
numerator_coefs = coefb(tf)
denominator_coefs = coefa(tf)

y = filt(tf, x)

plot(x, label="Original Signal", title="Signal")
plot!(y, label="Filtered Signal", xlabel="Sample", ylabel="Amplitude")
Out[0]:
In [ ]: 
spectrum_y, freqs_y = compute_spectrum(y, Fs)

plot(freqs_x, spectrum_x, label="Original Signal Spectrum", xlabel="Frequency (Hz)", ylabel="Amplitude", legend=:topright)
plot!(freqs_y, spectrum_y, label="Filtered Signal Spectrum", title="Signal Spectrum")
Out[0]:

Фильтр справился с задачей хуже, чем предыдущий, так как видны остаточные шумы в основном сигнале. Сглаживание есть, но его качество неудовлетворительно, поскольку помехи продолжают оказывать значительное влияние на итоговый результат.

Спектр обработанного сигнала показывает наличие пиков на частотах вплоть до 130 Гц, что свидетельствует о недостаточной фильтрации высокочастотных составляющих: фильтр неэффективно подавляет частоты выше своей частоты среза.

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

Частота среза (0.2 Гц) и выбранные параметры фильтра, возможно, не оптимальны для этого сигнала. Кроме того, порядок фильтра (4-й) может быть недостаточен для достижения требуемого подавления помех.

Высокочастотные фильтры (High-pass filters)

Эти фильтры пропускают сигналы с частотами выше заданного порога и ослабляют сигналы с более низкими частотами.

Тип фильтра, который мы здесь применили, – это высокочастотный фильтр (High-pass) с частотой среза 40 Гц, реализованный с помощью фильтра Баттерворта 4-го порядка.

In [ ]: 
responsetype = Highpass(40; fs=Fs)  # Высокочастотный фильтр, отсечка на 40 Гц, fs = 1000 Гц
designmethod = Butterworth(4)         # Баттерворт, порядок 4
highpass_filter = digitalfilter(responsetype, designmethod)

y = filt(highpass_filter, x)  # Применяем фильтр к сигналу x

plot(x, label="Original Signal", title="Signal")
plot!(y, label="Filtered Signal", xlabel="Sample", ylabel="Amplitude")
Out[0]:
In [ ]: 
spectrum_y, freqs_y = compute_spectrum(y, Fs)

plot(freqs_x, spectrum_x, label="Original Signal Spectrum", xlabel="Frequency (Hz)", ylabel="Amplitude", legend=:topright)
plot!(freqs_y, spectrum_y, label="Filtered Signal Spectrum", title="Signal Spectrum")
Out[0]:

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

Спектр показывает, что все частоты до 40 Гц (включая значимые низкочастотные компоненты до 24 Гц) полностью подавлены. Это привело к утрате основной информации о сигнале, оставив только высокочастотный шум.

Высокочастотный фильтр с отсечкой на 40 Гц совершенно не подходит для анализа данного сигнала, поскольку значимая частота сигнала (0.1 Гц) лежит в диапазоне низких частот, которые фильтр удаляет.

Полосовые фильтры (Band-pass filters)

Эти фильтры пропускают сигналы в пределах заданной полосы частот, ослабляя сигналы за пределами этой полосы. Тип фильтра, который мы здесь применяем, – это полосовой фильтр с полосой пропускания от 10 до 40 Гц, реализованный с использованием фильтра Баттерворта 4-го порядка.

Частота дискретизации: 1000 Гц.

In [ ]: 
responsetype = Bandpass(10, 40; fs=Fs)
designmethod = Butterworth(4)
y = filt(digitalfilter(responsetype, designmethod), x)

plot(x, label="Original Signal", title="Signal")
plot!(y, label="Filtered Signal", xlabel="Sample", ylabel="Amplitude")
Out[0]:
In [ ]: 
spectrum_y, freqs_y = compute_spectrum(y, Fs)

plot(freqs_x, spectrum_x, label="Original Signal Spectrum", xlabel="Frequency (Hz)", ylabel="Amplitude", legend=:topright)
plot!(freqs_y, spectrum_y, label="Filtered Signal Spectrum", title="Signal Spectrum")
Out[0]:

Фильтр сгладил шумы, что заметно при сравнении обработанного сигнала с исходным зашумлённым. Однако полезная информация о низкочастотных колебаниях сигнала была утрачена. Остался лишь «фрагмент» высокочастотной части сигнала, который не несёт значимой информации.

Спектр обработанного сигнала показывает, что фильтр эффективно подавил частоты ниже 10 Гц и выше 70 Гц, но сохранил диапазон от 10 до 40 Гц, который в данном случае не содержит полезной информации. Это привело к тому, что сигнал оказался представлен шумовыми компонентами в этой полосе, потеряв изначальный низкочастотный характер.

Полосозаграждающие фильтры (Band-stop filters)

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

В данном случае мы реализуем полосо-заградительный фильтр, подавляющий частоты в диапазоне 10–40 Гц и реализованный с использованием фильтра Баттерворта 4-го порядка.

In [ ]: 
responsetype = Bandstop(10, 40; fs=Fs)  # Полосо-заградительный фильтр, подавляет частоты от 10 до 40 Гц, fs = 1000 Гц
designmethod = Butterworth(4)            # Баттерворт, порядок 4

y = filt(digitalfilter(responsetype, designmethod), x)  # Применяем фильтр к сигналу x (определите x заранее)


plot(x, label="Original Signal", title="Signal")
plot!(y, label="Filtered Signal", xlabel="Sample", ylabel="Amplitude")
Out[0]:
In [ ]: 
spectrum_y, freqs_y = compute_spectrum(y, Fs)

plot(freqs_x, spectrum_x, label="Original Signal Spectrum", xlabel="Frequency (Hz)", ylabel="Amplitude", legend=:topright)
plot!(freqs_y, spectrum_y, label="Filtered Signal Spectrum", title="Signal Spectrum")
Out[0]:

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

Спектр показывает подавление частот в диапазоне от 12 до 32 Гц. Остальные частотные компоненты, включая шум, остались практически неизменными.

Это объясняет сохранение зашумленного характера сигнала.

Вывод

Мы рассмотрели различные фильтры для решения данной задачи и выяснили, что FIR-фильтр с Hanning окном оказался лучшим выбором. Его способность плавно и эффективно удалять шум выше заданной частоты среза позволила выделить полезную низкочастотную информацию (0.1 Гц), что было основной задачей.

Обобщая проведённые выще сравнения, можно сделать следующий вывод: для подобных задач стоит использовать низкочастотные фильтры с частотой среза, оптимально соответствующей диапазону полезной информации.