Design and Implementation of FIR Filter Based on MATLAB

I. Summary

Previous ArticleArticleThis article introduces how to design a filter using the fdatool toolbox. For details, see"Filter Design and file generation based on fdatool toolbox in MATLABHere are several examples to illustrate the process of designing an FIR Filter Using Matlab.

Ii. Experimental Platform

Http://dt.ap-southeast-1.maxcompute.aliyun-inc.com

Iii. Experiment Principles

Taking a low-pass filter as an example, its common design indicators include:

1. Band edge frequency FP (digital frequency: Ω P)
2. FST (digital frequency: Ω st)
3. Maximum ripple attenuation delta P =-20log10 (1-α P) in the passband, measured in dB
4. Minimum attenuation α S =-20log10 (α S) of the impedance band, in dB
5. Resistance band Fluctuation α S
6. Band Peak Fluctuation α P

Among them, 1, 2, 3, and 4 are the most commonly used. Article 5 and Article 6ProgramUsed to estimate the order of the filter.

Digital frequency = analog frequency/sampling frequency

Iv. instance analysis Example 1 uses the CASET window to design a FIR low-pass filter, with the boundary frequency Ω P = 0.3pi, the boundary frequency Ω S = 0.5pi, and the band attenuation delta s is not less than 50db. Method 1: manually calculate the Order N and β values of the filter, and then design the filter through the program.

Step 1: Calculate the order B and β values of the filter based on the transition band width and impedance band attenuation.

Step 2: design the filter through programming.

The procedure is as follows:

B = FIR1 (29, 0.4, Kaiser (30, 4.55 ));

[H1, W1] = freqz (B, 1 );

Plot (W1/PI, 20 * log10 (ABS (H1 )));

Axis ([0, 1,-80, 10]);

Grid;

Xlabel ('normalization Frequency/p ');

Ylabel ('amplitude/db ');

 

The waveforms are as follows:

 

Method 2:

Use the [N, Wn, beta, FTYPE] = kaiserord (F, A, Dev) function to estimate the order of the filter, and obtain the CASET window filter.

The functions kaiserord (F, A, Dev) or kaiserord (F, A, Dev, FS) are as follows ):

F indicates the corresponding frequency, and FS indicates the sampling frequency. When F is expressed as a digital frequency, FS does not need to be written.

A = [1 0] indicates the amplitude vector of each frequency band specified by F, which is generally expressed only by 0 and 1. The length relationship between A and F is (the length of 2 *) -2 = (F length)

Devs = [0.05 10 ^ (-2.5)] is used to specify the maximum output error or deviation between the frequency response of each band output filter and its expected amplitude, And the length is equal to the phase. The formula is as follows:

Band attenuation error = α S, band attenuation error = α P, 3 or 4 of the filter indicators can be obtained.

The default value of FS is 2Hz.

 

The procedure is as follows:

Fcuts = [0.3 0.5]; % normalized frequency Omega/PI, which indicates the ending frequency and the starting frequency of the blocked band.

Mags = [1 0];

Devs = [0.05 10 ^ (-2.5)];

[N, Wn, beta, FTYPE] = kaiserord (fcuts, Mags, Devs); % calculate the values of Kill window N and beta

HH = FIR1 (n, Wn, FTYPE, Kaiser (n + 1, beta), 'noscale ');

Freqz (HH );

 

The waveforms are as follows:

 

In practice, the MATLAB Signal Processing Toolbox function remezord is generally called to calculate the Order N and the weight function W (ω) of the same ripple filter, and the Remez function can be called to design the same ripple filter, directly obtain the filter coefficient. In the remezord function, the array fedge is the boundary frequency of the pass band and the drag band. The array mval is the amplitude of the two boundary, the array Dev is the fluctuation of the pass band and the drag band, and the FS is the sampling frequency unit is Hz.

 

Example 2 alternate use of lumizAlgorithmDesign an equi-ripple filter and design a linear phase low-pass FIR digital filter with the following indicators: Band-through boundary frequency fc = 800Hz, band-through boundary Fr = 1000Hz, the minimum attenuation at = 40db and the sampling frequency fs = 4000Hz.

 

