Finite Impulse Response Filters Part2: Sinc Functions and Windows

Introduction

In the last article, I introduced the subject of finite impulse response filters (FIRs) and had a look a convolution. An HDL implementation was then discussed showing just how easy it is to construct logic to implement a solution for convolving two signals—for our purposes a finite impulse response with a continuous stream of digital signal samples. However, I left things hanging as the implementation needs configuring its ‘taps’ (or coefficients) to program the impulse response into the logic, and we don’t know how to do that yet.

In this article I want to look more closely at the sinc function and explain why we’re messing around with it, and then window functions to turn the infinite nature of the sinc into a finite set of values we can configure into the logic implementation. Using a digital audio example to set some requirements on a practical design, we will then use winfilter to generate tap values, looking at the effects of the different window functions, quantisation bit width and the number of taps has on things like the filter’s roll-off rate, stop-band attenuation and ripple in the pass-band. We will be able to plot the frequency response, the phase, the impulse response, and even the window function to visualise the ?functions and the filter performance. Although we will explore a low pass filter in detail as a practical example, we will also look briefly at highpass, bandpass and bandstop filters as well, and how we can construct these. I mentioned in the first article that, actually, an arbitrary frequency response can be implemented using FIRs, but space restricts what we can discuss here, but I will give some references to allow adventurers to continue to explore on their own.

Once we have some taps values for our logic implementation, we can then test how the filter actually performs against our predictions for the digital audio example, and some simulation results will be analysed and cross-referenced. A simulation environment is provided with the Verilog implementation which should be adaptable for a preferred logic simulator.

The sinc Function and Impulse Response

I mentioned in the first article that the impulse response for a low pass filter is a sinc function of general form sine(x)/x. Firstly, this isn’t strictly true as it’s sine(πx)/πx, and secondly why is it a sinc function? In the introduction to the first article there was a diagram of an idealised low pass filter. This diagram is actually only half the story. It only plots the frequency axis from 0 to the sample frequency divided by 2 (fs/2). Without going into too much detail, when sampling data at a given sample frequency, only signals up to half that rate can be represented—this is the Nyquist limit. Any higher frequencies get ‘aliased’, and fold-back into the lower frequencies. What we see when we plot the entire frequency range is something like that shown in the diagram below:

If we simply extend the diagram from the first article to fs, we get a plot like that shown on the left of the diagram above, which is a reflection of the first half. If we kept plotting higher and higher frequencies this would actually repeat periodically, both in the positive and negative directions. The diagram on the right, then, shows what you get if plot from - fs/2 to + fs/2—i.e., a rectangle. Now, the fourier transform of a rectangle is a sinc function. I’m not going to prove this (others who know more mathematics than I do have done so already), but this is why we want a sinc function for our impulse response. The phase for a sinc function centred around 0, is 0 for all frequencies. Even if the sinc function is offset from 0 the phase remains linear, with a slope a function by the offset. This linear phase is one of the defining characteristics of an FIR, as mentioned in the first article. We will look at a phase plot when we discuss the audio example.

The diagrams above don’t specify what the sample or cut-off frequencies are, but as we increase and decrease the passband as a proportion of the sample frequency, we can get a ratio fc/fs. This can now be mapped into our sinc function to get the impulse response for the cut off frequency desired as a fraction of sample frequency:

So, fc/fs can be anything from 0 to 0.5 to determine the low pass filter pass- and stopbands. For example, if our sample rate was 100KHz, and the desired cut off frequency is, say, 12.5KHz, then fc/fs = 0.125. Plug this into the equation and we can calculate the impulse response for all values of i, to plus and minus infinity.

The value when i is zero isn’t just 1—previously, in the first article, we used nice round numbers. The sine function becomes linear the nearer we approach zero radians—i.e., sin(θ) tends to θ. Therefore, removing the sine, and dividing by iπ, leaves the value as shown in the equation for i at zero. So now we’d better deal with infinity.

Windowing and Window Functions

In the first article, I mentioned that truncating the infinite impulse response gave a ‘window’ on that response. I also mentioned that this simple truncation adds undesirable artifacts and that things can be improved with window functions that are a more rounded roll-off. The simple truncation, in fact, is known as a ‘uniform’ window. Over the years, researchers have experimented with various functions to improve the response from the basic uniform window, from a triangle (aka Bartlett window), a cosine window, a raised cosine window (aka Hamming window), through to some quite elaborate functions such as a Kaiser window. We can only look at a few, but the winfilter program supports quite a few, and the functions’ definitions are all in the source code (see windows.c). Starting from winfilter’s default settings, let’s look at a few example windows:

