Разработка цифровых фильтров на Python и C ++

Разработка цифровых фильтров на Python и C ++

В следующей статье я продемонстрирую общий подход к созданию цифровых фильтров. Цель фильтра - удалить из сигнала определенные частоты (шум). Мы увидим, как разработать проходной фильтр (фильтр нижних частот) и более продвинутый, узкополосный фильтр. Проектирование будет выполняться на Python с использованием в основном мощного пакета управления (пакет Python, который в достаточной степени реализует стандартные функции для анализа и проектирования систем управления - аналогично Matlab). Развертывание программного обеспечения и окончательная проверка конструкции фильтров будет производиться на языке C ++. Для построения графиков (на C ++) я буду использовать также обсуждавшуюся в моих предыдущих публикациях библиотеку matlablib для C ++. Эта библиотека будет использоваться только для построения графиков.

Заголовочный файл (для библиотеки печати) должен находиться в той же папке, что и ваш cpp.
Каждая программа, которую вы компилируете следующим образом,

//compile
g++ my_prog.cpp -o my_prog -I/usr/include/python3.8 -lpython3.8// //run
./my_prog//folder tree
├── my_prog
├── my_prog.cpp
├── matplotlibcpp.h

Введение

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

Для простоты мы можем предположить, что мы захватываем аудиосигнал (обычно это может быть сигнал от процесса, но здесь мы рассматриваем аудиосигнал). Сигнал очень зашумленный, поэтому его содержание очень трудно понять. Используя описанную в этой статье методологию проектирования фильтров и общие принципы обработки сигналов, мы постараемся удалить шум и сохранить сигнал контента без изменений.
Обратите внимание, что в этой статье я имею в виду значение Гц (например, 50 Гц), но все вычисление было выполнено на основе значений, умноженных на 2 * пи.

Ниже я изобразил «прямоугольный» сигнал, который состоит (мы обсудим позже) только из синусоидальных волн. Я также отобразил БПФ такого сигнала.

Исходный код для всех симуляций (Python и C ++) вы найдете здесь.

// signal definitions
sig1 = sin(t) + sin(3*t)/3 + sin(5*t)/5
sig2 = sin(t) + sin(3*t)/3 + sin(5*t)/5 + sin(7*t)/7 + sin(9*t)/9 + sin(11*t)/11 + sin(13*t)/13

FT для sig2 можно изобразить следующим образом:

Самая важная концепция в обработке сигналов - это Преобразование Фурье (FT), которое преобразует сигнал из временной области в частотную область (спектр). Используя FT, мы можем снова применить преобразование и восстановить сигнал в исходной области (времени).

Каждый сигнал во временной области мы обычно можем представить, используя два разных подхода.
Первый учитывает дискретные значения (амплитуды) сигнала в определенных точках. Далее эти дискретные значения (полученные / измеренные в определенный момент времени) суммируются. Конечное значение - это сигнал, который мы рассматриваем (пожалуйста, обратите внимание на изображение ниже, которое изображает принцип этого подхода).
Второй подход значительно расходится, поскольку вместо суммирования дискретных значений, взятых за определенное время, метод суммирует функцию это покрытие все время. Эти функции можно рассматривать как волны и математически представить как уравнение Эйлера, где омега выражает определенную частоту волны.

Как вы можете видеть на изображении выше и формуле FT, каждый сигнал может быть представлен суммированием волн, рассчитанных для разных частот.
Количество синусоидальных переплетений (рассчитанных для определенной частоты), которые суммируются, напрямую влияет на структуру сигнала (форма).
Вы можете себе представить, увеличение количества волн (пример выше) сделает окончательный сигнал более похожим на прямоугольный.
Обычно (это основной принцип теории преобразования Фурье) каждый сигнал может быть представлен суммой синусоидальных волн. Вычисляется (добавляется) больше волн, может быть достигнуто лучшее математическое представление («основного» / обработанного) сигнала.

