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

Пролог

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

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

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

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

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

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

Дана частота дискретизации, дана частота среза Fс. Дана топология фильтра первого порядка. Найти его коэффициенты 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)$

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

$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; nsize; 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; nsize; 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 фильтров нижних частот первого порядка.

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

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

Ссылки

Словарь

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