Синтез Цифрового БИХ Фильтра Низких Частот / forpes.ru

Главная
Синтез Цифрового БИХ Фильтра Низких Частот

Синтез Цифрового БИХ Фильтра Низких Частот +11

14.10.2024 13:50

aabzel 15 1600 Источник

Пролог

Недавно мне потребовалось синтезировать быстрый цифровой фильтр нижних частот. Причем этот фильтр должен работать в реальном времени на микроконтроллере. Тут я понял, что надо вспоминать с какой стороны следует подходить к цифровым IIR фильтрам.

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

Обратная задача

Дан цифровой БИХ фильтр первого порядка с известными коэффициентами a_0 и b_1.

Также известна частота дискретизации F_s. Требуется построить АЧХ, увидеть и оценить частоту среза.

Прямая задача

Дана частота дискретизации, дана частота среза F_сut. Дана топология фильтра первого порядка. Найти коэффициенты IIR фильтра a0 и b1.

IIR фильтры - это тот случай, когда обратная задача много проще прямой задачи.

Терминология

Частота среза — это такая частота, на которой ослабление фильтра равно -3 дБ в логарифмическом масштабе (в линейном это 0,707).

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

На схеме фильтра Z^(-1) это звено задержки семпла. Физически это регистр в RAM памяти, разрядностью равный разрядности семпла. Одна ячейка памяти. Каждый раз, когда число проходит через звено задержки, оно тратит на это T_a секунд. T_a - это время между семплами, величина обратная частоте дискретизации. Математически звено выглядит вот так.

$z^{-1}=e^{-j2\pi\frac{f}{f_a}} \qquad (1)$

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

$y_n=a_0x_n+b_1y_{n-1} \qquad (2)$

Это так называемое разностное уравнение. У звена запаздывания есть свойство (3)

$y_{n}z^{-1}=y_{n-1} \qquad (3)$

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

$\frac{y_n}{x_n}= \frac{a_0}{1-b_1z^{-1}} \qquad (4)$

Для перехода в частотную область в (4) подставляем (1). Получается

$A(j \omega) = \frac{a_0}{1-b_1e^{-j\omega T_a}} \qquad (5)$

Меня же интересует амплитудно-частотная характеристика в диапазоне от нуля до половины частоты дискретизации f_a (Далее просто зеркальное отражение будет).

$A(f) =\left| \frac{a_0}{1-b_1e^{-j 2 \pi f/f_a}}\right| \qquad (6)$

Обратите внимание на знак минус в знаменателе. При расчете более высоких порядков в знаменателе все слагаемые будут со знаком минус.

Эти вычисления также справедливы и для FIR фильтров, так как КИХ фильтры можно считать подмножеством БИХ фильтров, если обнулить коэффициенты в обратной связи.

Вычисления

Далее я буду численно вычислять по формуле (6) используя Cи в качестве калькулятора.

Си-код вычислений АЧХ

static double complex filter_delay_link(uint32_t n, double F_hz, 
                                        double F_sam, bool is_norm_freq){
    double complex exp = 0.0;
    double omega = 0,0;
    double T_a = 0,0;
    if (is_norm_freq) {
        T_a = 1.0;
        omega = 2.0*M_PI*F_hz;
    } else {
         T_a = 1.0/F_sam;
         omega = 2.0*M_PI*F_hz;
    }
    exp = cos(n*omega*T_a) - I *sin(n*omega*T_a);
    return exp;
}


static double complex  iir_calc_feed_forward_ll(IirHandle_t* Node,
                                                double f_hz, 
                                                bool is_norm_freq){
    double complex numerator = 0.0;
    uint32_t n = 0;
    for(n=0; n<Node->size; n++) {
         numerator += Node->b[n]*filter_delay_link(n,
                                                   f_hz,
                                                   Node->sample_rate_hz,
                                                   is_norm_freq);
    }
    return numerator;
}

static double complex iir_calc_feed_back_ll(IirHandle_t* Node,
                                            double f_hz, 
                                            bool is_norm_freq) {
    double complex denominator = 1.0;
    uint32_t n = 0;
    for(n=1; n<Node->size; n++) {
        denominator += Node->a[n]*filter_delay_link(n,
                                                    f_hz,
                                                    Node->sample_rate_hz,
                                                    is_norm_freq);
    }
    return denominator;
}