As you can see, the shapes of these windows more or less have a similar form, starting at 1 for the centre point and reducing either side towards 0, and it gets harder to spot the difference between them, though they do make a difference as we shall see. I haven’t plotted the uniform window, as I assume you can picture multiplying by 1 over the number of taps and 0 everywhere else.

If we now take the sinc impulse response that we calculated from the last section and multiply it by the selected window function, we get the finite windowed-sinc filter tap coefficients (sometimes know as the filter’s kernel) with which we need to program our HDL implementation. Plotting the frequency responses (in dBs) for these windowed-sinc impulse responses (including for the uniform window), we get the following results:

I’ve chosen the order of the types of windowed-sinc impulse responses with improved performance from first to last (left to right, then down). The uniform window has ripple in the passband and tails off slowly to an attenuation of -46dB ?or so by the fs/2 point. The Bartlett window is not much better, reaching around -55dB at fs/2, but with a slightly better roll-off. The cosine window improves things again, but with the Hamming window we now have a very decent response, reaching -60dB (for the given Alpha value) with a superior roll-off rate and maintaining that attenuation in the stopband. The Kaiser window, like the Hamming window, has an Alpha (α) value that can be manipulated to trade off passband ripple, roll-off rate, and stopband attenuation. With the values I’ve configured for the example, you can see that superior stopband attenuation can be achieved, with good roll-off and minimal passband ripple.

Phase Characteristics

The diagram below shows the phase plot for the example Kaiser windowed-sinc impulse response over the passband (0 to 24KHz).

The phase axis is in degrees and runs from -180° to + 180°, where the plot then wraps. If you imagine, though, that the axis goes to infinity, then you would see a straight line extending in a linear fashion.

This confirms what was suggested in the first article that the phase is linear. This remains true regardless of the number of taps, the window function used or the quantisation.

Window Mathematics

So far, I’ve taken the window functions as read without describing how they are calculated. The uniform window should be intuitive, as is, I hope, is the triangular window. For the rest, I want to look at the examples for those discussed above, but the comments above each window function in the winfilter source code details the mathematics for all the window functions not explored here. This section can be skipped for those not interested in the details, but just want the tap coefficient values, but I think it is interesting to have a look for informative purposes.

The cosine window isn’t as straight forward as plotting for ±?π, but gets raised to a power, defined by the Alpha setting (1.0 in the case of the example above). I.e.

The Hamming window is a raised cosine where the Alpha (α) value defines how much above 0 the cosine is raised:

A Kaiser window is now going to pile on the mathematics to some degree, using a ‘modified order 0 Bessel function’. The window is calculated as:

Where I0 is a modified order 0 Bessel function, I0,? is calculated as:

The Bessel function runs with k from 0 to infinity which, obviously, can’t be computed, but the terms converging towards 0 as k increases, and a value of 69 terms was heuristically chosen to be used in winfilter.

So you can see, things can get quite elaborate. The Chebyshev window, for example, is actually calculated in the frequency domain, and then converted into the time domain with an inverse discrete fourier transform. I’ve included some examples in this section so that there is an idea of what’s involved. The good thing, of course, is that these calculations do not have to be done as part of the implementation, but can be calculated offline, with the tap coefficients as the result to be configured into the FIR logic implementation. The source code documents all the windows, with C implementation for those that wish to study this more deeply.

Quantisation

All the plots we’ve seen so far all use double precision numbers for maximum accuracy but quantising the impulse response taps to a fixed number of bits (either signed integers or signed fixed-point values—it makes no odds) has an effect on the resultant filter performance. In particular, the stopband attenuation performance suffers. Taking the Kaiser windowed-sinc filter example from before, which reached better than -80dB in the stopband, we can change the quantisation setting and replot the frequency response.

With 16-bit quantisation, we still reached around -80dB attenuation, if with a little messier stopband ‘lobes’, but with 12-bit precision, this becomes around -60dB and with 8 bits around -45dB but rises above this further up the stopband frequencies. Thus the choice of window function and an Alpha value is not enough to reach a particular target performance, and the number of quantisation bits has to be chosen to sufficiently meet the requirements as well.

In theory, the achievable signal to noise ratio for a given number of bits (Q) is around 6Q + 1.76 dB. So a Q of 12 should yield -74dB. In the Q = 12 plot above, it actually does no better than -60dB across the stopband. What I see is that the initial transition will reach around the theoretical limit, but the uneven lobes will often bring this back up above this, so plotting the actual response is important to verify performance target have been met.

