Technical Article

Practical FIR Filter Design: Part 1 - Design with Octave or Matlab

January 24, 2016 by Tim Youngblood

A simple introduction to designing FIR filters in Octave or Matlab

This tutorial will focus on designing a finite impulse response (FIR) filter. As the series progresses, it will discuss the necessary steps to implement the filter on real hardware.

A FIR filter is a digital filter whose impulse response settles to zero in finite time as opposed to an infinite impulse response filter (IIR), which uses feedback and may respond indefinitely to an input signal.The great thing about FIR filters is that they are inherently stable and can easily be designed to have linear phase. I won't get into the details much further on FIR filters and their pro's and con's as this tutorial focuses more on designing filters fast and efficiently with the aid of Octave

Typically, in FIR filter design the length of the filter will need to be specified. You can guess and check until the filter matches your expected bandwidth and cutoff requirements, but this could be a long and tedious process. The equation below is an efficient way to compute a reasonable starting length. After trying the calculated N, one can then tweak N or parameters which make up N to meet filter specifications.

 

$$N \approx \frac{A_{dB} F_{s}}{22 \Delta f}$$

 

where

$$A_{dB}$$ is your stopband attenuation in dB

$$F_{s}$$ is your sampling frequency

and

$$\Delta f$$ is your transition bandwidth 

Lets say we want to filter an audio signal with the following characteristics and desired filter response:

- The samples may contain frequencies from 0-20kHz

- We wish to design a filter that passes only frequencies less than 10kHz

- We want a stopband attenuation of 40 dB at 15kHz

- Our system has a $$F_{s}$$ of 192kHz

With this information we can begin the process of design. Using the equation for N we estimate the filter length to approximately be:

 

$$N \approx \frac{(40) (192000Hz)}{22 (15kHz - 10kHz)}$$

$$N \approx  69.818$$

rounding to the nearest odd number

$$N \approx 69$$

*Designing an FIR filter length to be odd length will give the filter an integral delay of (N-1)/2. 

 

Using the Octave/Matlab code below, we can see how to design a lowpass filter with a bandwidth of 10kHz and a cutoff of 15kHz using Octave's built in fir1 function, which is well documented here

close all;
clear all;
clf;


f1 = 10000;
f2 = 15000;
delta_f = f2-f1;
Fs = 192000;
dB  = 40;
N = dB*Fs/(22*delta_f);

f =  [f1 ]/(Fs/2)
hc = fir1(round(N)-1, f,'low')


figure
plot((-0.5:1/4096:0.5-1/4096)*Fs,20*log10(abs(fftshift(fft(hc,4096)))))
axis([0 20000 -60 20])
title('Filter Frequency Response')
grid on

FIR_filter.zip

The code above gives us the following response:

Figure 1

But if we zoom in we will see that the attenuation at 10kHz is greater than 3dB:

\

Figure 2

The bandwidth of the filter is always specified to the -3dB point, so in the first iteration of design our filter has a smaller bandwidth than specified (somewhere less than 9kHz). We can see from the first figure that the attenuation in the stopband exceeded our specifications, perhaps we can tweak the attenuation and passband frequency to enhance the design. 

 

Tweaking the parameters in the code to

f1 = 11200;
f2 = 15000;
dB  = 30;

gives us a filter which closely matches our speicfications 

Figure 3

The filter design is now complete. Let's simulate how it works by adding the code below to the first bit of code we looked at.

x = sin(2*pi*[1:1000]*5000/Fs) +  sin(2*pi*[1:1000]*2000/Fs) + sin(2*pi*[1:1000]*13000/Fs)  + sin(2*pi*[1:1000]*18000/Fs);

sig = 20*log10(abs(fftshift(fft(x,4096))));
xf = filter(hc,1,x)

figure
subplot(211)
plot(x)
title('Sinusoid with frequency components 2000, 5000, 13000, and 18000 Hz')


subplot(212)
plot(xf)
title('Filtered Signal')
xlabel('time')
ylabel('amplitude')


