Design of IIR filter for hand-held teaching series

[guide reading]: it is often necessary to collect analog signals in embedded systems. It is inevitable to introduce interference into the signal chain of collecting analog signals, so how to filter the interference? Let's take a step-by-step description of how to design and deploy an IIR filter to your system.

What is an IIR filter?

IIR (Infinite Impulse Response) is an attribute applicable to many linear time invariant systems. The characteristic of these systems is that they have an impulse response h (T). The impulse response h (T) will not completely change to zero at a specific point, but will continue indefinitely. This is in contrast to the Finite Impulse Response system. In the Finite Impulse Response system, for a finite T, the impulse response does exactly turn to zero when time t > t. Common examples of linear time invariant systems are most electronic and digital filters. The system with this property is called IIR system or IIR filter.

This is a common textbook mathematical rigorous definition. Many people will be blinded when they see it. Can they speak human words?

Linear time invariant system theory, commonly known as LTI system theory, originates from applied mathematics and is directly applied in the fields of NMR spectroscopy, seismology, circuit, signal processing and control theory. It studies the response of a linear, time-varying system to any input signal. Although the trajectories of these systems are usually measured and tracked with time (such as acoustic waveform), when applied to image processing and field theory, LTI systems also have trajectories in spatial dimension. Therefore, these systems are also called linear time-varying translation, which is given in the most general range theory. In the discrete (i.e. sampling) system, the corresponding term is linear non time varying translation system. The circuit composed of resistance, capacitance and inductance is a good example of LTI system. For example, an operational amplifier system satisfies the time-domain superposition of signals in a certain frequency band, inputs a 100Hz and 200Hz sine signal, and the output frequency is the linear superposition of the two signals.

Describe the LTI system mathematically:

Sum {k = 0} ^ {n} a {K (T) \), generating a response \ (\ sum {k = 0} ^ {n} a {K (T) \)

Time invariance: if the input signal of the system is delayed by \ (\ tau \) seconds, the output response is also delayed by \ (\ tau \) seconds. Use mathematical description, that is, if you input \ (x_1(t)), the response \ (y_1(t)), and input \ (x_1(t+\tau)), the response \ (y_1(t+\tau) \). It's still hard to understand. Here's a picture. There's a picture and a truth:

It is assumed that a signal amplification circuit can amplify 100Hz sine signal twice:

The output is:

For 200 Hz sinusoidal signal, the amplification is 1.7 times. (friends who have done operational amplifier circuit design should have experience. In the same frequency band, the amplification factor is often not flat, that is, the amplitude frequency response is not flat in the frequency band, which is quite common).

That is, the input is:

The response is:

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

The response is:

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

In a real circuit, if the input is delayed for a certain time, the response will be delayed for the same time.

There are so many words above just to describe where the IIR filter can be used for digital filtering of signals. In summary, it is applicable to linear time invariant systems. In other words, in most circuit systems we can try to use IIR filter for digital filtering.

So what exactly is IIR filter? From the books of digital signal processing, we can see the Z-transform signal flow chart as follows:

It means one beat delay. In 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, for the above flow chart, the description in the time domain is as follows:

If you are familiar with Z transformation, the transfer function of Z transformation is:

From the programming point of view, the above digital filter, x(n-1), represents the last signal, which may be the last sampling from ADC, and y(n-1) is the output value of the last filter, so it is easy to understand that x(n-N) represents the previous nth input sample signal number, and y(n-M) is the output of the previous nth filter.

I have said so much just to better understand the concept. Only when the concept is understood correctly can we strive for it. Conceptual understanding is very important for engineers.

How to design?

Open fdatool

MATLAB provides FDATool, which is very easy to use, to help us design digital filter. The really wonderful place is beginning. Let's wait and see how to design and implement an IIR filter step by step.

First, open MATLAB and type fdatool on the command line

The pop-up window is fdatool, as follows:

In the design, there are several related concepts that need to be clarified:

Fs: the sampling rate, in Hz, is actually deployed in the system. Please make sure that the samples are sampled according to the constant sampling rate, otherwise the desired effect will not be obtained.

Fpass: passband, unit: Hz, that is, the highest frequency expected to pass in the medium term of the system.

Fstop: as of frequency, i.e. frequency at - 3dB of amplitude frequency response. If you don't understand this, please refer to relevant books by yourself.

dB: This is a term for the multiple of output and input without unit response. In electricity, the conversion relationship between decibel and magnification is as follows:

  • A(V)(dB)=20lg(Vo/Vi); voltage gain, Vo is the output voltage, Vi is the 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 type: Butterworth, Chebyshev Type I,Chebyshev Type II, (Chebyshev), Elipic, etc. are available.

  • Butterworth, also known as the maximum flat filter. The characteristic of Butterworth filter is that the frequency response curve in the passband is flat as much as possible without ripple.
  • Chebyshev is a kind of filter with ripple like frequency response amplitude on the passband or stopband. Chebyshev filter decays faster than Butterworth filter in the transition band, but the amplitude frequency characteristic of frequency response is not as flat as the latter.
  • Elliptic filter is a kind of filter with ripple of passband and stopband.
  • ... I will not introduce them here. If you are interested, you can check the signal processing books.

In terms of its characteristics, several of them are briefly introduced here:

  • Butterworth has the flattest passband.
  • Elliptical filter has the fastest attenuation, but there are ripples in the passband and stopband.
  • The attenuation of Chebyshev filter is faster than Butterworth filter, but slower than elliptical 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 carry on the digital filtering to 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 specified order is 8, the type is Butterworth IIR filter, the input order is 8, the sampling rate is 32000Hz, the cutoff frequency is set to 2000Hz, and then click Design Filter as shown below:

The phase frequency response curve is as follows:

Besides Chuci, we can also put the amplitude frequency and phase frequency curves on a frequency coordinate to see the design results:

To export the filter parameters, here we choose,

Select ASCII

Then you get a file. Save 2kHz_LPF.fcf. The file name is 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 stage is entered immediately.

Deploy test filter

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

Do you need to write the calculation formula described above? Of course, you can do the same. I won't write here. ARM's CMSIS library has helped you design a wide variety of digital signal processing functions, and after testing, you can use them directly here. It's not difficult to write by yourself. 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 filter
*S        :Points to an instance of a floating-point SOS cascade structure.
*numStages:Number of second order SOS in 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 status 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
{
  unsigned int numStages; /*2 The number of steps should be 2*numStages        */
  float *pState;          /*State coefficient array pointer, 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 specific filter function of the filter is arm ﹣ biquad ﹣ cascade ﹣ DF1 ﹣ F32

/**
 *  *S       : Points to an instance of the floating-point Biquad concatenation structure
 *  *pSrc    : Points to the input block.
 *  *pDst    : Points to the output 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;                 /*Status 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 variable*/
  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 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 every time the output is calculated */
      /* 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 every time the output is calculated */
      /* 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 every time the output is calculated */
      /* The status should be updated to: */
      /* Xn2 = Xn1    */
      /* Xn1 = Xn     */
      /* Yn2 = Yn1    */
      /* Yn1 = acc   */
      /* Read the 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 every time the output is calculated */
      /* 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 expansion 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 variable back into the pState array */
    *pState++ = Xn1;
    *pState++ = Xn2;
    *pState++ = Yn1;
    *pState++ = Yn2;

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

    /* Reset output pointer */
    pOut = pDst;

    /* Decrement cycle counter */
    stage--;

  } while(stage > 0u);
}

