Implementation of IIR Filter (C + +)

Iirthe implementation of the filter (C + +)

A program that is being written recently requires IIR filters, and The coefficients of IIR filters need to be adjusted dynamically. So it took a little time to study The implementation of IIR filter.

The parameters used in the IIR filter are determined beforehand, there is a website, as long as the filter parameter characteristics of the input, directly can generate the required C code.

Http://www-users.cs.york.ac.uk/~fisher/mkfilter/trad.html

have been lazy to use the results of this site directly, so there is no accumulation of useful code on the hand. This time, you need to start from scratch.

There are two issues that I face:

1. Based on the parameters of the filter (type of filter, cutoff frequency, order of filter, etc.), the coefficients of the difference equation corresponding to the filter are calculated.

2. Use the coefficients of the resulting difference equations to construct a working filter.

The first problem, for different types of filters, such as Butterworth type,Bessel type, and so on, the calculation method of the filter coefficients are different. I haven't done this part of the work yet, so I'm going to write an article after I have implemented the calculation methods of the coefficients of several common filter types.

Let's write the second question here. IIR the difference equation for the filter corresponds to:

The corresponding system functions are:

Here default a[0] = 1. In fact, you can always adjust the value of A[k] and B[k] to make a[0] = 1, so this condition always satisfies.

According to Oppenheimer's "discrete-time signal processing" above,IIR Filter has two basic realization forms, which become direct I and direct type Ii. I have written two classes to achieve both forms.

DirectItype

Class Iir_i{private:    double *m_pnum;    Double *m_pden;    Double *m_px;    Double *m_py;    int m_num_order;    int m_den_order;public:    iir_i ();    void Reset ();    void Setpara (double num[], int num_order, double den[], int den_order);    void resp (double data_in[], int m, double data_out[], int n);    Double filter (double data);    void filter (double data[], int len);    void filter (double data_in[], double data_out[], int len);};

whichm_pxStoreX[n-k]The value (M_px[0]Storex[n-0],M_px[1]StoreX[n-1], etc.),m_pyStoreY[n-k]The value (M_py[0]Storey[n-0],M_py[1]StoreY[n-1], etc.).

Three filter functions are used to do the actual filtering operation. Prior to this, the coefficients of the filter need to be initialized with the Setpara function.

Here is the implementation code:

