IIR filter Software implementation (matlab+c++)

Use C + + to write an IIR filter

We first design an IIR filter in MATLAB and generate a header file that reflects the frequency response characteristics of the IIR filter.

Theoretical support

IIR filter is called recursive filter, it is a kind of filter with feedback. When the order is large, it can be used in series with multiple second-section filter to improve the stability of the system.

A second-order nodal coefficient law:

You can write the difference equation for the K-Second second section

The cascade structure of N second-order sections is as follows:

According to the second-order section diagram, the output of the previous level as the input of the latter level, the function of IIR digital filtering can be realized by software.

Using MATLAB to generate header files

First open the filter Design & Analysis Tool in MATLAB

Here we first design a low-pass filter

FS represents the sampling frequency and the sampling frequency must be greater than twice times the maximum frequency of the original signal,

Otherwise, spectral aliasing is generated.

Fpass is the pass-band frequency and the fstop is the cutoff frequency of the resistance band

These parameters are set so that you can click on the design Filter

Generated is a second-order section filter combination, a total of 31 orders, that is, the combination of multiple second filter

The C Header File is then generated in the target option.

Numerator is the name of the array of molecular coefficients, and numerator length is an array of molecular coefficients,

Denominator is the denominator.

To parse the build header file

The following is a low-pass filter with a fpass of 10k,fstop of 12K 栵

Before using the header file, you need to convert the data type of MATLAB to the data type supported by C + + as appropriate, where we use the double type

Before analyzing the header file, look at the first section of the filter parameters provided by Matlab

Examples of data with the first second-order section:

• Numerator:1 2 1

• Denominator:1-0.55930961405944157 0.92579835996619642

• gain:0.34162218647668868

Numerator are the coefficients of the molecules, respectively, B0,B1,B2

Denominator are the coefficients of the denominator, respectively, of A0,A1,A2.

The gain is a gain for each section, which stabilizes the signal size in order to stabilize the sections

Then against the header file, the following is the main part of the header file interception:

#defineMwspt_nsec 41intNl[mwspt_nsec] = {1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,3,1,2,1 };Doublenum[mwspt_nsec][3] = {  {     0.3416221864767,0,0   },  {                   1,2,1   },  {     0.3180955154747,0,0   },  {                   1,2,1   },......

Mwspt_nsec is the filter order, the specific number of sections in the comments at the beginning of the header file

NL[MWSPT_NSEC] This array defines the number of useful data for each row of the Num[mwspt_nsec][3] array (can not be used)

In Num[mwspt_nsec][3] Array (molecular parameters) odd line The first item is the gain term, the even behavior of 3 coefficients, respectively, is b0,b1,b2.

Thus can find the law, define K for the current order, p is the first pointer to the array, then the gain term for each section is (P+6*K), the first coefficient is (p+3+6*k),

The second coefficient is (p+3+6*k+1), and the third factor is (p+3+6*k+2).

C + + Programming implementation

In the process of software design, the delay variable of each second-order section is only taken as an intermediate variable and is directly assigned in the process. This is because the delay variable for the next input data n+1 is the same as the previous input data and is designed in this way to save space on the register.

In order to improve the processing speed, it is necessary to use pointers for parameter passing in the program, paying special attention to the first address of the two-dimensional array is &a[0][0]->double* a.

Filter function

Double IIR (double *a, double *b,double* W, double xin, int n_iir)

{

int k;

Double temp = xin;

for (k = 0; k<n_iir; k++)

{

* (w+k*3) = temp-* (A + 3+6 * k + 1) * (* (w + k * 3+1))-* (A + 3 + 6 * k + 2) * (* (w + k * 3+2));

Here, temp is the input to this second-order section and the output of the last second section.

temp = * (b + 3 + 6 * k) * (* (w + k * 3)) + * (b + 3 + 6 * k + 1) * (* (w + k * 3+1)) + * (b + 3+6 * k + 2) * (* (w + k * 3+2));

Here, temp is the output of this second-order section and the input to the next second section.

* (W + k * 3 + 2) = * (W + k * 3 + 1);

* (W + k * 3 + 1) = * (W + k * 3);

temp = temp* (* (b + 6 * k));//magnification, stable signal

}

return temp;

}

Actual test

Test Fpass for a low-pass filter of 10k,fstop 12K

In the program input three frequency 2k,11k,20k signal, theoretically 2k completely through, 11k partial attenuation, 20K completely filtered.

is the original signal, the signal after filtering.

The actual test findings meet the design requirements, and the transition band signal is basically completely attenuated.

Testing in C + +programvoidMain () {Const intN = -; inti,j; DoubleXn[n]; Doublew[ -][3]; DoubleYn[n];  for(i =0; I < -; i++)//Initialize    {         for(j =0; J <3; J + +) W[i][j]=0; }     for(i =0; i < N; i++) {Xn[i]= Sin (2*3.1416* -/ -* i) + sin (2*3.1416*2/ -* i) + sin (2*3.1416* One/ -*i); Yn[i]=iir (&den[0][0], &num[0][0], &w[0][0],xn[i], -); } ofstream savefile_a ("Xn.txt");  for(i =0; i<n; i++) Savefile_a<<" "<< Xn[i] <<Endl;    Savefile_a.close (); Ofstream Savefile_b ("Yn.txt");  for(i =0; i<n; i++) Savefile_b<<" "<< Yn[i] <<Endl; Savefile_a.close ();} MATLAB program for analysis Xn1=fopen ('Xn.txt','R'); [Xn,count]=FSCANF (XN1,'%f'); fclose (XN1); N= Length (xn);Find sampling points xn_f= FFT (xn);%Fourier transform the signal to Xn_f=abs (Xn_f (1: n/2) ); F=50000/n* (0: n/2-1); subplot (211); Stem (f,abs (xn_f)); Xlabel ('Frequency/(s)'); Ylabel ('amplitude'); Title ('Original Signal Spectrum'); Grid;yn1=fopen ('Yn.txt','R'); [Yn,count]=FSCANF (YN1,'%f'); fclose (YN1); Yn_f= FFT (yn);%Fourier transform the signal to Yn_f=abs (Yn_f (1: n/2)); subplot (212); Stem (f,abs (yn_f)); Xlabel ('Frequency/(s)'); Ylabel ('amplitude'); Title ('signal spectrum after filtering'); grid;

IIR filter Software implementation (matlab+c++)