Synthesis of a Digital IIR Low Pass Filter

Prologue

Recently I needed to synthesize a fast digital low pass filter. Moreover, this filter must work in real time on a microcontroller. Then I realized that I need to remember from which side I should approach digital IIR filters.

Sometimes it happens that they send you a screenshot of a digital filter or you see something similar to an IIR filter in the code and you need to understand what it means. A kind of reverse engineering.

Inverse problem

A first-order digital IIR filter with known coefficients a_0 and b_1 is given.

The sampling frequency F_s is also known. You need to plot the frequency response and see the cutoff frequency.

Direct task

The sampling frequency is given, the cutoff frequency Fc is given. The topology of the first order filter is given. Find its coefficients a0 and b1.

IIR filters are the case when the inverse problem is simpler than the direct problem.

Terminology

Cutoff frequency – this is the frequency at which the filter attenuation is -3 dB on a logarithmic scale (on a linear scale it is 0.707).

Solution
As for me, the easiest way to solve is just the inverse problem. Using a well-known filter, you can understand what it is capable of. This is a kind of reverse engineering.

In the filter circuit Z^(-1) this is the sample delay section. Physically, this is a register in RAM memory, with a width equal to the width of the sample. One memory cell. Each time a number passes through the delay link, it spends T_a seconds doing so. T_a is the time between samples, the inverse of the sampling frequency. Mathematically, the link looks like this.

Looking at the filter topology, we can express the following relationship between inputs and outputs.

This is the so called difference equation. The link has the property of delay

As a result of algebraic transformations, it is possible to express the relationship between the output and input of this link. We got transfer function for the diagram.

To go to the frequency domain in (4) we substitute (1). It turns out

I’m interested in the amplitude-frequency response in the range from zero to half the sampling frequency f_a (Next it will just be a mirror image).

Note the minus sign in the denominator. When calculating higher orders, all terms in the denominator will have a minus sign.

These calculations are also valid for FIR filters, since FIR filters can be considered a subset of IIR filters if the coefficients in the feedback are set to zero.

Computations

Next, I will calculate numerically using formula (6) using C as a calculator.

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;
}

The result is this graph

Here the frequency is normalized. That is, X is not the frequency, but F/F_cut. So we solved the inverse problem. Based on this filter, its cutoff frequency was estimated.

In general, any digital filter is defined by two arrays of real numbers A and 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]

Next, running the same algorithm, we get a different frequency response.

There is also such a first-order IIR low-pass filter. The good thing about it is that there is also only one variable k to configure it.

For this IIR filter, the higher k, the lower the cutoff frequency. In this case, k starts from one.

Solution of a direct problem.

As for solving the direct problem, everything is more complicated in theory. But you can use the utility WinFilter. However, it should be noted that the utility will not help generate an IIR low-pass filter with a low cutoff frequency, for example 0.0002% of Fs. At low frequencies, for some reason the utility generates not filters, but amplifiers. However, if you are lucky, the utility will calculate the necessary coefficients for you.

The coefficients of stable IIR filters are sometimes even tabular values ​​that can be found in textbooks and reference books.

Advantages of digital IIR filters:

1–IIR filters require less computation compared to FIR filters

Disadvantages of digital IIR filters

–IIR filter sometimes enhances the signal. And this is not at all necessary. Here is an example of such a filter

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

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

And its frequency response

Applications of digital low pass filters.

1–Processing of signals from sensors of physical quantities: accelerometers, gyroscopes, distances.

2–SDR processing of samples from radio receivers. Double frequency filtering at the output of a quadrature mixer.

Results

We managed to learn how to use two types of IIR low-pass filters of the first order.

There is some understanding of how to approach digital IIR filters.

As you can see, to understand digital filters, it turned out that you need to understand at least trigonometry, difference equations, complex numbers, Euler's formula and algebra.

Links

Dictionary