Методология обработки сигналов для представления сигнала во временной области в виде композиции синусоидальных волн довольно сложна. Обычно неудобно визуализировать и обрабатывать.
Более подходящий подход - представить сигнал в частотной области, что может быть достигнуто с помощью преобразования Фурье. FT (как я уже упоминал, работает в обоих направлениях) может просто извлекать доминирующие частоты (основного / исходящего) сигнала. Мы можем сказать, что FT будет усиливать доминирующие частоты основного сигнала, поэтому будет проще использовать постобработку, чтобы (например) удалить определенную частоту, микшировать или усилить.
Применение FT для сигнала (в во временной области) мы просто переносим (отображаем) сигнал) из временной области в частотную, где основной сигнал представлен амплитудами композиционных волн основного сигнала (для определенной частоты).

ДИЗАЙН ЦИФРОВОГО ФИЛЬТРА

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

Обычно компьютерная система или DSP (цифровой сигнальный процессор) обрабатывают значения цифрового сигнала, взятые из аналогового сигнала в процессе выборки. Все операции, выполняемые в системе DSP, являются дискретными (это означает, что они происходят в определенное время - амплитуда аналогового сигнала снимается / измеряется с интервалами времени T).
См. Изображение ниже, которое представляет аналоговый и дискретный сигнал.

Теоретически все системы с дискретным временем могут быть описаны разностным уравнением N-го порядка, где x (n) - вход системы, а y (n) - выход. m и параметры являются фиктивными значениями и используются только для суммирования.

Уравнение имеет эквивалентный вид и может быть переписано как

Обратите внимание, что значения y (n-m) и x (n-p) являются отложенной функцией выхода y (n) и входа x (n). См. Пример ниже, который представляет дискретный сигнал задержки.

Цифровой фильтр описывается той же передаточной функцией, что и дискретная система, представленная ранее в этой статье. Filer, как мы можем себе представить, позволяет передавать определенную частоту (синусоидальную составляющую входного «основного» сигнала на определенной частоте) без изменений. Другие компоненты (частоты), на которые рассчитан фильтр, будут ослаблены.

Существует два типа цифровых фильтров: Фильтр с конечной импульсной характеристикой (FIR), где коэффициенты равны 0 (этот тип фильтра также называется авторегрессивным - AR ) и Фильтр бесконечной импульсной характеристики (IIR), где коэффициенты b передаточной функции равны 0 (также называемое скользящей средней MA). .
Как и следовало ожидать, мы можем объединить эти два типа фильтров, доработав ARMA фильтр.

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

Здесь стоит упомянуть, что в дискретной области каждый сигнал преобразуется (время в частотную область) с помощью дискретного преобразования Фурье.
DFT переводит сигнал во временной области в частотную, и этот процесс можно описать следующим образом

Мы будем использовать эту формулу и вычислить ДПФ для конкретных сигналов на C ++.

В следующем разделе будут показаны общие шаги, которые я применил для создания фильтра (низкий проход и выемка). Как я уже упоминал, для проектирования и проверки я использовал Python (и тоже рекомендую). Во-вторых, я применил C ++, где развернул дизайн фильтра (коэффициенты).

Действия, которые необходимо выполнить, можно резюмировать следующим образом:

1. Определите входной (основной) сигнал и график.
2. Определите частоту, которую вы хотите удалить (полоса пропускания или режекторный фильтр).
3. Вычислить ДПФ для входного сигнала.
4. Вычислить передаточную функцию фильтра.
5. Вычислить Body plot для передаточной функции фильтра и проверить характеристики (пропускание и полоса затухания).
6. Вычислите дискретную передаточную функцию для фильтра, используя метод Тастина.
7. Вычислить коэффициенты фильтра.
8. Примените коэффициенты фильтра к общей передаточной функции ARMA (проектный фильтр).
9. Вычислите ДПФ для выходного сигнала (проверьте фильтр).
10. Шаги 1, 3, 8, 9 и 10 реализованы на C ++.

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

В следующем разделе будет описано в нескольких шагах, как разработать фильтр нижних частот. Как я уже упоминал ранее, проектирование выполняется на Python с использованием пакета Control (аналогично Matlab) с поддержкой Scipy.

