Digital Filters: Understanding the Mechanics

Digital Filters: Understanding the Mechanics

16 Comments

Filters are perhaps the most widespread devices in signal processing. Applications cover analysis, signal enhancement, noise reduction, and much more. The goal of this article is to provide some basic insight to the inner workings of digital filters, and in particular to show the relation between filter coefficients and filter characteristics. To this end we will need some math, however I’ll try to restrict it to the absolutely necessary.

A Simple Filter

A digital filter acts on a digital signal, which is given as a sequence of incoming samples. The filter will modify the stream of samples and send out a resulting stream. Consider the following first order recursive filter:

filter

It operates on the actual input sample in and the previous input and output samples in1 and out1, respectively. These are multiplied by coefficients b0, b1, and a1, respectively, to yield the current output sample out. Then the actual samples in and out are assigned to the previous samples in1 and out1 for the next iteration.

Let’s have a look at how our filter modifies some common input signals. In the following graphs, the light blue curve represents the input while the filtered signal is shown in red.

Triangle                                                                                            White Noise
filteredTriangle         filteredNoise

Sine                                                                                                    Square
filteredSine        filteredSquare

In general, we see that the amplitude is reduced by the filter, and the shape of the waveform gets distorted – except for the sine wave, where the original shape is retained (apart from attenuation and a delay or, more precisely, a phase shift in the filtered signal). For this reason, sine waves are particularly suited to test or characterize filters.

Fortunately, any incoming signal can be decomposed into a combination of sine waves with different frequencies, amplitudes, and phases. Once we know how the filter acts on each of these sine waves, we can recombine the results to get the resulting filter output for any input signal.

Obviously, the filter behavior is determined by the coefficients b0, b1, and a1, in the above first order filter example. Different choices will result in a lowpass, hipass, low shelf, high shelf, or allpass filter. (Other types like bandpass, notch or peaking filter would require a slightly more complex filter, e.g. a biquad.) The relation between filter coefficients and filter characteristics is not obvious even for this simple filter, therefore we will try to gain insight in the following. Ultimately we want to be able to design a filter according to a given specification.

Figuring out the Filter

Let us begin by writing down the filter recursion equation. We will use letters x and y instead of words in and out for the input and output samples, respectively. At the nth iteration step, the current input sample will be denoted by xn while the previous sample will be denoted by xn-1. Likewise for the output samples. Using this notation, the first code line of our filter reads

FilterEqn.                                                                     (1)

This equation can be solved by direct iteration: each new output value is computed from the current and previous input and output values – but that’s exactly what our filter code does! There is a better way, a trick that allows us to solve for the entire output signal in one single step! Read on.

The Z-Transform . . .

The so-called Z-transform is given by multiplying each sample with the nth power of some number z and taking the sum over all nonzero samples:

zTransform.                                                                                           (2)

A similar transform may be applied to the output samples yn. The expression on the right hand side of equation (2) yields different values for different z, hence we obtain a function X(z).

Applying the Z-transform to the filter equation (1) is straight forward for the xn and yn terms; however, the other terms involving the previous samples xn-1 and yn-1 require special consideration. It can be shown that shifting the index towards the previous sample amounts to multiplying by z-1 in the transformed domain. That’s why z-1 is sometimes called the unit delay operator.

Putting it all together, we get the following transformed filter equation

transformedFilterEqn.                                          (3)

Equation (3) is much simpler that the original equation (1) because it establishes a relation between the input X and the output Y for the same z. Hence, with a modest amount of algebra, we can solve for Y to obtain the following simple expression:

transferEqn.                                                                                           (4)

In equation (4), H(z) denotes the transfer function and is given by

transferFunction.                                                                                         (5)

To summarize the result obtained so far: in the transformed domain, the action of the filter is a simple multiplication of the input signal X(z) by the transfer function H(z). The transfer function is a rational function of z.

In order to actually make use of this result, we need to elaborate a bit on the meaning of the z parameter in the following.

. . . Demystified!

There is a close relation between the Z-transform in equation (2) and the Fourier transform – so close that it is fair to call the Z-transform a disguised Fourier transform! To see this, substitute for z the expression ejω, where ω = 2πf /sample rate is the (dimensionless) angular frequency (in radians) and j is the imaginary unit (yes, we will use complex numbers here). Then equation (2) becomes

