Standard Band FIR Filter Design: fir1

fir1 uses a least-squares approximation to compute filter coefficients and then smooths the impulse response with a window. For an overview of windows and their properties, see Windows. fir1 resembles the IIR filter design functions in that it is formulated to design filters in standard band configurations: lowpass, bandpass, highpass, and bandstop.

The statements

n = 50;
Wn = 0.4;
b = fir1(n,Wn);

create row vector b containing the coefficients of the order n Hamming-windowed filter. This is a lowpass, linear phase FIR filter with cutoff frequency Wn. Wn is a number between 0 and 1, where 1 corresponds to the Nyquist frequency, half the sampling frequency. (Unlike other methods, here Wn corresponds to the 6 dB point.) For a highpass filter, simply append 'high' to the function's parameter list. For a bandpass or bandstop filter, specify Wn as a two-element vector containing the passband edge frequencies. Append 'stop' for the bandstop configuration.

b = fir1(n,Wn,window) uses the window specified in column vector window for the design. The vector window must be n+1 elements long. If you do not specify a window, fir1 applies a Hamming window.

Kaiser Window Order Estimation.  The kaiserord function estimates the filter order, cutoff frequency, and Kaiser window beta parameter needed to meet a given set of specifications. Given a vector of frequency band edges and a corresponding vector of magnitudes, as well as maximum allowable ripple, kaiserord returns appropriate input parameters for the fir1 function.

Multiband FIR Filter Design: fir2

The fir2 function also designs windowed FIR filters, but with an arbitrarily shaped piecewise linear frequency response. This is in contrast to fir1, which only designs filters in standard lowpass, highpass, bandpass, and bandstop configurations.

The commands

n = 50;
f = [0 .4 .5 1];
m = [1  1  0 0];
b = fir2(n,f,m);

return row vector b containing the n+1 coefficients of the order n FIR filter whose frequency-magnitude characteristics match those given by vectors f and m. f is a vector of frequency points ranging from 0 to 1, where 1 represents the Nyquist frequency. m is a vector containing the specified magnitude response at the points specified in f. (The IIR counterpart of this function is yulewalk, which also designs filters based on arbitrary piecewise linear magnitude responses. See IIR Filter Design for details.)

Basic Configurations

The default mode of operation of firls and firpm is to design type I or type II linear phase filters, depending on whether the order you want is even or odd, respectively. A lowpass example with approximate amplitude 1 from 0 to 0.4 Hz, and approximate amplitude 0 from 0.5 to 1.0 Hz is

n = 20;                    % Filter order
f = [0 0.4 0.5 1];         % Frequency band edges
a = [1  1  0 0];           % Amplitudes
b = firpm(n,f,a);

From 0.4 to 0.5 Hz, firpm performs no error minimization; this is a transition band or “don't care” region. A transition band minimizes the error more in the bands that you do care about, at the expense of a slower transition rate. In this way, these types of filters have an inherent trade-off similar to FIR design by windowing.

To compare least squares to equiripple filter design, use firls to create a similar filter. Type

bb = firls(n,f,a);

and compare their frequency responses using FVTool:

fvtool(b,1,bb,1)

Note that the y-axis shown in the figure below is in Magnitude Squared. You can set this by right-clicking on the axis label and selecting Magnitude Squared from the menu.

The filter designed with firpm exhibits equiripple behavior. Also note that the firls filter has a better response over most of the passband and stopband, but at the band edges (f = 0.4 and f = 0.5), the response is further away from the ideal than the firpm filter. This shows that the firpm filter's maximum error over the passband and stopband is smaller and, in fact, it is the smallest possible for this band edge configuration and filter length.

Think of frequency bands as lines over short frequency intervals. firpm and firls use this scheme to represent any piecewise linear frequency-response function with any transition bands. firls and firpm design lowpass, highpass, bandpass, and bandstop filters; a bandpass example is

f = [0 0.3  0.4  0.7  0.8  1];   % Band edges in pairs
a = [0  0    1    1    0   0];   % Bandpass filter amplitude

Technically, these f and a vectors define five bands:

Example highpass and bandstop filters are