x= (x/sum(x))/20
sig = 20*log10(abs(fftshift(fft(x,4096))));
xf = filter(hc,1,x)

figure
subplot(211)
plot((-0.5:1/4096:0.5-1/4096)*Fs,sig)
hold on
plot((-0.5:1/4096:0.5-1/4096)*Fs,20*log10(abs(fftshift(fft(hc,4096)))),'color','r')
hold off
axis([0 20000 -60 10])
title('Input to filter - 4 Sinusoids')
grid on
subplot(212)
plot((-0.5:1/4096:0.5-1/4096)*Fs,20*log10(abs(fftshift(fft(xf,4096)))))
axis([0 20000 -60 10])
title('Output from filter')
xlabel('Hz')
ylabel('dB')
grid on

matlab_code.zip

Figure 4

As we can see in Figure 4, we have the time domain signals on the left and the frequency domain on the right. The filter (in red) is overlaid onto the plot to show how the filter leaves the sinusoids in the passband and attenuates the signals in the transition and stopband. 

 

The fir1 function can also be used to produce notch filters, high pass filters, and bandpass filters by replacing these lines:

 

f =  [f1 ]/(Fs/2), may need to be specified with two arguments for bandpass and notch filters as such:

f = [f1 f2]/(Fs/2), where f1 is the left -3dB edge and f2 is the right -3dB edge 


hc = fir1(round(N)-1, f,'low') can be modified as such:

'low' can be replaced with 'stop' (notch), 'high' (highpass), 'bandpass' (bandpass)

 

From the lowpass filter demonstration above it should be easy to form the coefficients (this is the variable hc in the code) for any filter desired. Once you have calculated the coefficients it is important to scale and quantize them so you can implement your filter in a microprocessor. In part 2 we will get into scaling the coefficients and other steps necessary to put your coefficients into an N-bit microprocessor. This is important because without proper scaling you will experience quanitization noise that will affect the frequency response of your filter. Consequently, your design will not match your original specification or the output you simulated in Octave.

Next Article in Series: Practical FIR Filter Design: Part 2 - Implementing Your Filter

3 Comments
  • PeterCoxSmith January 28, 2016

    Thanks, I will give your article some in depth attention. I need to learn more about digital filters.

    Like. Reply
  • DrVigg January 29, 2016

    It is beneficial to work with the zeros of a discrete filter; the zeros of an FIR filter define it to within a gain factor. For a phase linear FIR filter, the impulse response will be symmetric. The zeros come in several categories (the following assumes a sampling rate of 1):

    1. Zeros at z=0. These just introduce delay.
    2. Zeros at +/- 1. These make the DC and Nyquist response, respectively, zero.
    3. Pairs of zeros in the real axis. These must come in pairs if phase linearity is to be maintained; if we locate a real zero at z_0 (and z_0 is not -1, 0, or +1), its reciprocal, 1/z_0, must also be a zero.
    4. Pairs of zeros on the unit circle. To keep the impulse response real, unless you’re including a zero in one of the categories above, you must also include its complex conjugate. So if we have a zero at z_0 = exp (j theta), and theta is not an integer multiple of pi, then we must also include its complex conjugate at exp (-j theta). These zeros will appear as deep notches in the Bode amplitude plot.
    5. Non-real roots off the unit circle. These come in fours. If we wish to include a non-real root at z_0, we must also include its complex conjugate to keep the impulse response real. We must include the reciprocals of z_0 and its complex conjugate, as well, if we wish to maintain phase linearity.

    A corollary to this is that the impulse response of any phase linear FIR filter may be obtained as the convolution of phase linear FIR filters with lengths no greater than 5.

    Like. Reply
    • mirceacluj July 30, 2023
      Dr. Vigg, ​​what do I have to do to understand your comment? In the practical application, I studied the mathematician Chebisev on the Internet with Euler's and De Moivre's formulas, and in addition I recapitulated trigonometry from college. Mathematically, I understand everything, but practically how to calculate a filter. Shall I resume the Matlab study again? and then let's practically make a PCB board with the KiCad software?
      Like. Reply