bool iir_calc_frequency_response(uint8_t num){
    bool res = false;
    char* file_name = "IIRFrequencyResponse.csv";
    file_pc_delete(file_name);
    IirHandle_t* Node = IirGetNode(num);
    if (Node) {
        LOG_INFO(IIR,"%s", IirNodeToStr(Node));
        double f_hz = 0;
        for(f_hz=0; f_hz < (Node->sample_rate_hz/2.0);) {
            double f_real_hz =(double) f_hz;
            double complex numerator   = iir_calc_feed_forward_ll(Node,
                                                                  f_real_hz,
                                                                  false);
            double complex denominator = iir_calc_feed_back_ll(Node,
                                                               f_real_hz,
                                                               false);
            double complex Amplitude = numerator/denominator;

            char text[300] = "?";
            strcpy(text, "");
            snprintf(text, sizeof(text), "%s%f,", text, f_real_hz);
            snprintf(text, sizeof(text), "%s%f", text, cabs(Amplitude));
            res = file_pc_print_line(file_name, text, strlen(text));
            if(f_hz<10.0){
               f_hz +=0.1;
            }else {
               f_hz +=2.0;
            }
        }
        char command[300]="";
        snprintf(command, sizeof(command),
                 "python.exe plot_csv_file.py %s frequency Amplitude",
                 file_name);
        res = win_cmd_run( command);
    }
    return res;
}


bool iir_calc_frequency_response_norm(uint8_t num){
    bool res = false;
    char* file_name = "IIRFrequencyResponse.csv";
    file_pc_delete(file_name);
    IirHandle_t* Node = IirGetNode(num);
    if (Node) {
        LOG_INFO(IIR,"%s", IirNodeToStr(Node));
        double f_norm = 0;
        for(f_norm=0; f_norm < (0.5);) {
            double complex numerator   = iir_calc_feed_forward_ll(Node,f_norm, true);
            double complex denominator = iir_calc_feed_back_ll(Node,f_norm, true);
            double complex Amplitude = numerator/denominator;
            char text[300] = "?";
            strcpy(text, "");
            snprintf(text, sizeof(text), "%s%f,", text, f_norm);
            snprintf(text, sizeof(text), "%s%f", text, cabs(Amplitude));
            res = file_pc_print_line(file_name, text, strlen(text));
            if(f_norm<0.01){
            	f_norm +=0.00001;
            }else {
            	f_norm +=0.001;
            }
        }
        char command[300]="";
        snprintf(command, sizeof(command),
                 "python.exe plot_csv_file.py %s frequency Amplitude", file_name);
        res = win_cmd_run( command);
    }
    return res;
}

Получился вот такой график

Тут частота нормирована. То есть по X не частота, а F/F_cut. Вот мы и решили обратную задачу. По данному фильтру оценили его частоту среза.

В общем виде любой цифровой фильтр определяется двумя массивами действительных чисел A и B.

IIR Low Pass Filter, Fc:100 Hz, Fs:48 kHz
#Feed Forward
A [  
  0.00000000000815456252,
  0.00000000004892737515,
  0.00000000012231843787,
  0.00000000016309125049,
  0.00000000012231843787,
  0.00000000004892737515,
  0.00000000000815456252]

#Feed Back
B [ 1.00000000000000000000,
    -5.87230126428934080000,
    14.36924807876324900000,
    -18.75369891209695100000,
    13.76862725126773400000,
    -5.39164367598673080000,
    0.87976852284442164000]

Далее, запустив тот же алгоритм, получается вот другая AЧХ.

Есть ещё вот такой БИХ фильтр низких частот первого порядка. Он хорош тем, что для его настройки тоже существует только одна переменная k.

Для этого IIR фильтра чем выше k, тем ниже частота среза. При этом k начинается с единицы.

Решение прямой задачи.

Что касается решения прямой задачи, то тут в теории всё сложнее. Зато можно воспользоваться утилитой WinFilter. Однако надо отметить, что утилита не поможет сгенерировать IIR фильтр нижних частот с малой частотой среза, например 0,0002% от Fs. На малых частотах утилита почему-то генерирует не фильтры, а усилители. Однако если повезёт то утилита рассчитает вам нужные коэффициенты.

Коэффициенты стабильных IIR фильтров порой и вовсе табличные величины которые можно найти в учебниках и справочниках.