ЦИФРОВОЙ ФИЛЬТР НИЗКОГО ПРОХОДА

Шаг 1.

В следующих примерах мы собираемся удалить частоту сигнала 50 Гц из входного сигнала. Удаление этой конкретной частоты является обычным явлением, поскольку есть частота источника питания и довольно часто присутствует, если мы проводим измерения.

Входной сигнал состоит из двух синусоидальных волн (5 Гц и 50 Гц), которые показаны на рисунке ниже.

//signal definition
A1 = 1
A2 = 0.2
f1 = 5
f2 = 50
y = A1*sin(2*pi*f1*t) + A2*sin(2*pi*f2*t)

Здесь мы используем БПФ для построения сигнала в частотной области. Как видите, есть спектр мощности сигнала с двумя доминирующими частотами (DFT мы используем в C ++).

Шаг 2.

Здесь мы используем фильтр нижних частот. Фильтр удалит все частоты выше частоты среза, которая в нашем случае составляет 5 Гц. Помните, наша цель - уменьшить частоту 50 Гц. Обратите внимание, что природа не «обеспечила» идеальный «дизайн», поэтому характеристика пропускания / затухания фильтра не будет идеальной. В зависимости от конструкции и физического материала фильтр состоит из наблюдаемой «переходной» полосы фильтра, где для определенных частот сигнал не будет полностью ослаблен. Пожалуйста, рассмотрите график Body и узнайте, как будет работать фильтр (Шаг 5).

Шаг 3.

Дискретное преобразование Фурье будет применяться при реализации проекта на C ++. В наших целях ДПФ или в целом (преобразование Фурье) используется для визуализации. Основное различие между FT и DTF связано с областью применения этих преобразований.
Помимо преобразования сигнала карты во временной области в частотную область, DFT работает с дискретными сигналами и отображает входной (дискретный) сигнал на определенные дискретные частоты. .