High Pass, Band Pass and Band Stop Filters

Until now only low pass filters have been discussed, with the sinc function being the fourier transform of the rectangular shape of the LPF, but what about highpass and bandpass filters?

For highpass filters we could think of this in terms of lowpass filters in one of two ways. Firstly, if we invert an LPF in the frequency domain, i.e. flip it in the y axis, it will suddenly become a high pass filter. Also, if we reversed the frequency response about the fs/2 midpoint, flipping it about the x axis, so that the stopband starts first ending in the pass band, that would also be a highpass filter. What does this mean for in the impulse responses?

For the inversion, in the frequency domain, if we subtracted the low pass filter response from an all-pass filter response (i.e.,? 1) this would achieve an inversion. We saw before that an all-pass filter in the frequency domain is a delta function (an impulse of magnitude 1) in the time domain. So subtracting our low pass impulse response from the delta function will invert the response. This means multiplying all the values by -1, except the centre point, where the delta occurs which is subtracted from 1. The diagram below shows the impulse response for an inverted LPF and its frequency response based on the Kaiser windowed examples from above.

To reverse an LPF, recalling the earlier idealised response diagram and the rectangular shape centred on 0 frequency, if we could shift this so that it was centred on fs/2 instead, that would reverse the response, with (from the earlier diagram) the response from - fs/2 to 0—a reflection of that from 0 to fs/2—now sitting from 0 to fs/2. We can do this if we convolve the LPF frequency response with a delta function at fs/2. A delta function in one domain is a 1 in the other if it is at point 0. If it is shifted, the signal in the other domain becomes rectangular between -1 and 1. A delta function at fs/2 produces a signal in the time domain at a frequency of fs/2, and we know that convolving in one domain means multiplying in the other domain. So if we multiply every other coefficient of an LPF impulse response by -1, we will reverse the frequency domain. This is shown below:

Now that we can make highpass filters, we can combine these to make bandpass and bandstop filters. For a highpass filter, imagine having an LPF and an HPF such that, in the frequency domain, they overlap. If you multiply these responses together then, where they overlap, there will be a pass band and elsewhere a stop band. Multiplying in the frequency domain is the same as convolving in the time domain, so we can take the impulses from the low- and highpass filters and convolve them together to generate the impulse response for the band pass filter. If we invert this band pass kernel, in the manner discussed above, this will become a band stop filter. So, by using convolution, multiplication, and lowpass filters, we can generate these other types of filters as well. This is why the bulk of these articles have concentrated on low pass filters, as the others are just a function of these.

Designing the Filter

To put what we’ve looked at so far together, we will work though an example filter design, using winfilter, from the world of digital audio. Digital audio sample rates perhaps more familiar to most people from consumer electronics and streaming services might be 44.1KHz from Compact Discs (CDs), or 48KHz from DVDs and (in an obsolete reference to my own past) Digital Audio Tape (DAT). Blu-ray supports many audio formats, with varying number of channels at rates from 48KHz to 192KHz. These audio formats have sample widths ranging from 16- to 24-bits. The sample rates for audio will not tax the HDL implementation, which can run at over 750KHz with a 100MHz clock, and so we’ll target an audio application.

As a side note, an FIR that runs at the higher rate of, say, 192KHz can still process signals at a lower (integer) rate such as 48KHz by ‘oversampling’ the signal. That is, inputting an individual sample multiple times (in this case four, for 4 times oversampling) and the filter will work just fine. The advantage of this is that the digital noise of the signal is spread over the larger frequency range. When the signal is converted back to the analogue domain and low pass filtered, the noise at the higher frequencies is attenuated, leaving less noise in the passband than would otherwise have been? without oversampling. The analogue anti-aliasing filter specification can also be relaxed as it has a much longer attenuation region with which to reach its attenuation target to filter aliased frequencies. A simplified system is shown in the diagram below:

In the diagram I’ve ignored any clock synchronisation between the sampling rate and the internal clock, but this would need considering—either the system clock rate can be a synchronous integer multiple of the sample rate, or CDC logic added. CDC can add some time domain jitter to the samples which would need analysing ?to ensure acceptable (beyond the scope of this article).

With oversampling, the specification of the analogue filter can also be relaxed, as mentioned above, as it can roll off more slowly from the cut off frequency, as the space between this and its reflected frequency response between fs/2 to fs is wider than without oversampling. But the analogue LPF must be present to remove the reflected aliased signals. A typical analogue LPF used for this might be a Bessel filter which has a linear phase delay (aka maximally flat group delay) preserving the phase linearity of the FIR. An example 4-pole Bessel filter is shown below, courtesy of the Texas Instruments filter design tool.