Достоинства цифровых IIR фильтров:

1--IIR фильтры требуют меньше вычислений в сравнении с FIR фильтрами

Недостатки цифровых IIR фильтров

--IIR фильтр порой усиливает сигнал. А это совсем не нужно. Вот пример такого фильтра

 A [] = {
        0.00000000816384086451,
        0.00000002449152259353,
        0.00000002449152259353,
        0.00000000816384086451
    };

B [] = {
        1.00000000000000000000,
        -2.99715048309490010000,
        2.99430502461701350000,
        -0.99715453863406001000
    };

И его АЧХ

Приложения цифровых фильтров нижних частот.

1--Обработка сигналов с датчиков физических величин: акселерометры, гироскопы, расстояния.

2--SDR обработка семплов с радиоприёмников. Фильтрация удвоенной частоты на выходе квадратурного смесителя.

Итоги

Удалось научиться пользоваться двумя разновидностями IIR фильтров нижних частот первого порядка.

Появилось некоторое понимание с какой стороны подходить к цифровым БИХ фильтрам.

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

Ссылки

#	Hyperlinks
1	Equation Editor
2	WinFilter
3	Свёртка
4	https://editor.codecogs.com/
5	ЦОС Окна (весовые функции).
10	Filter frequency response grapher
9	Программная реализация БИХ-фильтра в информационно-измерительном канале
6	https://www.controlpaths.com/2023/07/15/designing-iir-filters-in-python/
7	Цифровая обработка сигналов, Стивен Смит, 2018г, Москва, изд. Додэка
8	Полупроводниковая схемотехника, Том 2, 2007 г, Ульрих Титце, Кристоф Шенк

Словарь

Акроним	Расшифровка
ЦОС	Цифровая обработка сигналов
DSP	digital signal processor
AЧХ	Амплитудно-частотная характеристика
IIR	Infinite impulse response
БИХ	Фильтр с бесконечной импульсной характеристикой

Комментарии (15)

sinc_func
14.10.2024 14:23
#27414866
Фильтр всегда надо разбивать на биквадратнные секции, чтобы не иметь гемороя с уходом позиционирования полюсов и нулей. На практике желательно писать в целочисленном формате, а точнее - во fractional. Есть масса тонкостей с минимизацией шума вычислений.

Везде, где я видел описание фильтров, в - это FIR часть коэффициентов, а - это рекурсивная часть (a0 = 1.0).