FourierTransform,                                                                                      (6)

where the right hand side is an ordinary Fourier series. Thus the left hand side of equation (6), X(ω), is simply the frequency spectrum of the input sequence xn. The same interpretation applies to all other functions of z.

With this insight we can think of equation (4) in terms of frequency, if we substitute ejω for z. The frequency range from DC (ω = 0 or, equivalently, z = 1) to the Nyquist limit (ω = π or z = −1) maps, on the z-plane, onto the upper half of the unit circle. Formally we also include negative frequencies to complete the unit circle. The situation is depicted in the figure below.

unitCircle

In a nutshell, the z-transform is a disguised Fourier transform, and the z variable is an encoded frequency variable.

Designing Filters

We are now ready to design a filter. Our first example will be a DC blocker, i.e. a high pass filter which will be transparent for audio frequencies and will reject subsonic frequencies. Blocking DC means that the transfer function H(z) is zero for z = 1. From equation (5) we see that this requirement implies b1 = −b0, i.e. the two coefficients have equal magnitude but opposite sign. At the other end of the spectrum, the filter is supposed to be transparent, so H(z) = 1 for z = −1. Putting this into equation (5), we get

b0 = −b1 = (1 − a1)/2.

The remaining coefficient a1 determines the cutoff frequency ωC which separates the pass band from the stop band. You may derive the relation from equation (5) by setting |H(exp(jωC))|2 = 1/2; the result is

cutoffRelation                                                                                            (7)

For a 20 Hz cutoff frequency and assuming a sample rate of 44100 Hz we get

ωC = 2π∙20/44100 = 0.00285,

and from there the DC blocker filter coefficients:

a1 = −0.99715, b0 = 0.99858, and b1 = −0.99858.

As a second example, we will design a smoothing or low pass filter. We will want unity gain for DC, H(z) = 1 for z = 1, and full rejection at Nyquist, H(z) = 0 for z = −1. Putting this into equation (5) yields

b0 = b1 = (1 + a1)/2.

We choose a 740 Hz cutoff frequency and assume a 44100 Hz sample rate, so

ωC = 2π∙740/44100 = 0.1054.

Using this value in equation (7) yields the low pass filter coefficients:

a1 = −0.9, b0 = 0.05, and b1 = 0.05.

If you scroll up, you’ll notice that these are the same coefficients as those used in our code example.

So now you know how first order filters work. Wanna try biquads next?

;)

0 0 0 0 0

About the author:

16 Comments

  1. kohugaly
    kohugaly  - December 10, 2014 - 10:49 pm

    Thank you Martin!!! Finally, I understand how to work with Z-transform. I was trying to wrap my head around it for so long…

    • kohugaly
      kohugaly  - December 10, 2014 - 11:19 pm

      One thing is not clear to me though… how exactly did you got the equation (6)? (by the way there is equation (5) twice ;-), a typo?

  2. martinvicanek  - December 11, 2014 - 7:36 am

    Thanks, KG, I’ve sorted out the equation numbering, Eq. (7) may be derived from equation (5) by setting |H(exp(jω_C))|^2 = 1/2. Hint: use exp(jω_C) = cos(ω_C) + j*sin(ω_C).

    • kohugaly
      kohugaly  - December 11, 2014 - 10:27 am

      Aha… now I get it… |H(exp(jω_C))|^2 = 1/2 basically means that the gain at cutoff frequency should be sqrt(0.5)=0.702=-3dB. The thing is squared because we want to set only the magnitude at that frequency, phase being whatever it comes out. The value of H(z) is an imaginary number = frequency-domain representation, therefore it has magnitude (length from [0,0]) and phase (angle from [0,0]).

      I’m working on an article, that will hopefully demystify the frequency domain and discrete Fourier Transform, showing what frequency-domain representation actually means. So guys, if you are not familiar with why Martin used e^-i*2pi*f then stay tuned – you will get your answer soon with less math and more schematics to show, what’s happening.

    • drno  - December 17, 2015 - 1:57 pm

      i don‘t understad the equation 7. where come from those sin and cos?

      • martinvicanek  - December 17, 2015 - 11:26 pm

        In eq.(5) substitute exp(j*w_c) = cos(w_c) + j*sin(w_c). Then takte |H(z)|^2 and equate to 1/2.

        • drno  - December 20, 2015 - 2:01 am

          sorry, i’m litle slow with the mathematic,
          you cam explain me step by step thanks

          • martinvicanek  - December 20, 2015 - 2:47 pm
            /

            Start from eq. (5) with the expressions for b0 and b1:

            H(z) = ((1 – a)/2)*(1 – c + is)/(1 + ca – isa)

            where a is short for a1, c = cos(omega_C) and s = sin(omega_C). Take the moduls squared:

            |H(z)|^2 = ((1 – a)^2/4)*(1 – 2c + c^2 + s^2)/(1 + 2ca + c^2a^2 + s^2a^2)

            = (1 – 2a + a^2)/2*(1 – c)/(1 + 2ca + a^2)

            Equating this expression to 1/2 leads to a quadratic equation for a:

            a^2 + (2/c)*a + 1 = 0,

            which has two solutions

            a1 = -(1 + s)/c, a2 = -c/(1 + s).

            The first solution results in a pole outside the unit circle (unstable filter), hence we discard it and keep the second, which is the one in the text.

  3. Exo
    Exo  - December 11, 2014 - 4:58 pm

    Great Article Martin thanks! I feel like I understand a bit more. I am going to have to get some of the things I am working on finished before I can start experimenting with filters, hopefully soon!

  4. tulamide
    tulamide  - December 20, 2014 - 8:22 pm

    Thanks, Martin. I’m still trying to get at least the basics, but it’s so good to have such sources here!

    What I don’t quite understand is right from the beginning. It seems the filter adds 5% of the last input to 5% of the current input and then add 90% of the last output. But why is a1 then -0.9 and the equation …- a1 * out1? Wouldn’t it make more sense to set a1 to 0.9 and …+ a1 * out1? Also, why out1 and not just the current out? I don’t expect an answer, I know it can bug. But if there are easy answers and you don’t mind, feel free to answer xD

    • kohugaly
      kohugaly  - December 20, 2014 - 8:52 pm

      The minus sign for A coefficients phenomenon arises because the equation actually should look like this:
      a0*out0+a1*out1=b0*in0+b1*in1
      Obviously to calculate the value of “out0″ you need to put everything on the right-side and divide by a0. That’s why all coefficients A are with minus sign. The value of the actual A coefficients may end up positive or negative – in this case negative -> you are subtracting negative number = adding.

      “out1″ isn’t simply “out” for educational reasons. It might not be that obvious that out keeps it’s value, so martin emphases this fact by introducing out1 variable (this becomes necessary for higher order filters, like biquad, to have out1,out2,out3,..)

      • tulamide
        tulamide  - December 20, 2014 - 10:26 pm

        “The value of the actual A coefficients may end up positive or negative – in this case negative -> you are subtracting negative number = adding.”

        But that’s the point. I was asking, if I change anything in the filter’s behavior, if I use a1 = 0.9 and …+ a1 * out1 instead of how it is now? Both ways should lead to the same (mathematical) result, or not?

        Thanks for the out reminder! I wasn’t aware that out keeps the value just as in! Now it makes sense.

        • kohugaly
          kohugaly  - December 20, 2014 - 10:50 pm

          Yes, in this case it works. However, if you make the filter parametric (like the input is a cutoff frequency, gain etc.) and the coefficients are actually being calculated inside the code, you have to keep the scheme of minus before A coefficients.
          In some literature and sources, you may find version that add instead of subtract (so the A coefficients are actually inverted). It’s merely a convention. I even found some sources that used A to name coefficient multiplied by input and B for those multiplied by outputs (the opposite way then in this article).

          • tulamide
            tulamide  - December 21, 2014 - 5:47 pm
            /

            Thanks once again. It’s crazy how, after reading the answers, I think “what a obvious thing, how could you not have seen it?” :D

  5. great blog  - January 21, 2015 - 11:01 pm

    magnificent issues altogether, you just gained a logo new|a new} reader. What might you recommend in regards to your post that you just made some days in the past? Any certain?

Add Comment Register



Leave a comment

You must be logged in to post a comment.

Back to Top