Start test:

#include <stdio.h>
#include <math.h>
/*
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
 Get:
[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.Round off the fourth column of parameters
3.Multiply the last two columns by - 1, that is:
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 array of filter systems
*/
#Define IIR? Section 4 / * see the previous design output as 4 SOS blocks*/
static float iir_state[4*IIR_SECTION];/*History status 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;
/*Assume 512 sampling points*/
#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<BUF_SIZE;i++)
    {
    /*Analog 9000Hz input signal         */
        raw[i] = 0.5*1024.0/3*sin(2*PI*800*i/32000.0f)+rand()%50;
        raw_4k[i] = 0.5*1024.0/3*sin(2*PI*4000*i/32000.0f);
    /*Analog 9000Hz +11000Hz superimposed input*/
        raw_noise[i] = raw[i] + raw_4k[i];
    }

    arm_biquad_cascade_df1_init_f32(&S, IIR_SECTION, (float *)&iir_coeffs[0], (float *)&iir_state[0]);
    arm_biquad_cascade_df1_f32(&S, raw, raw_out, BUF_SIZE);

    for(int i=0;i<BUF_SIZE;i++)
    {
       fprintf(pFile,"%f,",raw[i]);
    }

    fprintf(pFile,"\n");
    for(int i=0;i<BUF_SIZE;i++)
    {
        fprintf(pFile,"%f,",raw_4k[i]);
    }
    fprintf(pFile,"\n");
    for(int i=0;i<BUF_SIZE;i++)
    {
        fprintf(pFile,"%f,",raw_out[i]);
    }

    /*Reinitialize*/
    arm_biquad_cascade_df1_init_f32(&S, IIR_SECTION, (float *)&iir_coeffs[0], (float *)&iir_state[0]);
    arm_biquad_cascade_df1_f32(&S, raw_noise, raw_noise_out, BUF_SIZE);

    fprintf(pFile,"\n");
    for(int i=0;i<BUF_SIZE;i++)
    {
        fprintf(pFile,"%f,",raw_noise[i]);
    }

    fprintf(pFile,"\n");
    for(int i=0;i<BUF_SIZE;i++)
    {
        fprintf(pFile,"%f,",raw_noise_out[i]);
    }

    fclose(pFile);

    return 0;
}

Using the csv file, the simulation data is stored and directly opened with excel, and the line data is generated into the curve chart as follows:

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

Summary:

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

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

  • When it needs to be reminded, the phase frequency response of IIR is not linear. If the system has strict requirements for phase frequency response, it needs to adopt other digital filter topology forms

  • In practical application, if the order is not high, the single-chip microcomputer or DSP with strong computing power can directly use floating-point processing.

  • If there is a strict real-time requirement for processing speed, it is necessary to filter in a very short time. It can be considered to reduce the order, or to use fixed-point IIR filtering algorithm. Or the functions in this paper can be optimized at assembly level.

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

Tags: C MATLAB Attribute Programming ascii

Posted on Sat, 16 May 2020 23:28:59 -0700 by Derek