f = [0 0.7  0.8  1];            % Band edges in pairs
a = [0  0    1   1];            % Highpass filter amplitude
f = [0 0.3  0.4  0.5  0.8  1];  % Band edges in pairs
a = [1  1    0    0    1   1];  % Bandstop filter amplitude

An example multiband bandpass filter is

f = [0 0.1 0.15 0.25 0.3 0.4 0.45 0.55 0.6 0.7 0.75 0.85 0.9 1];
a = [1  1   0    0    1   1   0    0    1   1   0    0    1  1];

Another possibility is a filter that has as a transition region the line connecting the passband with the stopband; this can help control “runaway” magnitude response in wide transition regions:

f = [0 0.4 0.42 0.48 0.5  1];
a = [1 1 0.8 0.2 0 0];   % Passband, linear transition, 
                         %  stopband

The Weight Vector

Both firls and firpm allow you to place more or less emphasis on minimizing the error in certain frequency bands relative to others. To do this, specify a weight vector following the frequency and amplitude vectors. An example lowpass equiripple filter with 10 times less ripple in the stopband than the passband is

n = 20;              % Filter order
f = [0 0.4 0.5 1];   % Frequency band edges
a = [1  1   0  0];   % Amplitudes
w = [1 10];          % Weight vector
b = firpm(n,f,a,w);

A legal weight vector is always half the length of the f and a vectors; there must be exactly one weight per band.

Anti-Symmetric Filters / Hilbert Transformers

When called with a trailing 'h' or 'Hilbert' option, firpm and firls design FIR filters with odd symmetry, that is, type III (for even order) or type IV (for odd order) linear phase filters. An ideal Hilbert transformer has this anti-symmetry property and an amplitude of 1 across the entire frequency range. Try the following approximate Hilbert transformers and plot them using FVTool:

b = firpm(21,[0.05 1],[1 1],'h');       % Highpass Hilbert
bb = firpm(20,[0.05 0.95],[1 1],'h');   % Bandpass Hilbert
fvtool(b,1,bb,1)

You can find the delayed Hilbert transform of a signal x by passing it through these filters.

fs = 1000;            % Sampling frequency
t = (0:1/fs:2)';      % Two second time vector
x = sin(2*pi*300*t);  % 300 Hz sine wave example signal
xh = filter(bb,1,x);  % Hilbert transform of x

The analytic signal corresponding to x is the complex signal that has x as its real part and the Hilbert transform of x as its imaginary part. For this FIR method (an alternative to the hilbert function), you must delay x by half the filter order to create the analytic signal:

xd = [zeros(10,1); x(1:length(x)-10)];	    % Delay 10 samples
xa = xd + j*xh;                            % Analytic signal

This method does not work directly for filters of odd order, which require a noninteger delay. In this case, the hilbert function, described in Hilbert Transform, estimates the analytic signal. Alternatively, use the resample function to delay the signal by a noninteger number of samples.

Differentiators

Differentiation of a signal in the time domain is equivalent to multiplication of the signal's Fourier transform by an imaginary ramp function. That is, to differentiate a signal, pass it through a filter that has a response H(ω) = jω. Approximate the ideal differentiator (with a delay) using firpm or firls with a 'd' or 'differentiator' option:

b = firpm(21,[0 1],[0 pi],'d');

For a type III filter, the differentiation band should stop short of the Nyquist frequency, and the amplitude vector must reflect that change to ensure the correct slope:

bb = firpm(20,[0 0.9],[0 0.9*pi],'d');

In the 'd' mode, firpm weights the error by 1/ω in nonzero amplitude bands to minimize the maximum relative error. firls weights the error by (1/ω)² in nonzero amplitude bands in the 'd' mode.

The following plots show the magnitude responses for the differentiators above.

fvtool(b,1,bb,1)
legend('Odd order','Even order','Location','best')

Basic Lowpass and Highpass CLS Filter Design

The most basic of the CLS design functions, fircls1, uses this technique to design lowpass and highpass FIR filters. As an example, consider designing a filter with order 61 impulse response and cutoff frequency of 0.3 (normalized). Further, define the upper and lower bounds that constrain the design process as:

To approach this design problem using fircls1, use the following commands:

n = 61;
wo = 0.3;
dp = 0.02;	
ds = 0.008;
h = fircls1(n,wo,dp,ds);
fvtool(h,1)