Если то что обозначено как B[] заценить с точки зрения деградации вычислений то на вскидку выглядит не очень хорошо (засовываем постоянку и смотрим чувствительность от микроприращений).
1. nikolz
  14.10.2024 14:23
  #27414876
  Фильтр всегда надо разбивать на биквадратнные секции, чтобы не иметь гемороя с уходом позиционирования полюсов и нулей.
  
  если правильно понял, то:
  
  Итоги
  
  Удалось научиться пользоваться двумя разновидностями IIR фильтров нижних частот первого порядка.
  
  разбивать пока нечего.
  1. sinc_func
    14.10.2024 14:23
    #27414978
    B [] = { 1.00000000000000000000, -2.99715048309490010000, 2.99430502461701350000, -0.99715453863406001000
    
    Коэффициент А1 я всегда держу меньше 2.0.
    
    Если мы скормим единицы в рекурсивную часть
    
    -2.99715048309490010000, 2.99430502461701350000, -0.99715453863406001000
    
    то для расчета в double все конечно прокатит, но во float все выглядит чуть проблемнее
    
    А у меня везде аккумулятор в биквадрате - это int32_t.
    
    b0, b1, b2, a1 и a2 - int16_t
    
    nikolz
    14.10.2024 14:23
    #27415082
    Вы задали полосу ФНЧ 0,0002% от частоты дискретизации Fs=44000 гц. Верно?
    
    Полоса ФНЧ равна 5 Гц. Длительность переходного процесса такого фильтра составит примерно 1 секунду. Для расчета АЧХ вам надо взять длительность сигнала примерно раз в 10 больше т е примерно 500 000 отсчетов.
    
    Программа Вам написала: уменьшите частоту дискретизации или увеличьте точность расчета коэффициентов. Т е на микропроцессоре вы такой фильтр так просто не сделаете. Да и сигнал такой Вы из-за шумов не увидите.
    
    Так не делается.
    
    ----------------
    
    Если Вас реально интересуют лишь частоты ниже 5 Гц, то зачем Вы квантуете сигнал с частотой 44 000 Гц?
    
    Можете попробовать сначала поставить ФНЧ скажем до 100 Гц потом понижаете в цифровом виде частоту дискретизации до 500 Гц и потом считаете свой фильтр да 5 Гц

nikolz
14.10.2024 14:23
#27414908
Частота среза — это такая частота, на которой ослабление фильтра равно -3 дБ в логарифмическом масштабе (в линейном это 0,707).

Это справедливо только для монотонной АЧХ в полосе пропускания.

Для фильтра с пульсациями, например:

такое определение ошибочное.

Полагаю, что правильно будет так:

Частота среза — это такая частота, после которой (для ФНЧ) ослабление фильтра меньше заданного (например -3 дБ в логарифмическом масштабе ,в линейном - 0,707).

nikolz
14.10.2024 14:23
#27414916
Однако надо отметить, что утилита не поможет сгенерировать IIR фильтр нижних частот с малой частотой среза, например 0,0002% от Fs.

Такой крутой срез будет у фильтра очень большого порядка. Очевидно у Вас программа не может рассчитать из-за потери точности вычислений либо неправильного выбора параметров.
1. aabzel Автор
  14.10.2024 14:23
  #27414970
  Такой крутой срез будет у фильтра очень большого порядка.
  
  Утилита уже на 4-м порядке отказывается вычисления делать.
  1. nikolz
    14.10.2024 14:23
    #27415002
    Ну программа вам же сказала, что фильтр не стабильный. А нестабильный фильтр - это что? Правильно это генератор. И она вам говорит что надо делать. Вам надо делать вычисления с очень высокой точностью и считать не гладкий фильтр а с пульсациями.
    
    Но по-моему, лучше взять книжку по цифровым фильтрам там обычно есть таблицы и формулы по которым можно определить какой порядок должен быть у фильтра и с какой погрешностью надо вычислять коэффициенты.

pvvv
14.10.2024 14:23
#27415376
собственно реализация простого ФНЧ фильтра выглядит как: Y += (X-Y)*K;

при желании более высоких порядков можно просто дальше каскадировать, даже не пересчитывая коэффициенты в отдельные a[],b[], Z += (Y-Z)*K;

beefdeadbeef
14.10.2024 14:23
#27415762
Две ссылки по теме:

https://www.earlevel.com/main/2021/09/02/biquad-calculator-v3/

https://www.earlevel.com/main/2016/09/29/cascading-filters/
1. aabzel Автор
  14.10.2024 14:23
  #27418622
  А что такое Q?
  1. beefdeadbeef
    14.10.2024 14:23
    #27419248
    "Добротность". Но вообще от контекста зависит, в комментариях ко второй ссылке есть рассуждения на этот счёт.
  1. Refridgerator
    14.10.2024 14:23
    #27421316
    Q по сути значит "резонанс", вот для примера НЧ-фильтр на 1кГц с разными Q:

Refridgerator
14.10.2024 14:23
#27416654
Недостатки цифровых IIR фильтров -IIR фильтр порой усиливает сигнал.

Это не недостаток IR фильтров, это кривая программа для расчёта. Коэффициент усиления фильтра регулируется просто масштабированием всех коэффициентов знаменателя. Действительные недостатки IIR фильтров это:

а) невозможность получить линейную ФЧХ;

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

ProLimit
14.10.2024 14:23
#27422084
Как то вы совсем примитивно изложили про фильтры. Обычно для разработчика задача сложнее. Что, если в программе не нужен заранее рассчитанный фильтр, а нужен настраиваемый через параметры? WinFilter и прочие калькуляторы вам тут не помогут. В своем проекте я делал через биквадратные FIR-фильтры где в базе 2nd order s-domain function (p0 + p1*s + p2*s^2)/(q0 + q1*s + q2*s^2) и разные типы фильтров заполняют коэффициенты p, q. Так можно сделать HPF, LPF, Notch, Lag compensator, band-pass и прочие.
Далее переводим в decrete-time domain (одинаково для всех типов) для выбранной частоты сэмплирования и получаем рекуррентные соотношения. Для бОльших порядков можно каскадировать, чтобы не усложнять математику и ограничиться float (на простых микропроцессорах это оптимально)