Solution: remezord and Remez functions can be used in MATLAB.

 

The procedure is as follows: 

 

Fedge = [800 1000];

 

Mval = [1 0];

 

Dev = [0.0559 0.01];

 

FS = 4000;

 

[N, fpts, Mag, Wt] = remezord (fedge, mval, Dev, FS );

 

B = Remez (n, fpts, Mag, wt );

 

[H, w] = freqz (B, 1,256 );

 

Plot (w * 2000/PI, 20 * log10 (ABS (H )));

 

Grid;

 

Xlabel ('frequency/Hz ');

 

Ylabel ('amplitude/db ');

 

The waveforms are as follows:

Example 3 using Matlab programming to design a digital band-pass filter, the index requirements are as follows: band edge frequency: Ω p1 = 0.45pi, Ω P2 = 0.65pi, band peak fluctuations: delta 1 <= 1 [dB]. Band edge frequency: Ω S1 = 0.3pi, Ω S2 = 0.8pi, minimum band attenuation: Delta 2> = 40 [dB]. Method 1: Window Function

 

The procedure is as follows:

[N, Wn, BTA, FTYPE] = kaiserord ([0.3 0.45 0.65 0.8], [0 1 0], [0.01 0.1087 0.01]); % use the kaiserord function to estimate the Order N and beta parameters of the filter.

H1 = FIR1 (n, Wn, FTYPE, Kaiser (n + 1, BTA), 'noscale ');

[HH1, W1] = freqz (H1, 1,256 );

Figure (1)

Subplot (2, 1, 1)

Plot (W1/PI, 20 * log10 (ABS (HH1 )))

Grid

Xlabel ('normalization frequency W'); ylabel ('amplitude/db ');

Subplot (2, 1, 2)

Plot (W1/PI, angle (HH1 ))

Grid

Xlabel ('normalization frequency W'); ylabel ('Phase/rad ');

 

The waveforms are as follows:

 

Filter coefficient:

H1 =

 

Columns 1 through 8

 

0.0041 0.0055-0.0091-0.0018-0.0056-0.0000 0.0391-0.0152

 

Columns 9 through 16

 

-0.0381 0.0077-0.0293 0.0940 0.0907-0.2630-0.0517 0.3500

 

Columns 17 through 24

 

-0.0517-0.2630 0.0907 0.0940-0.0293 0.0077-0.0381-0.0152

 

Columns 25 through 31

 

0.0391-0.0000-0.0056-0.0018-0.0091 0.0055 0.0041

 

If freqz (H1, 1,256) is used directly, the amplitude-frequency and phase-frequency curves are obtained:

  Method 2:Same ripple Method Design

 

The procedure is as follows:

[N, fpts, Mag, Wt] = remezord ([0.3 0.45 0.65 0.8], [0 1 0], [0.01 0.1087 0.01]); % use the remezord function to estimate the Order N, normalized band edge vector fpts, within-band amplitude response vector mag, and Weighted Vector W required by the Remez function, the filter designed by the Remez function meets the performance requirements specified by F, A, and Dev.
H2 = Remez (n, fpts, Mag, wt); % design the same ripple Filter
[HH2, W2] = freqz (H2, 1,256 );
Figure (2)
Subplot (2, 1, 1)
Plot (W2/PI, 20 * log10 (ABS (HH2 )))
Grid
Xlabel ('normalization frequency W'); ylabel ('amplitude/db ');
Subplot (2, 1, 2)
Plot (W2/PI, angle (HH2 ))
Grid
Xlabel ('normalization frequency W'); ylabel ('Phase/rad ');
H2

 

The waveforms are as follows:

The filter coefficients are as follows:

H2 =

 

Columns 1 through 9

 

-0.0013 0.0092-0.0255-0.0642 0.1177-0.0922-0.2466-0.0466 0.3116

 

Columns 10 through 17

 

-0.0466-0.2466 0.0922 0.1177-0.0642-0.0255 0.0092-0.0013

 

If freqz (H2, 1,256) is used directly, the amplitude-frequency and phase-frequency curves are obtained:

 

Method 3: Use fdatool

 In this method, the order and bate values of the filter need to be calculated in advance, and then the corresponding parameters are set to generate the filter.