Реализацию DFT (на C ++) вы найдете в файле cpp (см. Также в и этой статьи. Формула, используемая для вычисления DFT (амплитуды составляющих сигнала для определенной частоты), может быть записана следующим образом:

Вычисляя ДПФ (в C ++) для определенного входного сигнала, вы ожидаете получить следующие результаты:

Шаг 4.

Фильтр нижних частот описывается следующей простой передаточной функцией

Частота среза в нашем случае 5 Гц. Для этого фильтра и частоты среза мы можем вычислить с помощью пакета Control (метод Pole-Zero Map - который используется для анализа стабильности системы = ›система может быть ассоциирована как стабильная, если все полюса расположены слева от реальной ось. Стабильная система - это система, которая возвращается в нормальное состояние после приложения силы или других возмущений). Фильтр не имеет нулей, но имеет один полюс, равный… .. (

Шаг 5.

Теперь для передаточной функции фильтра мы вычисляем Body plot, см. Результаты ниже,

Шаг 6.

Выполнение следующего шага дает нам обратную связь о коэффициенте цифрового фильтра. Чтобы преобразовать нашу непрерывную передаточную функцию в дискретную (в связи с тем, что мы разрабатываем цифровой фильтр), мы используем Билинейное преобразование, также называемое методом Тастина.

Преобразование простое. Заменим s следующей функцией:

Шаг 7.

Вычисление коэффициента дискретной передаточной функции фильтра можно производить вручную, однако мы будем использовать Python. Применяя метод to_discrete (), Python возвращает значения коэффициентов.

Шаг 8.

На следующем шаге мы вычислим функцию окончательного фильтра, которая будет проверена на Python, а затем на C ++.

Вычисленные на предыдущем шаге коэффициенты позволяют нам написать окончательную передаточную функцию фильтра, которая соответствует нашим требованиям (ослабить все частоты выше 5 Гц). (подробности см. в Jupyter Notebook)

Шаг 9

Запуск нашего моделирования, при котором входной сигнал (5 Гц и 50 Гц) проходит через разработанный фильтр. Результаты ожидаются. Синусоидальная волна 50 Гц была удалена (ослаблена). На следующих рисунках показан выходной сигнал фильтра в частотной области (FT в Python, DFT в C ++) и сигнал (после фильтра) во временной области.

In C++.

ДИЗАЙН ФИЛЬТРА С ЗАЖИМОМ

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

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

В предыдущем разделе мы обсуждали фильтр нижних частот. Эта концепция будет использоваться повторно, однако с важными изменениями. Вы могли видеть, что фильтр нижних частот удаляет волну 50 Гц. Но как быть в ситуации, когда мы бы хотели сохранить частоту 80 Гц !?
Применение для этой задачи стандартного фильтра нижних частот уберет и эту частоту. Нужно ли применять альтернативное решение? да.
Режекторный фильтр удалит (в данном конкретном случае) только одну определенную частоту - 50 Гц.

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

График тела низких частот.

С другой стороны, как я уже упоминал выше, мы ищем «узкополосный» фильтр, который ослабляет только определенную частоту. Для этого нам нужно взглянуть на фильтр высоких частот (например, 20 Гц) - это фактическое значение 20рад / с, как я упоминал ранее.

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

Для нашего конкретного примера график Body для этого типа фильтра может быть построен и изображен:

Сюжет High Band Body.

Теперь мы объединим передаточную функцию и график Body plot.

График тела для полосового фильтра.

Как видите, существует своего рода узкополосный фильтр, но у нас нет надлежащего контроля над глубиной фильтра, поэтому мы должны немного переосмыслить наше решение.

Другой и достаточный подход - снова взять фильтр нижних частот, но на этот раз второго порядка (двухполюсный) и понизить демпфирование. Посмотрите, что происходит с сюжетом Body.

Как вы видете. Уменьшая демпфирование почти до нуля, мы получаем почти бесконечное усиление для частоты среза (Python, к сожалению, здесь не вычисляет должным образом).

Теперь применим трюк и перевернем передаточную функцию.

График Body для этой новой (перевернутой передаточной функции) выглядит следующим образом:

Это почти идеальная передаточная функция, которую мы ищем, однако у нас есть еще две проблемы, с которыми нам нужно справиться:
1. высокочастотные сигналы, которые усиливаются на 40 дБ / декаду, поскольку у нас сейчас только два нуля (без полюсов) и,
2. передаточная функция не может быть физически реализована, так как числитель передаточной функции фильтра больше знаменателя

Мы решим эти задачи, добавив два полюса. Каждый полюс будет снижать амплитуду высоких частот на 20 дБ / декаду и выравнивать нашу частотную характеристику по сравнению с собственной частотой.

КАК ВЫБРАТЬ ЭТИ ПОЛЮСА?

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

pole_above = (a * wn)/(s+a* wn)
pole_lower = (a * wn)/(s+(wn /a))

Применяя их к полюсам для нашей перевернутой передаточной функции, мы получим:

Мы можем проверить наш дизайн фильтра, снова построив график Body.

График тела для окончательного режекторного фильтра,

Теперь фильтр намного глубже - для нашей частоты сигнал будет ослаблен примерно на 28 дБ.

Мы только что закончили проектирование нашего фильтра Notch. Теперь мы должны применить те же шаги проектирования, которые я указал и объяснил выше, пока мы разрабатывали фильтр нижних частот. Принцип тот же. Подробности вы найдете в блокноте jupyter.

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

Для первого сигнала, как и ожидалось, 50 Гц были удалены. Однако наиболее интересным является второй сигнал, который состоит из трех частот: 20, 50 и 80 Гц. Notch-фильтр в этом случае работает отлично. Снимает почти грамотно 50 Гц, не влияя на боковые частоты 20 и 80 Гц.

Пожалуйста, обратите внимание на результаты ниже (на Python и C ++).

Рассматриваемый сигнал после Notch-фильтра (во временной области).

Спектральный ДПФ относительного сигнала (50 Гц) удален.

Спасибо за чтение,

Если эта статья была вам интересна, я действительно рекомендую вам проверить эти два впечатляющих репозитория:

// https://github.com/capitanov/dsp-theory
// https://github.com/AllenDowney/ThinkDSP