Design and Implementation of High-Order IIR filter

When I first left school, many of the things I learned at school could not be used at work. Over time, when I needed these things, I found that all the knowledge I had learned had been returned to the teacher. For me, the design of digital filter is an example.

 

Some time ago, it took about two months to review the original digital signal processing content, from the Euler's formula to the transfer function. The following are my work achievements over the past one or two months. ^_^

 

Digital IIR filter designed by analog filter: Impulse Response invariant method and bilinear transform method

 

I mainly focus on bilinear transformation.

 

When selecting a simulated filter prototype, the following options are available: Butterworth/Chebyshev I/Chebyshev II/elliptic /...

 

N-order iir lpf:

 

First, we need to find the Amplitude Square Function of the simulated filter prototype. This can be easily found on Wiki. Don't ask me how I got here. I don't know either.

 

Take Chebyshev I as an example:

* 1
* | H (JW) | ^ 2 = --------------------------
* SQRT (1 + e ^ 2 * tn ^ 2 (w/WC ))

Because S = JW

 

* 1
* H (s) * H (-S) = -----------------------
* 1 + e ^ 2 * tn ^ 2 (-J * s/WC)

 

The above formula can be used to obtain the pole of the simulated filter. Don't ask me how to solve it. Wiki has

 

* PK = (+-sinh (a) * sin (B) + jcosh (a) * Cos (B) * WC
*
* = RP + J * IP
*
* 1 1
* Where a = --- * arsinh (---),
* N e
*
* Pi 2 * K-1
* B = --- * -------, k = 1, 2,... n/2
* 2 N

 

Of course, only the n/2 poles shown in the above formula are on the real axis, and the other half of the Poles are under the real axis, and the two are bounded. For the convenience of writing comments in the code, I wrote the PK ~

 

By now, all the parts related to the analog filter type have been processed. The following processing applies to various analog filter prototypes:

 

The key to implementing a digital filter in the code is to convert a high-order IIR filter into a first-or second-order filter, here we can use the addition and multiplication of the bounded complex numbers to remove the features of the virtual part and convert the n-order IIR filter to N/2 2.

 

Continue to calculate

 

* N/2 1
* H (S) = Mul -----------------------
* K = 1 (1-S/PK) * (1-S/PK ~)
*
* N/2 1
* = Mul ------------------------------------------
* K = 1 1-(PK + PK ~) /(Pk * PK ~)) * S + s ^ 2/(Pk * PK ~)
*
* N/2 1
* = Mul ---------------------------------------------
* K = 1 1-(2 * RP)/(RP ^ 2 + IP ^ 2) * s + s ^ 2/(RP ^ 2 + IP ^ 2)
*
* N/2 1
* = Mul -------------------
* K = 1 1 + A1 * s + A2 * s ^ 2

 

Now it's easy. The next step is to convert the analog filter into a digital filter with a dual-line variation:

* N/2 1 + 2 * Z ^-1 + Z ^-2
* H (z) = Mul --------------------------------------------
* K = 1 (1 + A1 + A2) + (2-2 * A2) * Z ^-1 + (1-a1 + A2) * Z ^-2

 

The following code is easy to implement.

 

LPF is OK, and the frequency conversion method for HPF is also very easy. BPF is a little troublesome and requires a level 4 knot.

 

======================================

 

The following are some concepts:

 

Natural Frequency: Hz

Normalized digital frequency: Natural Frequency/sampling rate

Digital Frequency: 2 * pI * normalized digital frequency

Analog frequency: Tan (pI * normalized digital frequency)

 

Let's guess: What frequency is the WC above?

 

==================================

The last is a little bit of experience:

 

During formula derivation, multiple variables are used to replace the results of a specific stage, so that the code can be clearly written. In addition, if an error is found in formula derivation, the entire code will not be overwritten. Of course, it is another question when optimizing performance.