Shows the setting interface:

 

After the area of the above circle is set, the filter is generated. The analysis menu can be used to observe the various characteristic curves and filter coefficients of the generated filter. The filter coefficients here are the same as those in method 1.

 

The waveforms are as follows:

 

V. Result Analysis 5.1 filter design summary

FIR filters generally adopt the window function method and the same ripple design method. The window function method also contains two branches. One is to use the formula to manually calculateNValues correspond to other window function parameter values, and are then substituted into the window function andFIR1One is to use the function * rord to estimateNAnd use corresponding parametersFIR1. Note that* RordWill underestimate or overestimate the orderNMay cause the filter to fail to reach the specified performance. In this case, the order should be slightly increased or decreased. If the cutoff frequency is0OrNyquestFrequency nearby or setDevIf the value is large, no correct result is returned.

Filter Implementation form and features: Because the ordinary filter in the use of the window function is its through the belt ripple and the drag band ripple is different (generally the first drag band ripple is the largest) Therefore, when the first impedance band is used to attenuation the side lobe, their attenuation is greatly higher than its frequency. Based on the approximate relationship between impedance band attenuation and number of items n = P (delta 2) * fs/TW, the larger the impedance band attenuation, the more items required.

5.2 differences between the window function method and the same ripple Design

The design of the window function is based on the Least Square integral method, that is, the error of the filter is:

That is, the minimum method is required to design the filter, which is more loyal to the ideal filter (that is, the filter coefficient is closer to the ideal filter ).

The proof is as follows:

Therefore, the smaller the amplitude spectrum difference, the closer the actual filter is to the ideal filter.

The same ripple filter is implemented by minimizing the maximum weighted error. Its error is:

The minimum error is required to implement the filter. The obtained filter coefficient is far different from the window function design.

In Example 3H1AndH2For comparison.

% Sigsum is used to sum all elements of the array.

Function Y = sigsum (N1, N2, n, x );

Y = 0;

For I = N1 + 1-min (n): N2 + 1-min (N)

Y = Y + X (I );

End

 

N = 0.001: 30.001;

H = 2 * Cos (0.55 * pI * (n-15). * sin (0.175 * pI * (n-15)./(pI * (n-15 ));

Delta1 = h-h1;

N = 0.001: 16.001;

H = 2 * Cos (0.55 * pI * (n-15). * sin (0.175 * pI * (n-15)./(pI * (n-15 ));

Delta2 = h-h2;

Y1 = sigsum (0, 30, [0: 30], (ABS (delta1). ^ 2)/31;

Y2 = sigsum (0, 16, [0: 16], (ABS (delta2). ^ 2)/17;

 

The result is as follows:

Y1 =

 

1.9099e-004

 

Y2 =

 

0.0278

 

In this way, the filter coefficient implemented by the window function is closer to that of the ideal filter (Y1For the difference between the window function and the ideal filter,Y2For the difference between the same ripple filter and the ideal filter );

 

 

Comparing the amplitude spectrum of the two, we can see that the anti-band edge of the same corrugated filter is smoother than that of the window function (the ideal filter is vertical descent ).

From the design point of view, because the window function design method is to transform the ideal filter through the existing window function, you can use a hand-calculated method to conveniently design the filter.

The implementation of the ripple filter is achieved through a large number of iterative operations. Such a method can only be designed through software.

The problem of the number of items because the Equi-corrugated filter can have an average distribution error, so for the same impedance band attenuation, the filtering coefficient is less than the window function.

5.3 notes

1. Description of different shapes of Phase-Frequency Characteristic Curves

The first figure above is drawn in angle units. The following figure uses Rad Unit. From the graph, we can see that 0.3 To 0.8 The two diagrams of the digital frequency are strictly linear. Why is there a hop in the following figure? Rad Only -Pi -- Pi When the phase -Pi You can only jump Pi But cannot be greater Pi And the angle can be increased continuously.

 

2. CallFirlOrRemeUseScale(Default) Normalize the filter, that is, the response amplitude at the center frequency of the filter is0db. UseNoscaleThe filter is not normalized.