/** \brief The internal state of the filter is zeroed, the coefficient of the filter remains * \return */void Iir_i::reset () {for (int I = 0; I <= m_num_order; i++) {m    _pnum[i] = 0.0;    } for (int i = 0; I <= m_den_order; i++) {m_pden[i] = 0.0;    }}iir_i::iir_i () {m_pnum = NULL;    M_pden = NULL;    M_PX = NULL;    M_py = NULL;    M_num_order =-1; M_den_order = -1;};/ * * \brief * * \param num molecular polynomial coefficients, ascending order, num[0] is constant term * \param m molecular polynomial order * \param den denominator polynomial coefficients, ascending order, den[0] is constant term * \param m min    Order of the female polynomial * \return */void Iir_i::setpara (double num[], int num_order, double den[], int den_order) {delete[] m_pnum;    Delete[] M_pden;    Delete[] m_px;    Delete[] m_py;    M_pnum = new Double[num_order + 1];    M_pden = new Double[den_order + 1];    M_num_order = Num_order;    M_den_order = Den_order;    M_PX = new Double[num_order + 1];    M_py = new Double[den_order + 1];        for (int i = 0; I <= m_num_order; i++) {m_pnum[i] = Num[i];    M_px[i] = 0.0; } for (int i = 0; I <= m_den_order; i++)    {M_pden[i] = Den[i];    M_py[i] = 0.0; }}/** \brief filter function, using direct I-type structure * * \param data incoming input data * \return filtered results */double iir_i::filter (double data) {m_py[0] = 0.0;    The result of storing filter m_px[0] = data;    for (int i = 0; I <= m_num_order; i++) {m_py[0] = M_py[0] + m_pnum[i] * M_px[i];    } for (int i = 1; I <= m_den_order; i++) {m_py[0] = m_py[0]-m_pden[i] * m_py[i];    } for (int i = M_num_order; I >= 1; i--) {m_px[i] = m_px[i-1];    } for (int i = M_den_order; I >= 1; i--) {m_py[i] = m_py[i-1]; } return m_py[0];} /** \brief filter function, using direct I-type structure * * \param data[] Incoming input data, return the result of the filter * \param len data[] array length * \return */void iir_i::filter (dou    ble data[], int len) {int i, K;        for (k = 0; k < Len; k++) {m_px[0] = data[k];        DATA[K] = 0.0;        for (i = 0; I <= m_num_order; i++) {data[k] = Data[k] + m_pnum[i] * M_px[i];    } for (i = 1; I <= m_den_order; i++)    {Data[k] = Data[k]-m_pden[i] * m_py[i];        }//We get the Y value now m_py[0] = data[k];        for (i = m_num_order; I >= 1; i--) {m_px[i] = m_px[i-1];        } for (i = m_den_order; I >= 1; i--) {m_py[i] = m_py[i-1];  }}}/** \brief filter function, using direct I-structure * * \param data_in[] Input data * \param data_out[] Save filtered data * \param len array length * \return */void    Iir_i::filter (double data_in[], double data_out[], int len) {int I, K;        for (k = 0; k < Len; k++) {m_px[0] = data_in[k];        M_py[0] = 0.0;        for (i = 0; I <= m_num_order; i++) {m_py[0] = M_py[0] + m_pnum[i] * M_px[i];        } for (i = 1; I <= m_den_order; i++) {m_py[0] = m_py[0]-m_pden[i] * m_py[i];        } for (i = m_num_order; I >= 1; i--) {m_px[i] = m_px[i-1]; } for (i = m_den_order; I >= 1; i--) {M_py[i] = M_py[i-1];    } Data_out[k] = m_py[0]; }}

In addition, a resp function is provided to calculate the time domain response of the filter. And does not require the data_in to be the same length as the data_out array. the data before default data_in[0] is all 0, andthe number afterdata_in[m-1] is all DATA_IN[M-1]. Therefore, using this function to calculate the step response input data only needs to provide a data point on the line. And this function does not change the internal state of the filter.

/** \brief calculates the time domain response of the IIR filter without affecting the internal state of the filter * \param data_in is the input of the filter, input defaults to 0 hours before the default is 0,data_in[m] and subsequent input defaults to data_in[m-1] * \param D Output of the ata_out filter * \param M length of input data * \param N length of output data * \return */void IIR_I::RESP (double data_in[], int M, double data_out [], int N) {    int i, k, IL;    for (k = 0; k < N; k++)    {        data_out[k] = 0.0;        for (i = 0; I <= m_num_order; i++)        {            if (k-i >= 0)            {                il = ((k-i) < m)? (k-i): (M-1);                DATA_OUT[K] = Data_out[k] + m_pnum[i] * Data_in[il];            }        }        for (i = 1; I <= m_den_order; i++)        {            if (k-i >= 0)            {                data_out[k] = data_out[k]-m_pden[i] * da Ta_out[k-i];}}}    

DirectIItype

/**< IIR Filter Direct type II realization */class iir_ii{public:    Iir_ii ();    void Reset ();    void Setpara (double num[], int num_order, double den[], int den_order);    void resp (double data_in[], int m, double data_out[], int n);    Double filter (double data);    void filter (double data[], int len);    void filter (double data_in[], double data_out[], int len);p rotected:    private:    double *m_pnum;    Double *m_pden;    Double *M_PW;    int m_num_order;    int m_den_order;    int m_n;}; Class Iir_bode{private:    double *m_pnum;    Double *m_pden;    int m_num_order;    int m_den_order;    Complex<double> Poly_val (double p[], int order, double Omega);p ublic:    iir_bode ();    void Setpara (double num[], int num_order, double den[], int den_order);    Complex<double> Bode (double omega);    void Bode (double omega[], int n, complex<double> resp[]);};

The storage units required in the direct type II implementation are much less. In addition, the external interfaces of the two implementations are identical.

Iir_ii::iir_ii () {//ctor m_pnum = NULL;    M_pden = NULL;    M_PW = NULL;    M_num_order =-1;    M_den_order =-1; M_n = 0;};/ * * \brief The internal state of the filter is zeroed, the coefficient of the filter remains * \return */void Iir_ii::reset () {for (int i = 0; i < m_n; i++) {M_pw[i] =    0.0; }}/** \brief * * \param num molecular polynomial coefficients, ascending order, num[0] is constant term * \param m the order of the polynomial * \param den denominator polynomial coefficients, ascending order, den[0] is a constant term * \param     The order of the M-denominator polynomial * \return */void Iir_ii::setpara (double num[], int num_order, double den[], int den_order) {delete[] m_pnum;    Delete[] M_pden;    Delete[] M_PW;    M_num_order = Num_order;    M_den_order = Den_order;    M_n = Max (Num_order, Den_order) + 1;    M_pnum = new Double[m_n];    M_pden = new Double[m_n];    M_PW = new Double[m_n];        for (int i = 0; i < m_n; i++) {m_pnum[i] = 0.0;        M_pden[i] = 0.0;    M_pw[i] = 0.0;    } for (int i = 0; I <= num_order; i++) {m_pnum[i] = Num[i];  } for (int i = 0; I <= den_order; i++) {m_pden[i] = Den[i];  }}/** \brief calculates the time domain response of the IIR filter without affecting the internal state of the filter * \param data_in is the input of the filter, input defaults to 0 hours before the default is 0,data_in[m] and subsequent input defaults to data_in[m-1] * \par Output of the AM data_out filter * \param M length of input data * \param N length of output data * \return */void IIR_II::RESP (double data_in[], int M, double dat    a_out[], int N) {int i, k, IL;        for (k = 0; k < N; k++) {data_out[k] = 0.0;  for (i = 0; I <= m_num_order; i++) {if (k-i >= 0) {il = ((k-i) < M)?                (k-i): (M-1);            DATA_OUT[K] = Data_out[k] + m_pnum[i] * Data_in[il];                }} for (i = 1; I <= m_den_order; i++) {if (k-i >= 0) {            DATA_OUT[K] = data_out[k]-m_pden[i] * data_out[k-i]; }}}}/** \brief filter function with direct type II structure * * \param data input * \return filtered results */double iir_ii::filter (double data) {m    _PW[0] = data; for (int i = 1; I <= m_den_order; i++)//update w[n] Status {m_pw[0] = m_pw[0]-m_pden[i]* M_pw[i];    } data = 0.0;    for (int i = 0; I <= m_num_order; i++) {data = data + m_pnum[i] * M_pw[i];    } for (int i = m_n-1; I >= 1; i--) {m_pw[i] = m_pw[i-1]; } return data; /** \brief filter function, using the direct type II structure * * \param data[] Incoming input data, the return of the result of the filter * \param len data[] array length * \return */void iir_ii::filter (d    Ouble data[], int len) {int i, K;        for (k = 0; k < Len; k++) {m_pw[0] = data[k];        for (i = 1; I <= m_den_order; i++)//update w[n] Status {m_pw[0] = m_pw[0]-m_pden[i] * m_pw[i];        } Data[k] = 0.0;        for (i = 0; I <= m_num_order; i++) {data[k] = Data[k] + m_pnum[i] * M_pw[i];        } for (i = m_n-1; I >= 1; i--) {m_pw[i] = m_pw[i-1]; }}}/** \brief filter function, using direct type II structure * * \param data_in[] Input data * \param data_out[] Save filtered data * \param len array length * \return */voi    D iir_ii::filter (double data_in[], double data_out[], int len) {int i, K; FoR (k = 0; k < Len; k++) {m_pw[0] = data_in[k];        for (i = 1; I <= m_den_order; i++)//update w[n] Status {m_pw[0] = m_pw[0]-m_pden[i] * m_pw[i];        } Data_out[k] = 0.0;        for (i = 0; I <= m_num_order; i++) {data_out[k] = Data_out[k] + m_pnum[i] * M_pw[i];        } for (i = m_n-1; I >= 1; i--) {m_pw[i] = m_pw[i-1]; }    }}

Test

Here are a few test examples that first calculate The step response of a 4-step Chebyshev low-pass filter.

void Resp_test (void) {    double b[5] = {0.001836, 0.007344, 0.011016, 0.007344, 0.001836};    Double A[5] = {1.0, -3.0544, 3.8291, -2.2925, 0.55075};    Double X[2] = {1.0, 1.0};    Double y[100];    Iir_ii filter;    Filter.setpara (b, 4, A, 4);    FILTER.RESP (x, 2, y, +);    for (int i= 0; i<; i++)    {        cout << y[i] << Endl;    }}

The results are as follows:

It is also the filter that calculates the result of the input signal as the Delta function.

void Filter_test (void) {    double b[5] = {0.001836, 0.007344, 0.011016, 0.007344, 0.001836};    Double A[5] = {1.0, -3.0544, 3.8291, -2.2925, 0.55075};    Double x[100], y[100];    int len = +;    Iir_i filter;    Iir_ii filter;    Filter.setpara (b, 4, A, 4);    for (int i = 0; i < len; i++)    {        x[i] = 0.0;        Y[i] = 0.0;    }    X[0] = 1.0;    Filter.filter (x, y, +);    Filter.reset ();    Filter.filter (x, n);    for (int i = 0; i < len; i++)    {        cout << x[i] << "," << y[i]<<  Endl;    }}

The results are as follows:

Implementation of IIR Filter (c + +)