Design of IIR filter for hand to hand teaching series

Time:2020-11-22

[guide]: in embedded system, it is often necessary to collect analog signals. It is inevitable to introduce interference into the signal chain of collecting analog signals. So how to filter out the interference? Today, let’s describe step by step how to design and deploy an IIR filter to your system.

What is IIR filter?

Infinite impulse response(IIRInfinite Impulse RResponse) is an attribute applicable to many linear time invariant systems. These systems are characterized by an impulse response H (T), which does not completely change to zero at a specific point, but continues indefinitely. This is related to the finite impulse response(FIRFinite Impulse RIn the FIR system, for a finite T, the impulse response is exactly zero when time t > t. Common examples of linear time invariant systems are most electronic and digital filters. Systems with this property are called IIR systems or IIR filters.

This is a common textbook type strict definition of mathematics, many people see this will be confused, can you speak people?

Linear time invariant system theory, commonly known as LTI system theory, originates from applied mathematics and is directly applied in nuclear magnetic resonance spectroscopy, seismology, circuit, signal processing and control theory. It studies the response of linear, non time-varying systems to arbitrary input signals. Although the trajectories of these systems are usually measured and tracked over time (for example, acoustic waveforms), LTI systems also have trajectories in the spatial dimension when applied to image processing and field theory. Therefore, these systems are also known as linear time invariant translation, which is given in the most general range theory. In discrete (i.e., sampled) systems, the corresponding term is linear time-varying translational systems. The circuit composed of resistor, capacitor and inductor is a good example of LTI system. For example, an operational amplifier system can satisfy the time-domain superposition of signals in a certain frequency band, and input a 100Hz and 200Hz sinusoidal signal, and the output frequency is the linear superposition of these two signals.

Describe the LTI system with mathematics

linear: input\(x_1(t)\)To generate a response\(y_1(t)\), and input\(x_2(t)\)To generate a response\(y_2(t)\)Then zoom and add and input\(a_1x_1(t)+a_2x_2(t)\)To generate a scaling, plus and response\(a_1y_1(t)+a_2y_2(t)\)In which\(a_1\)and\(a_2\)Is a scalar. For any, there is: input\(\sum_{k=0}^{N}a_kx_k(t)\)To generate a response\(\sum_{k=0}^{N}a_ky_k(t)\)

Time invariance : if the input signal of the system is delayed\(\tau\)Second, then the output response is delayed accordingly\(\tau\)Seconds. Describe it mathematically, i.e. if input\(x_1(t)\)To generate a response\(y_1(t)\), and input\(x_1(t+\tau)\)To generate a response\(y_1(t+\tau)\)。 This description is still not easy to understand. Let’s take a picture. There are pictures with truth

It is assumed that a signal amplifying circuit can amplify 100Hz sinusoidal signal twice

image-20200326205428645

Then the output is:

image-20200326205811819

For 200Hz sinusoidal signal, the amplification factor is 1.7 times. (those who have done op amp circuit design should have experience. In the same frequency band, the amplification factor is often uneven, that is, the amplitude frequency response is not flat in the frequency band, which is relatively common).

That is, the input is:

image-20200326210319201

The response was:

image-20200326210332571

If 100Hz and 200Hz time-domain superimposed signals are input, the input is:

image-20200326210550908

Then the response is as follows:

image-20200326210603870

It can be seen from these figures that the shape of the input signal remains unchanged and the output is a linear time-domain superposition of the corresponding input.

For time invariantIn a real circuit, if the input is delayed for a certain time, the response corresponds to the output with the same delay.

So much of the above text is just to describe where the IIR filter can be used for digital filtering of signals. In summary, it is suitable for linear time invariant systems. In other words, in most circuit systems, we can try to use IIR filter for digital filtering.

So what is IIR filter? From the books on digital signal processing, we can see such a Z-transform signal flow diagram:

image-20200326221206401

image-20200326214242009In the digital system, it means the last sampling value for the input signal and the last output value for the output.

In the time domain, the flow diagram described in time domain is as follows:

image-20200326221258617

If you are familiar with the z-transform, the z-transform transfer function is:

image-20200326221348262

From the programming point of view, if x (n-1) represents the last signal, it may be the last sampling from ADC, while y (n-1) is the output value of the last filter, which is better understood. X (N-N) represents the signal of the nth input sample before, and Y (n-m) is the output of the previous m-th filter.

I said so much just to better understand the concept, onlyOnly when the concept is understood correctly can we use it to fight for it。 Concept understanding is very important for engineers.

How to design?

Open FDATool

Matlab provides a very easy-to-use FDATool to help us design digital filters. The really wonderful place begins. Let’s wait and see how to design and implement an IIR filter step by step.

First, open MATLAB and type FDATool in the command line

image-20200326233121867

The pop-up window is FDATool, as follows:

In the specific design, several related concepts need to be clarified

Fs: sampling rate, unit: Hz. If it is deployed in the system, make sure that the sample is sampled according to the constant sampling rate, otherwise the desired effect will not be obtained.

Fpass: passband, unit: Hz, i.e. the highest frequency expected to pass through the system.

Fstop: cutoff frequency, that is, the frequency at – 3dB of amplitude frequency response. If you don’t understand this, please refer to relevant books.

DB: This is a term for the output and input multiples of a unitless response. In electricity, the conversion relationship between decibel and magnification is as follows:

  • A (V) (DB) = 20lg (VO / VI); voltage gain, VO is output voltage, VI is input voltage
  • A (I) (DB) = 20lg (IO / II); current gain, IO is output current, II is input current
  • A (P) (DB) = 10lg (PO / PI); power gain, Po is output power, PI is input power

Filter typeThere are Butterworth (Butterworth), Chebyshev type I, Chebyshev type II, (Chebyshev) and elipic.

image-20200327002802811

  • ButterworthAlso known as the maximum flat filter. The characteristic of Butterworth filter is that the frequency response curve in the passband is flat as far as possible without ripple.
  • ChebyshevIt is a filter with ripple fluctuation of frequency response amplitude on passband or stopband. The attenuation of Chebyshev filter is faster than that of Butterworth filter in transition band, but the amplitude frequency characteristic of frequency response is not as flat as that of Butterworth filter.
  • EllipseElliptic filter is a kind of filter with ripple in passband and stopband.
  • I will not introduce them one by one. If you are interested, you can check the signal processing books.

In terms of their characteristics, several of them are briefly introduced here

  • Butterworth has the flattest passband.
  • The elliptic filter has the fastest attenuation, but the passband and stopband have ripple.
  • The attenuation of Chebyshev filter is faster than that of Butterworth filter, but slower than that of elliptic filter, and the ripple area can be selected.

Suppose we need to design an IIR filter, the sampling rate is 32000hz, the useful signal frequency is within 2000Hz, design IIR filter to digital filter the signal. In order to save computational power, we specify the order of the filter, that is, the maximum value of N / m in the transfer function. Generally speaking, n is greater than m.

Here, the order is 8, the type is Butterworth IIR filter, the input order is 8, the sampling rate is 32000hz, and the cutoff frequency is set to 2000Hz, and then click design filter, as shown in the following figure:

The phase frequency response curve is as follows:

Besides the songs of Chu, we can also put the amplitude frequency and phase frequency curves on a frequency coordinate to see the design results

The filter parameters are derived. Here we choose,

Select ASCII

Then you get a file and save 2KHz_ LPF.fcf , file name as you like.

The contents of the document are as follows:

% Generated by MATLAB(R) 8.4 and the Signal Processing Toolbox 6.22.
% Generated on: 27-Mar-2020 21:27:06

% Coefficient Format: Decimal

% Discrete-Time IIR Filter (real)                            
% -------------------------------                            
% Filter Structure    : Direct-Form II, Second-Order Sections
% Number of Sections  : 4                                    
% Stable              : Yes                                  
% Linear Phase        : No                                   

                                                            
SOS Matrix:                                                  
1  2  1  1  -1.7193929141691948  0.8610574795347461          
1  2  1  1  -1.5237898734101736  0.64933827386370635         
1  2  1  1  -1.4017399331200424  0.51723237044751591         
1  2  1  1  -1.3435020629061745  0.45419615396638446         
                                                             
Scale Values:                                                
0.035416141341387819                                         
0.031387100113383172                                         
0.028873109331868367                                         
0.027673522765052503

At this point, the design work is finished, and the filter deployment test phase is immediately entered.

Deploy test filter

Here, inexperienced friends may say, how can I use such a pile of parameters?

Do you need to write the formula described above? Of course, you can also do this, not here, arm’s CMSIS library has helped you design a wide range of digital signal processing functions, and after testing, you can use it directly here. It’s not difficult to write by yourself if you are interested. As long as you understand the concept of Z-transfer function, it’s very easy to realize. Here we use 32-bit floating-point implementation function: arm_ biquad_ cascade_ df1_ f32。

This function is located in CMSIS / DSP / source / filteringfunctions / arm_ biquad_ cascade_ df1_ init_ F32. C

And CMSIS, DSP, source, filteringfunctions, arm_ biquad_ cascade_ df1_ f32.c

Let’s take a look at this function

arm_biquad_cascade_df1_init_f32.c:

/*
*Function: initialize the filter
*S: points to an instance of the floating-point SOS concatenation structure.
*Numstages: the number of second-order SOS in the filter
*Pcoeffs: filter parameter pointer, the parameters are stored in the following order
*          {b10, b11, b12, a11, a12, b20, b21, b22, a21, a22, ...}
*Pstate: history state buffer pointer
*/
void arm_biquad_cascade_df1_init_f32(
        arm_biquad_casd_df1_inst_f32 * S,
        uint8_t numStages,
  const float32_t * pCoeffs,
        float32_t * pState)
{
  /* Assign filter stages */
  S->numStages = numStages;

  /* Assign coefficient pointer */
  S->pCoeffs = pCoeffs;

  /* Clear state buffer and size is always 4 * numStages */
  memset(pState, 0, (4U * (uint32_t) numStages) * sizeof(float32_t));

  /* Assign state pointer */
  S->pState = pState;
}

arm_ Math. H defines the structure to be used. For this example, the related structure is arm_ biquad_ casd_ df1_ inst_ F32

typedef struct
{
  The number of nodes of order / * 2 should be 2 * numstages*/
  Float * pstate; / * state coefficient array pointer, the array length is 4 * numstages*/
  Float * pcoeffs; / * coefficient array pointer, the length of the array is 5 * numstages*/
} arm_biquad_casd_df1_inst_f32;

The filter function is arm_ biquad_ cascade_ df1_ F32

/**
 ** s: points to an instance of the floating-point biquad concatenation structure
 ** PSRC: points to the input data block.
 ** PDST: points to the output data block.
 *Blocksize: the number of samples to process per call.
 *Return value: none
 */
void arm_biquad_cascade_df1_f32(
  const arm_biquad_casd_df1_inst_f32 * S,
  float * pSrc,
  float * pDst,
  unsigned int blockSize)
{
  Float * pin = PSRC; / * source pointer*/
  Float * pout = PDST; / * destination pointer*/
  Float * pstate = s - > pstate; / * state pointer*/
  Float * pcoeffs = s - > pcoeffs; / * parameter pointer*/
  Float ACC; / * accumulator*/
  Float B0, B1, B2, A1, A2; / * filter parameters*/
  Float xN1, xn2, yn1, YN2; / * filter state variables*/
  Float xn; / * temporary input*/
  Unsigned int sample, stage = s - > numstages; / * cycle count*/

  do
  {
    /* Reading the coefficients */
    b0 = *pCoeffs++;
    b1 = *pCoeffs++;
    b2 = *pCoeffs++;
    a1 = *pCoeffs++;
    a2 = *pCoeffs++;

    Xn1 = pState[0];
    Xn2 = pState[1];
    Yn1 = pState[2];
    Yn2 = pState[3];

    sample = blockSize >> 2u;

    while(sample > 0u)
    {
      /*Read the first input*/
      Xn = *pIn++;

      /* acc =  b0 * x[n] + b1 * x[n-1] + b2 * x[n-2] + a1 * y[n-1] + a2 * y[n-2] */
      Yn2 = (b0 * Xn) + (b1 * Xn1) + (b2 * Xn2) + (a1 * Yn1) + (a2 * Yn2);

      /* Store the result in the accumulator in the destination buffer. */
      *pOut++ = Yn2;

      /*The status should be updated after each output calculation*/
      /*The status should be updated to:*/
      /* Xn2 = Xn1    */
      /* Xn1 = Xn     */
      /* Yn2 = Yn1    */
      /* Yn1 = acc   */

      /* Read the second input */
      Xn2 = *pIn++;

      /* acc =  b0 * x[n] + b1 * x[n-1] + b2 * x[n-2] + a1 * y[n-1] + a2 * y[n-2] */
      Yn1 = (b0 * Xn2) + (b1 * Xn) + (b2 * Xn1) + (a1 * Yn2) + (a2 * Yn1);

      /*Store the result in the accumulator of the destination buffer*/
      *pOut++ = Yn1;

      /*The status should be updated after each output calculation*/
      /*The status should be updated to:*/
      /* Xn2 = Xn1    */
      /* Xn1 = Xn     */
      /* Yn2 = Yn1    */
      /* Yn1 = acc   */

      /*Read the third input*/
      Xn1 = *pIn++;

      /* acc =  b0 * x[n] + b1 * x[n-1] + b2 * x[n-2] + a1 * y[n-1] + a2 * y[n-2] */
      Yn2 = (b0 * Xn1) + (b1 * Xn2) + (b2 * Xn) + (a1 * Yn1) + (a2 * Yn2);

      /*Store the result in the accumulator of the destination buffer*/
      *pOut++ = Yn2;

      /*The status should be updated after each output calculation*/
      /*The status should be updated to:*/
      /* Xn2 = Xn1    */
      /* Xn1 = Xn     */
      /* Yn2 = Yn1    */
      /* Yn1 = acc   */
      /*Fourth input*/
      Xn = *pIn++;

      /* acc =  b0 * x[n] + b1 * x[n-1] + b2 * x[n-2] + a1 * y[n-1] + a2 * y[n-2] */
      Yn1 = (b0 * Xn) + (b1 * Xn1) + (b2 * Xn2) + (a1 * Yn2) + (a2 * Yn1);

      /*Store the result in the accumulator of the destination buffer*/
      *pOut++ = Yn1;

      /*The status should be updated after each output calculation*/
      /*The status should be updated to:*/
      /* Xn2 = Xn1    */
      /* Xn1 = Xn     */
      /* Yn2 = Yn1    */
      /* Yn1 = acc   */
      Xn2 = Xn1;
      Xn1 = Xn;

      /*Decrement cycle counter*/
      sample--;
    }

    /*If blocksize is not a multiple of 4,
    *Please calculate any remaining output samples here.
    *Loop unrolling is not used*/
    sample = blockSize & 0x3u;

    while(sample > 0u)
    {
      /*Read input*/
      Xn = *pIn++;

      /* acc =  b0 * x[n] + b1 * x[n-1] + b2 * x[n-2] + a1 * y[n-1] + a2 * y[n-2] */
      acc = (b0 * Xn) + (b1 * Xn1) + (b2 * Xn2) + (a1 * Yn1) + (a2 * Yn2);

      /*Store the result in the accumulator of the destination buffer*/
      *pOut++ = acc;

      /*The status should be updated each time the output is calculated. * /
      /*The status should be updated to:*/
      /* Xn2 = Xn1    */
      /* Xn1 = Xn     */
      /* Yn2 = Yn1    */
      /* Yn1 = acc   */
      Xn2 = Xn1;
      Xn1 = Xn;
      Yn2 = Yn1;
      Yn1 = acc;

      /*D decrement cycle counter*/
      sample--;
    }

    /*Store the updated state variables back to the pstate array*/
    *pState++ = Xn1;
    *pState++ = Xn2;
    *pState++ = Yn1;
    *pState++ = Yn2;

    /*The first stage is from input buffer to output buffer*/
    /*Subsequent numstages occur in place in the output buffer*/
    pIn = pDst;

    /*Reset output pointer*/
    pOut = pDst;

    /*Decrement cycle counter*/
    stage--;

  } while(stage > 0u);
}

Start test:

#include 
#include 
/*
SOS Matrix:                                                  
1  2  1  1  -1.7193929141691948  0.8610574795347461          
1  2  1  1  -1.5237898734101736  0.64933827386370635         
1  2  1  1  -1.4017399331200424  0.51723237044751591         
1  2  1  1  -1.3435020629061745  0.45419615396638446         
                                                             
Scale Values:                                                
0.035416141341387819                                         
0.031387100113383172                                         
0.028873109331868367                                         
0.027673522765052503  
Do the following conversion:
1. Zoom
[1  2  1] * 0.035416141341387819
[1  2  1] * 0.031387100113383172
[1  2  1] * 0.028873109331868367
[1  2  1] * 0.027673522765052503
The results are as follows:
[0.035416141341387819  2*0.035416141341387819  0.035416141341387819]
[0.031387100113383172  2*0.031387100113383172  0.031387100113383172] 
[0.028873109331868367  2*0.028873109331868367  0.028873109331868367] 
[0.027673522765052503  2*0.027673522765052503  0.027673522765052503]
2. Omit the fourth column parameters
3. Multiply the last two columns by - 1, namely:
0.035416141341387819  2*0.035416141341387819  0.035416141341387819  -1.7193929141691948  0.8610574795347461          
0.031387100113383172  2*0.031387100113383172  0.031387100113383172  -1.5237898734101736  0.64933827386370635         
0.028873109331868367  2*0.028873109331868367  0.028873109331868367  -1.4017399331200424  0.51723237044751591         
0.027673522765052503  2*0.027673522765052503  0.027673522765052503  -1.3435020629061745  0.45419615396638446 
So we get the filter array
*/
#define IIR_ Section 4 / * see the previous design output is 4 SOS blocks*/
static float iir_ state[4*IIR_ Section]; / * historical state buffer*/
const float iir_coeffs[5*IIR_SECTION]={
    0.035416141341387819,2*0.035416141341387819,0.035416141341387819,1.7193929141691948,-0.8610574795347461,    0.031387100113383172,2*0.031387100113383172,0.031387100113383172,1.5237898734101736,-0.64933827386370635,    0.028873109331868367,2*0.028873109331868367,0.028873109331868367,1.4017399331200424,-0.51723237044751591,    0.027673522765052503,2*0.027673522765052503,0.027673522765052503,1.3435020629061745,-0.45419615396638446
};
static arm_biquad_casd_df1_inst_f32 S;
/*It is assumed that 512 points are sampled*/
#define BUF_SIZE 512
#define PI 3.1415926
#define SAMPLE_RATE  32000 /*32000Hz*/
int main()
{
    float raw[BUF_SIZE];
    float raw_4k[BUF_SIZE];
    float raw_out[BUF_SIZE];

    float raw_noise[BUF_SIZE];
    float raw_noise_out[BUF_SIZE];

    arm_biquad_casd_df1_inst_f32 S;
    FILE *pFile=fopen("./simulation.csv","wt+");
    if(pFile==NULL)
    {
        printf("file opened failed");
        return -1;
    }

    for(int i=0;i

Use the CSV file to store the simulation data, open it directly with Excel, and generate the curve chart of the line data as follows:

image-20200517125200485

  • The first picture shows the waveform of 800Hz signal mixed with random noise
  • The second picture shows 4000Hz signal, which is useless interference signal for the assumed system
  • The third picture shows that the frequency of useful signal has been well restored after filtering with 800Hz random noise
  • The fourth picture shows the waveform of 800Hz signal mixed with random noise and 4000Hz interference. For the system, from the time domain, it is obvious that the useful signal has been completely distorted
  • In the fifth picture, 800Hz signal is mixed with random noise, and 4000Hz interference input is superimposed. After passing through the low-pass filter, the waveform is basically the same as the third picture, and the interference signal has been filtered out very well.

summary

  • IIR filter in linear time invariant system can solve the general noise problem in engineering

  • If it is necessary to design a band-pass filter and a high pass filter, the steps are basically similar, but the parameters of the filter and the number of SOS blocks may be different

  • If the system has strict requirements on phase frequency response, other digital filter topologies should be used

  • In practical application, if the order is not high, the powerful MCU or DSP can use floating point processing directly.

  • If there are strict real-time requirements for processing speed and need to filter in a very short time, we can consider reducing the order or using fixed-point IIR filtering algorithm. Or we can optimize the function at assembly level.

The article is from WeChat official account: embedded Inn, focusing on official account, and getting more updates.

Design of IIR filter for hand to hand teaching series

Recommended Today

Interviewer: young man, what do you think of the principle of distributed system

1 Concept 1.1 model 1.2 copies 1.3 indicators for measuring distributed systems 2. Principle of distributed system 2.1 data distribution 2.2 basic copy agreement 2.3 lease mechanism 2.4 quorum mechanism 2.5 log technology 2.6 two phase submission protocol 2.7 MVCC 2.8 Paxos protocol 2.9 CAP 1 Concept 1.1 model node In a specific project, a […]