Note that the y-axis shown below is in Magnitude Squared. You can set this by right-clicking on the axis label and selecting Magnitude Squared from the menu.

Multiband CLS Filter Design

fircls uses the same technique to design FIR filters with a specified piecewise constant magnitude response. In this case, you can specify a vector of band edges and a corresponding vector of band amplitudes. In addition, you can specify the maximum amount of ripple for each band.

For example, assume the specifications for a filter call for:

Design a CLS filter with impulse response order 129 that meets these specifications:

n = 129;
f = [0 0.3 0.5 0.7 0.9 1];
a = [0 0.5 0 1 0];
up = [0.005 0.51 0.03 1.02 0.05];
lo = [-0.005 0.49 -0.03 0.98 -0.05];
h = fircls(n,f,a,up,lo);
fvtool(h,1)

Note that the y-axis shown below is in Magnitude Squared. You can set this by right-clicking on the axis label and selecting Magnitude Squared from the menu.

Weighted CLS Filter Design

Weighted CLS filter design lets you design lowpass or highpass FIR filters with relative weighting of the error minimization in each band. The fircls1 function enables you to specify the passband and stopband edges for the least squares weighting function, as well as a constant k that specifies the ratio of the stopband to passband weighting.

For example, consider specifications that call for an FIR filter with impulse response order of 55 and cutoff frequency of 0.3 (normalized). Also assume maximum allowable passband ripple of 0.02 and maximum allowable stopband ripple of 0.004. In addition, add weighting requirements:

To approach this using fircls1, type

n = 55;
wo = 0.3;
dp = 0.02;
ds = 0.004;
wp = 0.28;
ws = 0.32;
k = 10;
h = fircls1(n,wo,dp,ds,wp,ws,k);
fvtool(h,1)

Note that the y-axis shown below is in Magnitude Squared. You can set this by right-clicking on the axis label and selecting Magnitude Squared from the menu.

Multiband Filter Design

Consider a multiband filter with the following special frequency-domain characteristics.

A linear-phase multiband filter may be designed using the predefined frequency-response function multiband, as follows:

b = cfirpm(38, [-1 -0.5 -0.4 0.3 0.4 0.8], ...
               {'multiband', [5 1 2 2 2 1]}, [1 10 5]);

For the specific case of a multiband filter, we can use a shorthand filter design notation similar to the syntax for firpm:

b = cfirpm(38,[-1 -0.5 -0.4 0.3 0.4 0.8], ...
              [5 1 2 2 2 1], [1 10 5]);

As with firpm, a vector of band edges is passed to cfirpm. This vector defines the frequency bands over which optimization is performed; note that there are two transition bands, from –0.5 to –0.4 and from 0.3 to 0.4.

In either case, the frequency response is obtained and plotted using linear scale in FVTool:

fvtool(b,1)

Note that the range of data shown below is (-pi,pi).

The filter response for this multiband filter is complex, which is expected because of the asymmetry in the frequency domain. The impulse response, which you can select from the FVTool toolbar, is shown below.

Filter Design with Reduced Delay

Consider the design of a 62-tap lowpass filter with a half-Nyquist cutoff. If we specify a negative offset value to the lowpass filter design function, the group delay offset for the design is significantly less than that obtained for a standard linear-phase design. This filter design may be computed as follows:

b = cfirpm(61,[0 0.5 0.55 1],{'lowpass',-16});

The resulting magnitude response is

fvtool(b,1)

The y-axis is in Magnitude Squared, which you can set by right-clicking on the axis label and selecting Magnitude Squared from the menu.

The group delay of the filter reveals that the offset has been reduced from N/2 to N/2-16 (i.e., from 30.5 to 14.5). Now, however, the group delay is no longer flat in the passband region. To create this plot, click the Group Delay Response button on the toolbar.

If we compare this nonlinear-phase filter to a linear-phase filter that has exactly 14.5 samples of group delay, the resulting filter is of order 2*14.5, or 29. Using b = cfirpm(29,[0 0.5 0.55 1],'lowpass'), the passband and stopband ripple is much greater for the order 29 filter. These comparisons can assist you in deciding which filter is more appropriate for a specific application.