Before setting out on the specific audio example, a quick summary of the winfilter program usage is in order.

Note on Using winfilter

Firstly, a note on the winfilter source code. The winfilter source code (originally xfilter for Linux and X-windows) has its origins from more than twenty years ago. In fact, some parts of the code, particularly the DFT/FFT parts, originate from my undergraduate dissertation work done in 1988/’89. Therefore, it’s written in C, rather than C++, and its style is not necessarily how I’d construct code nowadays. What I was gratified to see, when revisiting the code, was that there are plenty of comments, including the mathematics for things such as the window functions alongside the code that implements them. So, if a bit rough around the edges in terms of coding best practices, it is still a useful reference for all the details we can’t cover in this article (we didn’t explore every window function or the DFT and FFTs functionality—this latter we’ll save for another article).

As the name implies, winfilter is a Windows based graphical program as shown in the diagram below with the default values:

The documentation that comes with winfilter explains in detail each of the functions on the GUI and what they are all for, but a brief summary here is worthwhile, I think.

My intention for the program was to be able to explore FIRs, all the different windows functions, and the effects of the different parameters on filter performance. In terms of actually generating some tap coefficients we can actually use, my intention is to do a two phase approach. The first is to select parameters and alter them until we meet the specification required, minimising parameters that cost in terms of logic, such as number of quantisation bits and the number of taps. For this we can plot frequency magnitudes to verify predicted performance. Once satisfied with these results, we can then plot the impulse response. The program will always output the plot values to the file specified in the ‘Configure Output’ box. This would normally be both X and Y values so the glGraph library I’m using can read these and display the graph. A new checkbox has been added, labelled ‘Symmetric Impulse’, which can be checked and, so long as the quantisation value, Q, is non-zero, the output will just be a list of signed hexadecimal numbers of the correct width for the Q value, suitable for configuring the memory in the HDL implementation.

The program does come with an automode design function. This is only available for the Kaiser window function and so this is automatically selected. In this mode, instead of setting? the number of taps and an ‘Alpha’ value and seeing what performance is achieved, one can specify the pass- to stopband transition (the roll off) and a target attenuation in dBs. When executed it will update the Taps and Alpha boxes that achieve that specification, which then informs you the number of taps the HDL implementation needs to be configured for. Note that the Q value for quantisation will also have an effect.

By default, the quantisation value Q is set to 0. This means no quantisation, but actually the program will use double precision floating point values, which is the best we can do. Using Automode when Q is 0 will, indeed achieve the specification configured. However, if you specify a target of -100dB attenuation with a 1KHz roll-off period (other parameters being default), and then set the quantisation to 8 bits, you have no chance of achieving that target, and this parameter must be experimented with as well.

Example Specification

So, for our audio example let’s set some parameters.

• 192KHz sampling rate
• Kaiser windowed-sinc filter
• Pass- to stopband transition < 3KHz
• -60dB attenuation in the stopband

The first parameter we can set from this specification in the quantisation Q. For an attenuation of -60dB we saw from our quantisation plots from before that 12 bits is the smallest we can practically get away with. The number of taps will determine the frequency transition width and we can use the Kaiser window automode to set the 3KHz transition and the target attenuation. This gives 116 taps and an Alpha of 5.6533. If we make the taps 127 (the implementation module’s default parameter setting) and set an Alpha value of 6.0, along with a quantisation of 12 bits, we get a frequency response as shown below:

This appears to meet the brief well. At some of the higher frequencies the lobes do push slightly above -60dB, but the analogue anti-aliasing LPF will be attenuating by those frequencies, bringing them back into line. Now we have our settings we can generate the tap coefficients by selecting ‘Symmetrical Impulse’ in the Configure Output box. The winfilter GUI will look like the following:

If the program is executed with these settings, the generated filter.dat file will now contain 127 signed hexadecimal 13 bits numbers—from the last article we made the HDL LUT values one bit larger that Q so that an equivalent 1.0 can be stored (0x800). The table below shows the first 64 hexadecimal values, remembering we will only store half the table, with the impulse response reflected about the value at index 63:

These numbers have been scaled by the winfilter program to give maximum range, and thus resolution, to the impulse response, with the peak of the impulse response at 2048 (at index 63). This will add a gain to the filter. When winfilter uses floating point values (Q set to 0) the impulse response generated gives unity gain at DC, but when quantised this would reduce the resolution, and so the numbers are rescaled to make the peak the maximum quantised value. This can be set back to unity, if needed, by external means, such as the gain of the anti-aliasing filter, or reduce the filter kernel values to normalise them. This latter, though, will reduce the filter’s response performance.

So, after all this effort, we simply need put these numbers into the TAP lookup table of the logic implementation and we have a filter to our specification. Let’s check if it works.

Performance

To test the FIR implementation, a simple test bench (tb) is provided in the github repository to simulate the module and check the results. This test bench has a set of user definable parameters to configure the test.

The input signal to the UUT can be selected to be either an impulse (IP_FREQ_KHZ is 0), or a sine wave with a frequency defined by IP_FREQ_KHZ, and a peak-to-peak value defined by IP_PK_PK. The system clock rate is defined with CLK_FREQ_HZ, and the input signal sample rate is defined with SMPL_RATE_KHZ. The default values are set up so that the clock rate is? a multiple of the sample rate. I.e., a 96MHz clock and a 192KHz sample rate. The filter’s parameters, SMPL_BITS and TAPS are also exposed for configuration with the test bench parameters. The illustration below shows a simplified block diagram of the test bench setup.

Simulation Results

For the purposes of this document, three tests were performed for the example configuration discussed in the last section: an impulse input, a sinewave at 15KHz in the passband, and a sine wave at 24KHz in the stopband.

The diagram below shows the waveform trace for the impulse response. The signals displayed (from top to bottom) are a clock and reset followed by the smplvalid pulses at the sample rate. The impulse signal, on sample, is next (look to the left), followed by opvalid and then out from the FIR filter. The last signal uses the simulators ‘interpolation’ mode to emulate a DAC and analogue anit-aliasing filter to remove the steps. As can be seen, we do, indeed, get the impulse response at the FIR’s output. This also illustrates the latency? of the filter.

The next plot is for the 15KHz signal, which is at a frequency just before the predicted filter response starts to attenuate. Here, a sine wave is generated and then sampled to input to the filter. As expected, we do see this signal passing through the filter after the latency period.

The last plot does the same thing as the previous test, but the input signal frequency is now 24KHz, the lowest frequency in the predicted stopband. Here we see that the signal is filtered and does not pass through. There is some initial response after the filter latency, and this is due to the non-linearity of the transition from no input samples to a sine wave input temporarily generating frequencies in the passband, but it soon settles down to attenuating the input signal.

Conclusions

In this article, we built on the foundations of the previous article, discovering ?convolution, and with an HDL implementation ready to program with some numbers. In this article was discussed how to generate those numbers for a lowpass filter, from the sinc function being the digital fourier transform of the rectangular frequency response to give an infinite impulse response. This was turned into a finite impulse response using ‘window’ functions of varying degrees of complexity and performance. To make practical for the Verilog implementation the finite impulse response was then quantised, and the limitations explored by this step.

Having? generated a practical set of coefficients for a low pass filter, we then looked at how to manipulate and combine the tap coefficient to produce high pass and band pass filters.

With all this knowledge a practical example from the world of digital audio was specified and a simulation constructed to test the FIR implementation with some filter tap coefficients to match that specification. The results did indeed confirm the filter’s response to impulses, passband frequencies and stopband frequencies to be as predicted.

The convolution/FIR logic can be seen, I hope to be fairly simple, with the complexity in the calculation of the tap coefficients. Although more complex, the mathematics of this is not, I’d like to think, at the realm of professors of pure mathematicians, and at least means none of it ends up in the logic implementation and can be done offline. The resultant FIR, though, has impressive characteristics for such simple functionality, and is useful in a wide range of applications where stable and linear filtering is required.

More Information

For a very large section of my career, the foundation for the DSP knowledge I have has come from a book “The Scientist and Engineer's Guide to Digital Signal Processing”, By Steven W. Smith. I cannot recommend this book highly enough for anyone starting out in digital signal processing. It’s writing style and clarity is something that I have tried to emulate in my own documentation and these articles. Even, better, you can get hold it for free, on Steven Smith’s website at www.dspguide.com. This book includes way more information and detail that I can fit in a few articles, but includes much more than FIR filters. I still refer to this text (I actually bought a hard copy, despite it being free). As far as implementation goes, it focuses more on digital signal processors, where I aim towards logic implementation, but is very much relevant for both disciplines. So, if you enjoyed these articles, and want to explore DSP more, check out this book.